地下水水源地保護區(qū)劃分方法研究
——以大汶口盆地為例

王 昕 翌，王 敏，馮 建 國，高 宗 軍

(山東科技大學 地球科學與工程學院，山東 青島 266510)

地下水是我國重要的供水資源，地下水水質(zhì)的優(yōu)劣直接關系到人類的飲水安全和健康狀況。目前，我國居民飲水安全狀況不容樂觀，飲用地下水水源污染已成為一個較突出的社會問題，地下水保護的一個直接有效的途徑就是建立地下水水源地保護區(qū)[1-4]。國內(nèi)外已有許多關于地下水水源地保護區(qū)劃分方面的研究與討論，如田文英[5]通過解析模型法對遼寧省柳河平原區(qū)水源地進行了保護區(qū)劃分；王麗娟等[6]選用穩(wěn)定流和非穩(wěn)定流模型法進行了保護區(qū)劃分。但是，地下水水源地的復雜性導致運用單一的方法進行區(qū)劃劃分會存在劃分結(jié)果與實際不符、劃分精度不高等問題。因此，學者們對區(qū)劃方法做出了許多改進：如趙紅梅等[7]選用公式法和數(shù)值模擬法2種方法對成都平原某水源地進行保護區(qū)劃分，并比較2種方法的區(qū)劃結(jié)果，認為在實際應用中2種方法相結(jié)合所得結(jié)果更為合理；江廣長等[8]對多種保護區(qū)劃分方法進行對比分析，認為各類劃分方法的不足之處是未能綜合考慮水源地內(nèi)各項影響因素，區(qū)劃方法與地下水水源地間的聯(lián)系不甚緊密，并提出在數(shù)值模擬法中引入不確定分析的觀點。

空間疊加技術在水文地質(zhì)研究中的廣泛應用，也為改進區(qū)劃方法提供了思路。如Dawoud[9]通過空間疊加手段建立了覆蓋埃及全國的地下水質(zhì)量監(jiān)測網(wǎng)，將多組包含了不同含水層類型和水化學特征的空間數(shù)據(jù)圖層疊加，實現(xiàn)了如人口密度、水文地質(zhì)條件、鹽漬化風險等眾多影響因素的信息整合。將空間疊加技術應用到區(qū)劃工作中，有助于提高劃分精度，使區(qū)劃結(jié)果能夠全面地反映出水源地的特征。

結(jié)合前人的研究，本文以大汶口盆地為例，首先選用數(shù)值模擬法對盆地內(nèi)水源地進行保護區(qū)劃分，在此基礎上，應用空間疊加分析，將數(shù)值模擬法與地下水污染風險評價、防污性能評價結(jié)合，修正數(shù)值模擬法的區(qū)劃結(jié)果，并將2種方法所得結(jié)果進行對比討論，以期完善保護區(qū)劃分方法，為地下水安全保護提供依據(jù)。

1 大汶口盆地概況

大汶口盆地位于山東省西南部，地跨泰安、新泰、肥城等縣市，城鎮(zhèn)眾多，交通便利，農(nóng)業(yè)發(fā)達，礦產(chǎn)豐富。盆地內(nèi)地表水系發(fā)育，區(qū)域性河流為汶河水系，位于盆地中南部，主要支流有大汶河和漕河。整體地勢東高西低，南高北低，總體形態(tài)為單斷凹陷，是一個較完整的水文地質(zhì)單元。

盆地內(nèi)地下水類型有巖漿巖風化裂隙水、碎屑巖類孔隙裂隙水、松散巖類孔隙水和碳酸鹽巖裂隙巖溶水4種類型。地下水主要依靠大氣降水補給，汶河的側(cè)向補給和滲漏對大汶口至西界一帶的孔隙水及巖溶水起著重要作用。由于盆地受構(gòu)造控制，各時代沉積物差異很大，地下水的開發(fā)利用受地層及富水性控制明顯，地下水開采以巖溶水為主，其次為孔隙水和裂隙水。

盆地內(nèi)有東武、樓德2處地下水水源地，是當?shù)鼐用竦闹饕┧吹亍?處水源地均以開采地下巖溶水為主，水位動態(tài)主要受大氣降水、汶河徑流量以及地下水開采量的影響，排水途徑均為補給上覆第四系孔隙水，側(cè)向徑流排泄及人工開采，二者間的水力聯(lián)系十分密切。其中東武地下水水源地分布在東武駕莊-大汶口鎮(zhèn)一帶，開采井深度為100～200 m，單井涌水量一般在5 000 m3/d以上，地下水開采量2.2萬m3/d；樓德地下水水源地分布在樓德-甘露村-西杜村一帶，開采井深度為130～300 m，單井涌水量一般在3 000 m3/d 以上，地下水開采量為2.0萬m3/d。自1989年以來，東武、樓德2處水源地長期處于開采狀態(tài)，2006年12月至2018年12月，東武地下水水源地平均水位下降了41.58 m，樓德地下水水源地平均水位下降了38.97 m，巖溶水水位呈逐漸下降趨勢。同時工業(yè)廢水、農(nóng)藥化肥等污染源不斷增多，滲入地下導致地下水污染，地下水水質(zhì)急劇惡化。根據(jù)1996年水質(zhì)監(jiān)測資料顯示，東武、樓德水源地水化學類型均以HCO3—Ca型水為主，經(jīng)2018年調(diào)查，東武水源地水化學類型以HCO3·SO4·Cl—Ca型水為主，可溶性固體平均為1 291 mg/L，總硬度平均為889 mg/L，SO42-濃度平均為314 mg/L；樓德地下水水源地水化學類型為HCO3·SO4—Ca型水，總硬度平均為330 mg/L，SO42-濃度平均為30 mg/L，HCO3-濃度為230 mg/L。與開采之初相比，2處水源地地下水各項指標均有所增加，SO42-濃度明顯升高，水質(zhì)狀況以Ⅲ類較差水和Ⅳ類極差水為主，地下水環(huán)境保護刻不容緩。

2 地下水水源地保護區(qū)劃分方法

2.1 劃分方法與思路

地下水水源地并非獨立存在的，開展地下水水源地保護區(qū)劃分工作，應在全面了解整個水文地質(zhì)單元特征的前提下，以地質(zhì)單元內(nèi)現(xiàn)有水源地為核心進行區(qū)劃，將保護區(qū)的范圍涉及整個水文地質(zhì)單元。傳統(tǒng)的數(shù)值模擬法雖然能夠較客觀地刻畫含水層結(jié)構(gòu)及水文地質(zhì)條件，但依然存在參數(shù)的不確定性、模型概化過程中易出現(xiàn)誤差等問題[8]。因此，本文首先應用傳統(tǒng)的數(shù)值模擬法，在建立研究區(qū)地下水水流模型和水質(zhì)模型的基礎上，依據(jù)各級保護區(qū)對污染溶質(zhì)運移時間的不同要求，來初步劃定不同級別保護區(qū)的邊界范圍。在此基礎上，將地下水防污性能評價和污染風險評價引入?yún)^(qū)劃方法中，把研究區(qū)的地質(zhì)條件、人類活動、污染源分布情況等因素納入?yún)^(qū)劃標準中，綜合評估整個水文地質(zhì)單元的地下水防污能力和污染風險，然后利用空間疊加技術，整合數(shù)值模擬法和2項評價所得的空間數(shù)據(jù)圖層，同時結(jié)合研究區(qū)自然地理和水文地質(zhì)條件對保護區(qū)邊界進行調(diào)整，得出最終的地下水水源地保護區(qū)范圍(見圖1)。

圖1 保護區(qū)劃分方法流程示意Fig.1 Structure diagram of protection zone division method

2.2 數(shù)值模擬法—保護區(qū)初步劃分

應用數(shù)值模擬法劃分地下水水源地保護區(qū)范圍的原理，是根據(jù)研究區(qū)的含水層結(jié)構(gòu)、地質(zhì)構(gòu)造等基礎數(shù)據(jù)，建立研究區(qū)的水文地質(zhì)概念模型，在此基礎上，利用模擬軟件如GMS(Groundwater Modeling System)、Visual Modflow等來模擬研究區(qū)的地下水流場。根據(jù)研究區(qū)內(nèi)污染源的空間分布情況和水質(zhì)資料，進行地下水溶質(zhì)運移模擬，對地下水水源地周邊的污染物入滲點位置進行預測，在水源地周邊不同位置設置入滲點，其中那些污染溶質(zhì)運移100 d和1 000 d后剛好到達水源地的入滲點，分別為水源地一級和二級保護區(qū)的邊界點[10-11]。

2.3 地下水防污性能評價

目前國內(nèi)外一般采用DRASTIC法進行地下水防污性能評價。DRASTIC法通常選取地下水埋深D、地下水的凈補給量R、含水層介質(zhì)A、土壤類別S、地形坡度T、包氣帶介質(zhì)I和含水層水力傳導系數(shù)C這7項影響因子，通過對指標進行量化評分來代表各項因子的數(shù)值范圍或類別(見表1)，并給予各項因子不同的權重參數(shù)[12]，DRASTIC地下水防污性能DI值由公式

表1 評價因子評分體系Tab.1 Evaluation factor scoring system

(1)確定：

DRASTIC=5D+4R+3A+2S+1T+5I+3C

(1)

7項因子的權重賦值可根據(jù)研究區(qū)實際情況進行調(diào)整，使評價結(jié)果更為合理。在一般條件下，地下水防污性能指數(shù)為23～226，根據(jù)計算后所得DI值，將地下水防污性能指數(shù)DI共分為5個等級(見表2)，DI值越高，則地下水防污能力越差。

表2 地下水防污性能分級標準Tab.2 Classification standards of groundwater anti-fouling performance

2.4 地下水污染風險評價

地下水污染風險評價是在污染載荷風險(P)評價、含水層固有脆弱性(U)評價、污染危害性(V)評價的基礎上進行的。污染源荷載風險(P)的評價指標為：研究區(qū)內(nèi)污染源種類K、污染物產(chǎn)生量Q、污染物釋放可能性L和調(diào)查點與污染源間的距離D，并依據(jù)研究區(qū)實際情況對4項評價指標分別打分[13-14]，按公式(2)計算單個污染載荷P′的值。

P′=KQLD

(2)

將單個污染載荷P′的計算結(jié)果按照表3的分級標準重新分類，即為污染載荷風險(P)評價的結(jié)果。

表3 污染載荷風險評價結(jié)果重新分級標準Tab.3 Pollution load risk assessment results reclassification criteria

含水層脆弱性的高低與地下水防污能力的強弱密不可分，地下水防污能力強，則表示脆弱性低；地下水防污能力弱，即脆弱性高[15]。故含水層固有脆弱性評價(U)將地下水防污性能評價的結(jié)果進行重新分級并賦值，賦值標準見表4。

污染危害性(V)一般根據(jù)研究區(qū)內(nèi)地下水的用途來進行分級，因此本文據(jù)此將污染危害性分為3級(見表5)。最終按照表6的評價標準得到污染風險評價結(jié)果R，其中“0”代表低污染風險；“1”表示中等風險；“2”表示高污染風險[14]。

表4 含水層固有脆弱性重新分級標準Tab.4 Reclassification criteria for the inherent vulnerability of aquifers

表5 污染危害性分級標準Tab.5 Pollution hazard classification standards

表6 地下水污染風險評價標準Tab.6 Risk assessment criteria for groundwater pollution

2.5 空間疊加分析

通過對不同圖層的疊加分析，可以整合多層數(shù)據(jù)間的聯(lián)系與變化特征[9，16]。應用GIS軟件如MapGIS、ArcGIS平臺上的空間分析功能，首先進行空間疊加求交運算，分離出地下水防污能力低且污染風險高的區(qū)域和地下水防污能力低且污染風險中等、地下水防污能力中等且污染風險高的區(qū)域；再通過空間疊加求并運算，將上述區(qū)域分別納入數(shù)值模擬法擬定的一、二級保護區(qū)范圍內(nèi)。空間疊加分析分級標準列于表7。

3 大汶口盆地地下水水源地保護區(qū)劃分

本文將大汶口盆地整個水文地質(zhì)單元作為研究區(qū)，以盆地內(nèi)東武、樓德2處地下水水源地為核心，根據(jù)2018年調(diào)查獲得的45個承壓水水質(zhì)監(jiān)測點及盆地基礎數(shù)據(jù)，利用數(shù)值模擬法和基于數(shù)值模擬法的空間疊加方法對大汶口盆地進行保護區(qū)劃分，并對2種方法所得結(jié)果進行對比討論。

表7 空間疊加分析分級標準Tab.7 Spatial superposition analysis division criteria

3.1 大汶口盆地數(shù)值模擬結(jié)果及保護區(qū)初步劃分

根據(jù)大汶口盆地水文地質(zhì)條件，盆地內(nèi)地下水開采以巖溶水為主，且?guī)r溶含水層與上覆孔隙含水層間有穩(wěn)定分布的巨厚第三系泥頁巖阻隔，故本次模擬以巖溶含水層為目的含水層，將其概化為非均質(zhì)各向同性的承壓含水層，水流運動為平面二維流。模擬區(qū)包括一類水頭邊界和二類流量邊界，其中盆地北側(cè)為多組斷裂視為隔水邊界，大汶河流域處理為水頭邊界。模擬時僅考慮對流、彌散作用對于溶質(zhì)濃度變化的影響，依據(jù)盆地經(jīng)驗數(shù)據(jù)，孔隙度取值為0.07，彌散系數(shù)取值為15 m2/d。

數(shù)值模擬模型包括水流模型和水質(zhì)模型2個部分，模擬選用GMS6.0軟件中的MT3D模塊，在運用ModFlow模塊模擬出地下水流場的基礎上運用MT3D模塊進行地下水溶質(zhì)運移模擬。模擬區(qū)內(nèi)地下水流特征用如下數(shù)學模型表示：

(3)

式中：K為巖溶含水層滲透系數(shù)，m/d；W為含水層的源匯項，m2/d；Ss為承壓含水層的貯水率；τ為流量邊界；h0為初始水位。

模擬區(qū)地下水溶質(zhì)運移用如下數(shù)學模型表示：

(4)

結(jié)合大汶口盆地內(nèi)污染源的空間分布情況，分別在東武和樓德水源地周圍不同位置設置污染物入滲點進行試算，其中那些污染質(zhì)運移100 d和1 000 d后剛好分別到達東武、樓德地下水水源地的入滲點，分別為其一級和二級保護區(qū)的邊界點(見圖2)。將這些邊界點圈閉，即為水源地一級、二級保護區(qū)的理論范圍(見圖3)。其中，東武水源地一級保護區(qū)面積為17.2 km2，樓德水源地一級保護區(qū)面積為15.3 km2，模擬結(jié)果顯示2處水源地的二級保護區(qū)范圍互相交疊，故將二者的二級保護區(qū)范圍劃定為一個整體，二級保護區(qū)面積為456.3 km2。

圖2 溶質(zhì)運移模擬Fig.2 Solute transport simulation diagram

圖3 數(shù)值模擬法初步劃定保護區(qū)范圍Fig.3 The scope of protected areas is preliminarily determined by numerical simulation method

3.2 大汶口盆地地下水防污性能評價

由于大汶口盆地內(nèi)農(nóng)田大量噴灑農(nóng)藥，所以為貼合盆地實際情況，對地下水防污性能評價的權重稍作調(diào)整，適當增加土壤類別S、地形坡度T和包氣帶介質(zhì)I這3項影響因子的權重，調(diào)整后的地下水防污性能DI值由公式(5)確定：

DRASTIC=5D+4R+3A+5S+3T+4I+3C

(5)

根據(jù)盆地內(nèi)各巖溶水點的調(diào)查數(shù)據(jù)，結(jié)合表1的評分體系和表2的分級標準，將大汶口盆地地下水防污性能分為四級(見表8)。防污性能高和非常高的地區(qū)呈帶狀分布在盆地北部、東北部和南部山區(qū)，富水性小，第三系蓋層較厚，地下水防污能力強；防污性能低區(qū)主要分布在大汶口鎮(zhèn)、樓德鎮(zhèn)周邊地區(qū)以及大汶河富水性大的地段，河流對地下水的補給量大，局部地區(qū)第三系蓋層缺失，地下水易受地表水的影響；防污性能中等的地區(qū)主要分布在防污性能低區(qū)的外圍。大汶口盆地地下水防污性能評價空間分布如圖4所示。

表8 地下水防污性能評價綜合成果Tab.8 Comprehensive evaluation results of groundwater anti-fouling performance

圖4 大汶口盆地地下水防污性能評價空間分布Fig.4 Spatial distribution of groundwater anti-fouling performance evaluation

3.3 大汶口盆地地下水污染風險評價

綜合大汶口盆地內(nèi)各調(diào)查點污染載荷風險評價(P)、含水層固有脆弱性評價(U)、污染危害性評價(V)的最終評分，對應表6的評價標準，將大汶口盆地地下水污染風險分為三級(見表9)。污染風險高的地區(qū)主要分布在汶河北岸、大汶口鎮(zhèn)和樓德鎮(zhèn)周邊地區(qū)，污染渠道較多且人口眾多，地下水開采量大且多用于飲用，另外汶河南岸第三系缺失地區(qū)，地下水受污染河流直接補給，地下水受污染可能性極大；污染風險中等的地區(qū)主要分布在高污染風險地區(qū)的外圍；污染風險低的地區(qū)主要分布在盆地北部、東部和南部的基巖山區(qū)，人類活動稀少，對地下水的開采利用較少，并且遠離污染源。大汶口盆地污染風險評價空間分布如圖5所示。

表9 地下水污染風險評價綜合成果Tab.9 Comprehensive results of groundwater pollution risk assessment

圖5 大汶口盆地地下水污染風險評價空間分布Fig.3 Spatial distribution of groundwater pollution risk assessment

3.4 空間疊加及保護區(qū)范圍確定

在MapGIS軟件中的空間分析平臺上，將地下水防污性能評價結(jié)果(見圖4)、地下水污染風險評價結(jié)果(見圖5)和數(shù)值模擬法劃定的保護區(qū)范圍(見圖3)進行空間疊加(見圖6)。依據(jù)表7的疊加分級標準，把大汶口盆地內(nèi)地下水防污能力低且污染風險高的區(qū)域加入一級保護區(qū)范圍內(nèi)；把地下水防污能力中等且污染風險高、地下水防污能力低且污染風險中等的區(qū)域加入二級保護區(qū)范圍內(nèi)。

圖6 空間疊加過程示意Fig.6 Schematic diagram of spatial superposition process

根據(jù)疊加整合結(jié)果，雖然盆地西部汶河下游區(qū)域，地下水防污能力弱且污染風險高，但距離東武、樓德水源地極遠，是地下水的排泄區(qū)，故不納入地下水水源地一、二級保護區(qū)范圍內(nèi)。最后參照盆地的自然環(huán)境和水文地質(zhì)條件對保護區(qū)邊界進行適當調(diào)整：盡量利用盆地內(nèi)的交通道路、隔水斷層等界線作為保護區(qū)邊界，結(jié)合盆地的地形地貌對保護區(qū)邊界進行適當增加或縮減，準保護區(qū)范圍擴展至整個盆地，即以水文地質(zhì)單元的地表分水嶺為邊界。疊加后東武水源地一級保護區(qū)面積為33.6 km2，位于盆地中部；樓德水源地一級保護區(qū)面積為30.4 km2，位于盆地東南部；二級保護區(qū)面積為569.5 km2(見圖7)。

圖7 大汶口盆地地下水水源地保護區(qū)Fig.7 Groundwater source protection area of Dawenkou basin

與數(shù)值模擬法的劃分結(jié)果相比，東武水源地一級保護區(qū)邊界向東北方向有延伸，增加面積為13.2 km2；樓德水源地一級保護區(qū)邊界向東南方向有延伸，增加面積為15.1 km2；二級保護區(qū)向西側(cè)、南側(cè)以及東南側(cè)延伸，面積增加113.2 km2，一級、二級保護區(qū)邊界均向地下水補給源方向擴展。根據(jù)結(jié)果顯示，經(jīng)過空間疊加后得到的大汶口盆地地下水水源地保護區(qū)，包括了盆地內(nèi)地下水補給區(qū)中污染風險高和防污能力弱的區(qū)域，得到的結(jié)果更為貼近大汶口盆地實際情況。

2種方法的劃分結(jié)果存在差異，主要原因是傳統(tǒng)的數(shù)值模擬法雖然考慮了地下水徑流方向等因素的影響，但由于數(shù)值模擬過程中存在的種種誤差，使得劃分出的保護區(qū)與實際情況存在出入。而在對盆地進行地下水污染風險和防污性能評價時，能夠充分利用滲透系數(shù)、包氣帶厚度等地質(zhì)點量化信息，同時亦將盆地內(nèi)污染源的空間分布及地下水用途等因素納入了區(qū)劃標準，客觀反映了盆地的水文地質(zhì)條件和地下水環(huán)境狀況，最終通過空間疊加手段，將綜合了多項指標的地下水污染風險和防污性能評價結(jié)果整合入數(shù)值模擬法中，優(yōu)化數(shù)值模擬法的區(qū)劃結(jié)果，彌補了數(shù)值模擬法的不足。

基于數(shù)值模擬的空間疊加法的優(yōu)點是思路清晰，適用性強，可以根據(jù)研究區(qū)的實際情況來選擇評價方法和需要疊加整合的要素，結(jié)果貼合實際；本次針對實際水源地的區(qū)劃過程也反映出該方法工作量較大且需要足夠多的基礎資料支撐，需要以研究區(qū)準確的水文地質(zhì)信息及充足的調(diào)查資料為依托。本次區(qū)劃由于盆地資料限制，僅考慮了地下水受污染的可能性和防污性能等要素，而如地下水水質(zhì)、水量的動態(tài)變化等因素未能綜合處理，在后期工作中，可以將更多可能對區(qū)劃結(jié)果造成影響的因素納入考量，優(yōu)化疊加分級標準，進一步完善數(shù)字化的綜合評價體系。該方法適用于水文地質(zhì)背景復雜的研究區(qū)，尤其對于以開采巖溶水為主的大型水源地效果更佳，在實際應用中，可針對不同研究區(qū)的具體情況來調(diào)整評價手段和打分體系，并與其他區(qū)劃方法和評價手段相結(jié)合，使保護區(qū)劃分結(jié)果更加準確合理。

4 結(jié) 論

(1) 在數(shù)值模擬法所得保護區(qū)范圍基礎上，通過空間疊加技術整合地下水防污性能和污染風險評價的結(jié)果，得到大汶口盆地東武地下水水源地一級保護區(qū)面積為33.6 km2，樓德地下水水源地一級保護區(qū)面積為30.4 km2，二級保護區(qū)面積為569.5 km2，準保護區(qū)面積為1 490.0 km2。區(qū)劃結(jié)果更為貼合大汶口盆地實際情況，證明基于數(shù)值模擬的空間疊加法具有可行性。

(2) 地下水防污性能和污染風險評價可以較全面地綜合研究區(qū)的水文地質(zhì)特征、人類活動及污染源的空間分布等環(huán)境要素，而空間疊加技術有效地將2項評價結(jié)果與數(shù)值模擬法的區(qū)劃結(jié)果結(jié)合，修正傳統(tǒng)數(shù)值模擬法存在的誤差，使區(qū)劃結(jié)果能更加全面地反映出水源地的特征，為保護區(qū)內(nèi)污染防治、土地利用、設施建設等做出相關規(guī)定提供依據(jù)，可更好地實現(xiàn)地下水資源的保護，提高地下水資源的利用效率。