GMS在莒縣城區(qū)地下水水源地保護(hù)區(qū)劃分中的應(yīng)用

(吉林大學(xué) 新能源與環(huán)境學(xué)院，吉林 長(zhǎng)春 130021)

飲用水水質(zhì)直接關(guān)系到人們的生命健康，習(xí)近平總書(shū)記“4·13”重要講話中提到：全面開(kāi)展集中式飲用水源地專項(xiàng)整治，加強(qiáng)近岸海域污染治理，盡快實(shí)現(xiàn)城鎮(zhèn)污水管網(wǎng)全覆蓋、全收集、全處理[1]。地下水水源地保護(hù)區(qū)劃分的主要目的是為了避免水源地附近的人類活動(dòng)及其所排放污染物對(duì)水源造成污染而劃定的保護(hù)區(qū)域，因此在理想條件下，應(yīng)在人類活動(dòng)區(qū)域與地下水水源地之間預(yù)留足夠的空間，在人為和自然條件下使地下水得到適當(dāng)處理，從而達(dá)到人們生活和飲用的標(biāo)準(zhǔn)。此外，在水源地出現(xiàn)意外污染時(shí)，也可作為處理和修復(fù)的區(qū)域。建立水源保護(hù)區(qū)，是為了更好地保護(hù)水源地水質(zhì)安全和應(yīng)對(duì)突發(fā)情況做準(zhǔn)備[2]。

水源地保護(hù)區(qū)劃分的常用方法有公式法、簡(jiǎn)化圖形法、分析法、解析模型法和數(shù)值法等[3-6]，其中數(shù)值法通過(guò)構(gòu)建地下水模型，并計(jì)入溶質(zhì)運(yùn)移模型，以達(dá)到劃分保護(hù)區(qū)的目的。已有的地下水保護(hù)區(qū)劃分?jǐn)?shù)值法方面的研究[7-10]大多只建立了水流模型，然后利用粒子示蹤來(lái)劃分保護(hù)區(qū)范圍，比如王涵等[10]利用MODFLOW建立地下水滲流場(chǎng)，利用MODPATH進(jìn)行粒子反向示蹤模擬，對(duì)呼和浩特市城市的水源地進(jìn)行了保護(hù)區(qū)劃分；韓京龍等[11]使用同樣的方法對(duì)磐石市地下水水源地進(jìn)行了保護(hù)區(qū)劃分?？紤]到污染物在隨地下水進(jìn)行遷移的過(guò)程中，由于自身分子擴(kuò)散的作用，污染物的運(yùn)移速度要快于地下水流的運(yùn)移速度，如果僅進(jìn)行地下水流模擬，并不能真正模擬出污染物的運(yùn)移范圍。因此，以山東省莒縣城區(qū)相近的三個(gè)水源地為研究對(duì)象，根據(jù)研究區(qū)的水文地質(zhì)條件建立地下水的水流模型，在水流模型的基礎(chǔ)上進(jìn)行溶質(zhì)運(yùn)移模擬[11-18]，得到三個(gè)水源地的保護(hù)區(qū)范圍，為莒縣城區(qū)的地下水水源地保護(hù)和開(kāi)發(fā)利用提供了可靠的依據(jù)。

1 研究區(qū)概況

研究區(qū)位于山東省東南部日照市莒縣，南北長(zhǎng)16 km，東西寬10 km，面積約160 km2，如圖1所示。地形較為簡(jiǎn)單，東、西部為丘陵，中部為沖積平原區(qū)，地形總趨勢(shì)為北部高、南部低。平原地形標(biāo)高95～110 m，地形坡度2～14°。地下水主要以松散巖類孔隙裂隙水為主，由第四系砂及砂礫石組成，厚度20～30 m，含水量大，單井出水量最低可達(dá)1 000 m3/d以上，淺層地下水的埋深平水期為3～5 m，豐水期更淺，為2～4 m。包氣帶巖性以粉質(zhì)粘土為主，并有粉土、粉砂夾層。地下水主要為降水入滲補(bǔ)給和河流滲漏補(bǔ)給，農(nóng)業(yè)開(kāi)采和工業(yè)開(kāi)采為主要排泄項(xiàng)，以地下水作為飲用水源，地下水水源地主要有3個(gè)，由北向南分別為凈水廠、第一水廠和文心水廠。過(guò)凈水廠由南向北的剖面圖如圖2所示。三個(gè)水廠水井井徑約8 m，井深60 m，動(dòng)態(tài)水位范圍為80～100 m。

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

根據(jù)研究區(qū)基礎(chǔ)水文地質(zhì)條件、水動(dòng)力場(chǎng)、水化學(xué)場(chǎng)分析確定概念模型。整個(gè)莒縣沭河盆地的上層除河谷地區(qū)外均為弱透水性的黏土層，中間層為第四系孔隙潛水砂層，底部為隔水的厚層砂巖。

2.1 確定模型模擬區(qū)域

確定模擬區(qū)域要把可能影響到的區(qū)域盡量包括進(jìn)來(lái)，故模擬區(qū)東到第四系與基巖的分界處，西到西部沭河和柳清河之間的分水嶺處，北到具有相對(duì)穩(wěn)定地下水位的定水頭邊界處，南到沭河及其一條支流的定水頭邊界處。

2.2 地層概化

根據(jù)研究區(qū)的水文地質(zhì)概況，整個(gè)莒縣沭河盆地的上層除河谷地區(qū)外為弱透水性的黏土層，中間層為第四系孔隙潛水砂層，底部為隔水的厚層砂巖。因此，把地層概化為3層：第一層為透水性弱的第四系孔隙潛水層，第二層為透水性強(qiáng)的第四系孔隙潛水層，第三層為隔水層。

3 地下水水流數(shù)值模擬

根據(jù)莒縣的水文地質(zhì)條件對(duì)研究區(qū)進(jìn)行模型概化，利用GMS軟件中的MODFLOW模塊建立地下水水流的數(shù)值模型。

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

圖2 地層剖面圖Fig. 2 Stratigraphic profile

3.1 地下水水流數(shù)學(xué)模型

根據(jù)地質(zhì)條件和水文地質(zhì)條件分析，莒縣城區(qū)地下水流可視作平面二維流，對(duì)應(yīng)的地下水流數(shù)學(xué)模型為：

(1)

其中：h—地下水水頭，m；Kx、Ky—x、y方向滲透系數(shù)，m/d；h1—含水層第一類邊界上的水頭，m；q—含水層第二類邊界上的流量，m3/d；h0—含水層初始水頭，m；W—源匯項(xiàng)強(qiáng)度(包括開(kāi)采強(qiáng)度等)，m3/d；Σ1—含水層第一類邊界；Σ2—含水層第二類邊界；D—研究區(qū)；μ—給水系數(shù)。

3.2 邊界條件處理

研究區(qū)存在一類水頭邊界和二類流量邊界。南部為一條河流，常年有水，與地下水存在水力聯(lián)系，故處理為水頭邊界。北部有多個(gè)地下水位觀測(cè)孔，根據(jù)長(zhǎng)期水位觀測(cè)資料顯示，在沭河流域進(jìn)行地下水開(kāi)采對(duì)該處地下水位無(wú)明顯影響，因此可處理為水頭邊界。在研究區(qū)的西部有一地表分水嶺，位于沭河和柳清河之間，兩側(cè)的地下水水位高于河流兩側(cè)水位，所以該處可作為研究區(qū)的地下隔水邊界。研究區(qū)東部的定流量邊界是第四系和基巖的交界線，由于基巖區(qū)地勢(shì)高且透水性較差，大氣降水沿基巖表面經(jīng)地表徑流匯入研究區(qū)內(nèi)的含水層內(nèi)，故定義為給定流量邊界。

3.3 模型區(qū)域剖分

根據(jù)研究區(qū)的水文地質(zhì)條件，研究區(qū)在垂向上分為3層，水平方向分為128×187的網(wǎng)格，每個(gè)網(wǎng)格代表100×100的區(qū)域。

3.4 源匯項(xiàng)的處理

研究區(qū)內(nèi)潛水含水層的補(bǔ)給方式主要包括降水入滲補(bǔ)給、農(nóng)業(yè)水灌溉補(bǔ)給和河流側(cè)向徑流補(bǔ)給，排泄方式主要為蒸發(fā)排泄和人工開(kāi)采。人工開(kāi)采按不同開(kāi)采分區(qū)的開(kāi)采總量平均分布到每個(gè)開(kāi)采分區(qū)上，開(kāi)采分區(qū)圖如圖3所示，每個(gè)分區(qū)的開(kāi)采量如表1所示；根據(jù)研究區(qū)降雨情況，降水入滲系數(shù)取0.125～0.3，得到不同時(shí)間的降雨補(bǔ)給強(qiáng)度，見(jiàn)圖4；灌溉回滲強(qiáng)度為0.000 022 m/d；多年平均蒸發(fā)為1 410 mm，蒸發(fā)系數(shù)為0.04，蒸發(fā)強(qiáng)度為0.000 016 m/d，蒸發(fā)極限埋深為2.5 m。

圖3 研究區(qū)開(kāi)采分區(qū)圖Fig. 3 Zoning map of mining in the study area

圖4 不同時(shí)間段下的降雨補(bǔ)給強(qiáng)度Fig. 4 Rainfall recharge intensity at different time periods

表1 不同開(kāi)采分區(qū)的開(kāi)采強(qiáng)度

Tab. 1 Mining intensity of different mining zones m/d

分區(qū)一區(qū)二區(qū)三區(qū)四區(qū)開(kāi)采量0.000 024 70.000 140 00.000 411 00.000 685 0

3.5 水文地質(zhì)參數(shù)

主要為滲透系數(shù)K、潛水的給水度、承壓水的貯水系數(shù)和降雨入滲系數(shù)α等，各水文地質(zhì)參數(shù)按相應(yīng)試驗(yàn)、地層巖性和參考經(jīng)驗(yàn)值等途徑綜合給定初始值，并通過(guò)模型調(diào)參，進(jìn)而得到研究區(qū)模型的水文地質(zhì)參數(shù)。

3.6 模型的識(shí)別驗(yàn)證

在上述條件下運(yùn)行模型，可得到研究區(qū)的模擬地下水流場(chǎng)圖。通過(guò)與相同時(shí)期的真實(shí)地下水流場(chǎng)圖進(jìn)行對(duì)比擬合，對(duì)相關(guān)水文地質(zhì)參數(shù)進(jìn)行合理化調(diào)整，使模型與實(shí)際情況更加契合。該模型進(jìn)行擬合的主要原則為：在水文地質(zhì)參數(shù)滿足研究區(qū)水文地質(zhì)條件的基礎(chǔ)上使模型的等水位線和水均衡中各源匯項(xiàng)與實(shí)際相符。取兩個(gè)長(zhǎng)期觀測(cè)點(diǎn)對(duì)2016年的觀測(cè)水位進(jìn)行模型識(shí)別，取2017年豐水期的地下水等水位線進(jìn)行驗(yàn)證。經(jīng)過(guò)識(shí)別和驗(yàn)證，得到研究區(qū)水文地質(zhì)參數(shù)值如表2所示。其中第一層的滲透系數(shù)(巖性)分區(qū)如圖5所示，第二層和第三層滲透系數(shù)未分區(qū)，第二層、第三層均為砂巖，每個(gè)分區(qū)的滲透系數(shù)見(jiàn)表2。

取2017年豐水期的觀測(cè)水位進(jìn)行模型驗(yàn)證。圖6為模型計(jì)算的2017年豐水期流場(chǎng)與實(shí)測(cè)流場(chǎng)的對(duì)比，模型模擬的水位和模擬的流場(chǎng)與實(shí)際水位和實(shí)際流場(chǎng)也比較接近。經(jīng)過(guò)模型驗(yàn)證，表明本次模擬具有較高的精度和可信性。

表2 研究區(qū)地層水文地質(zhì)參數(shù)

單位：k—m/d；μ—m3/d

圖5 第一層含水層滲透系數(shù)分區(qū)圖Fig. 5 Zoning map of permeability coefficient of the first aquifer

圖6 模型驗(yàn)證期實(shí)測(cè)水位與 計(jì)算流場(chǎng)的對(duì)比圖Fig. 6 Comparison of measured water level and calculated flow field during model validation period

圖7 保護(hù)區(qū)范圍圖Fig. 7 Scope map of the protected areas

4 溶質(zhì)運(yùn)移模擬

在上述地下水水流模型的基礎(chǔ)上，應(yīng)用GMS軟件中的MT3DMS模塊進(jìn)行地下水的溶質(zhì)運(yùn)移模擬，用以計(jì)算污染物的捕獲區(qū)范圍為保護(hù)區(qū)范圍。由于NaCl性質(zhì)穩(wěn)定，一般不參與物理化學(xué)反應(yīng)，本身對(duì)環(huán)境影響小，在保護(hù)區(qū)劃分過(guò)程中更常用。依據(jù)《飲用水水源保護(hù)區(qū)劃分技術(shù)規(guī)范》(HJ/T338—2018)[4],模型以地下水開(kāi)采井為中心，NaCl做為溶質(zhì)，以NaCl遷移100 d的范圍做為水源地一級(jí)保護(hù)區(qū)范圍；在一級(jí)保護(hù)區(qū)外，以NaCl遷移1 000 d的范圍做為二級(jí)保護(hù)區(qū)范圍。

4.1 地下水溶質(zhì)運(yùn)移數(shù)學(xué)模型

地下水溶質(zhì)運(yùn)移的對(duì)流-彌散方程：

(2)

其中：Rd—阻滯因子；c—地下水因子濃度，mg/L；t—時(shí)間，d；xi—坐標(biāo)軸方向上的距離，m；Dij—地下水水動(dòng)力彌散系數(shù)；vi—地下水滲流速度，m/d；qs—源和匯的單位流量，m3/d；cs—源和匯的濃度，mg/L；θ—含水層孔隙率，%。

4.2 模擬過(guò)程

保護(hù)區(qū)范圍應(yīng)以開(kāi)采井為中心，地下水在開(kāi)采過(guò)程中，周?chē)廴疚锵蜷_(kāi)采井中心移動(dòng)，以100 d和1 000 d污染物能運(yùn)移到開(kāi)采井中心的位置點(diǎn)組成的范圍為保護(hù)區(qū)邊界，但在實(shí)際應(yīng)用過(guò)程中不可能對(duì)每個(gè)點(diǎn)都進(jìn)行模擬預(yù)測(cè)來(lái)尋找保護(hù)區(qū)邊界，因此在應(yīng)用過(guò)程中轉(zhuǎn)換思路，將水源地的開(kāi)采井變換為注水井，注水量與開(kāi)采量相同，給開(kāi)采井設(shè)置一個(gè)初始污染物濃度，模擬污染物從開(kāi)采井向四周運(yùn)移的過(guò)程，以污染物從開(kāi)采井向四周運(yùn)移100 d和1 000 d形成的污染暈為保護(hù)區(qū)邊界，本次模擬過(guò)程中注水量設(shè)置為10 000 m3/d，污染源濃度設(shè)置為10 000 mg/L,根據(jù)場(chǎng)地的水文地質(zhì)情況，參照吉林省雙陽(yáng)區(qū)彌散試驗(yàn)場(chǎng)地所得出的結(jié)果，對(duì)比地層巖性、水流速度等參數(shù)，取此研究區(qū)彌散系數(shù)為50 m，模擬總時(shí)長(zhǎng)為1 440 d(4年)。本次研究中的三個(gè)水源地均為孔隙潛水型水源地，主要取水層位均為第二層的砂層，因此主要關(guān)注污染物在模型第二層中的運(yùn)移范圍。最終模擬的溶質(zhì)運(yùn)移100 d和1 000 d的污染暈范圍如圖7和表3所示。

表3 水源地溶質(zhì)運(yùn)移各方向位移表

由模擬結(jié)果可以看出：隨著時(shí)間的增加，污染物的污染暈逐漸擴(kuò)大，而且沿水流方向污染暈的遷移速度更快，說(shuō)明污染物的運(yùn)移擴(kuò)散主要受地下水流動(dòng)的影響，同時(shí)在其他方向污染物也發(fā)生了遷移，說(shuō)明污染物的運(yùn)移擴(kuò)散除了受水流作用的影響外還受污染物自身分子擴(kuò)散作用的影響。

在保護(hù)區(qū)實(shí)際管理操作中，考慮到三個(gè)水源地相距較近，參照計(jì)算結(jié)果，可以將獨(dú)立的三個(gè)水源地二級(jí)保護(hù)區(qū)適當(dāng)外擴(kuò)邊界，整體劃入一個(gè)沿河谷的條帶區(qū)域內(nèi)(可作為準(zhǔn)保護(hù)區(qū))進(jìn)行保護(hù)管理。

5 結(jié)論

以莒縣城區(qū)的三個(gè)水源地為研究對(duì)象，利用GMS中MODFLOW和MT3DMS模塊，對(duì)水源地進(jìn)行保護(hù)區(qū)劃分，得出結(jié)論如下：

1) 首先應(yīng)用GMS軟件的MODFLOW模塊建立了莒縣城區(qū)的地下水水流數(shù)值模型，通過(guò)等水位線的擬合情況驗(yàn)證了模型的合理性，說(shuō)明建立的模型能夠準(zhǔn)確刻畫(huà)保護(hù)區(qū)的范圍。

2) 在水流模型的基礎(chǔ)上再通過(guò)MT3DMS模塊建立了研究區(qū)的溶質(zhì)運(yùn)移模型，得到溶質(zhì)運(yùn)移100 d和1 000 d的污染暈范圍：依據(jù)《飲用水水源保護(hù)區(qū)劃分技術(shù)規(guī)范》(HJ/T338—2018)，將溶質(zhì)運(yùn)移100 d和1 000 d的范圍劃分為一級(jí)保護(hù)區(qū)和二級(jí)保護(hù)區(qū)，凈水廠、第一水廠和文心水廠的一級(jí)保護(hù)區(qū)和二級(jí)保護(hù)區(qū)的污染暈范圍面積分別為34 617和150 800 m2、54 467和203 000 m2、32 200和164 000 m2。

3) 在實(shí)施保護(hù)區(qū)劃分時(shí)，考慮到三個(gè)水源地相距較近，參照計(jì)算結(jié)果，實(shí)際操作時(shí)，可以將三個(gè)獨(dú)立的二級(jí)保護(hù)區(qū)適當(dāng)外擴(kuò)將三個(gè)水源地整體劃入一個(gè)沿河谷的條帶區(qū)域內(nèi)進(jìn)行保護(hù)管理。