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

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

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

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

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

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

1 大汶口盆地概況

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

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

盆地內有東武、樓德2處地下水水源地，是當地居民的主要供水水源地。2處水源地均以開采地下巖溶水為主，水位動態主要受大氣降水、汶河徑流量以及地下水開采量的影響，排水途徑均為補給上覆第四系孔隙水，側向徑流排泄及人工開采，二者間的水力聯系十分密切。其中東武地下水水源地分布在東武駕莊-大汶口鎮一帶，開采井深度為100～200 m，單井涌水量一般在5 000 m3/d以上，地下水開采量2.2萬m3/d；樓德地下水水源地分布在樓德-甘露村-西杜村一帶，開采井深度為130～300 m，單井涌水量一般在3 000 m3/d 以上，地下水開采量為2.0萬m3/d。自1989年以來，東武、樓德2處水源地長期處于開采狀態，2006年12月至2018年12月，東武地下水水源地平均水位下降了41.58 m，樓德地下水水源地平均水位下降了38.97 m，巖溶水水位呈逐漸下降趨勢。同時工業廢水、農藥化肥等污染源不斷增多，滲入地下導致地下水污染，地下水水質急劇惡化。根據1996年水質監測資料顯示，東武、樓德水源地水化學類型均以HCO3—Ca型水為主，經2018年調查，東武水源地水化學類型以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-濃度明顯升高，水質狀況以Ⅲ類較差水和Ⅳ類極差水為主，地下水環境保護刻不容緩。

2 地下水水源地保護區劃分方法

2.1 劃分方法與思路

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

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

2.2 數值模擬法—保護區初步劃分

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

2.3 地下水防污性能評價

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

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

(1)確定：

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

(1)

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

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

2.4 地下水污染風險評價

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

P′=KQLD

(2)

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

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

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

污染危害性(V)一般根據研究區內地下水的用途來進行分級，因此本文據此將污染危害性分為3級(見表5)。最終按照表6的評價標準得到污染風險評價結果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 空間疊加分析

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

3 大汶口盆地地下水水源地保護區劃分

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

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

3.1 大汶口盆地數值模擬結果及保護區初步劃分

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

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

(3)

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

模擬區地下水溶質運移用如下數學模型表示：

(4)

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

圖2 溶質運移模擬Fig.2 Solute transport simulation diagram

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

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

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

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

(5)

根據盆地內各巖溶水點的調查數據，結合表1的評分體系和表2的分級標準，將大汶口盆地地下水防污性能分為四級(見表8)。防污性能高和非常高的地區呈帶狀分布在盆地北部、東北部和南部山區，富水性小，第三系蓋層較厚，地下水防污能力強；防污性能低區主要分布在大汶口鎮、樓德鎮周邊地區以及大汶河富水性大的地段，河流對地下水的補給量大，局部地區第三系蓋層缺失，地下水易受地表水的影響；防污性能中等的地區主要分布在防污性能低區的外圍。大汶口盆地地下水防污性能評價空間分布如圖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 大汶口盆地地下水污染風險評價

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

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

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

3.4 空間疊加及保護區范圍確定

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

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

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

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

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

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

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

4 結 論

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

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