基于數(shù)值模型的肥城市地下水源地保護(hù)區(qū)劃分

王曉瑋，王延嶺，趙志偉，陳偉清

(1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083；2.山東省第五地質(zhì)礦產(chǎn)勘查院，山東 泰安 271000)

?

基于數(shù)值模型的肥城市地下水源地保護(hù)區(qū)劃分

王曉瑋1，王延嶺2，趙志偉2，陳偉清2

(1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083；2.山東省第五地質(zhì)礦產(chǎn)勘查院，山東 泰安 271000)

在北方隱伏巖溶分布地區(qū)，許多城市開發(fā)利用巖溶地下水資源作為飲用水水源，但利用數(shù)值模型方法進(jìn)行巖溶水飲用水源保護(hù)區(qū)劃分的工作較少。以山東省肥城市城區(qū)水源地為例，采用數(shù)值模型方法，構(gòu)建地下水流數(shù)值模型，并采用質(zhì)點(diǎn)追蹤技術(shù)對巖溶水飲用水源保護(hù)區(qū)進(jìn)行了劃分。一級保護(hù)區(qū)面積1.53km2，二級保護(hù)區(qū)面積28.4km2，準(zhǔn)保護(hù)區(qū)面積約272.9km2，全部位于肥城市轄區(qū)內(nèi)。此次劃定的保護(hù)區(qū)范圍比之前劃定的保護(hù)區(qū)范圍更合理、準(zhǔn)確，有利于加強(qiáng)地下水資源的管理和保護(hù)。

巖溶水；飲用水源保護(hù)區(qū)；數(shù)值模型；質(zhì)點(diǎn)追蹤；肥城市

0 引言

北方地區(qū)廣泛分布的隱伏巖溶地下水水資源是很多城市重要甚至唯一的飲用水水源。由于巖溶含水介質(zhì)的隱蔽性和分布的復(fù)雜性，污染的途徑和運(yùn)移時(shí)間長，一旦污染不易被及時(shí)發(fā)現(xiàn)[1]。因此該類水源地地下水資源的管理應(yīng)以預(yù)防為主，劃定合理的保護(hù)區(qū)范圍，建立完善的監(jiān)督和管理措施。根據(jù)國家環(huán)境保護(hù)總局發(fā)布的《飲用水水源保護(hù)區(qū)劃分技術(shù)規(guī)范》(HJ/T338-2007)[2]，巖溶水飲用水水源保護(hù)區(qū)的劃分辦法包括經(jīng)驗(yàn)公式法和數(shù)值模型法。與經(jīng)驗(yàn)公式法相比，數(shù)值模型方法可以刻畫復(fù)雜的水文地質(zhì)條件，且能夠較為真實(shí)的還原地下水系統(tǒng)的特征和運(yùn)動(dòng)狀態(tài)[3]，是較為科學(xué)的水源地保護(hù)區(qū)劃分方法[4-6]，通過數(shù)值模擬刻畫巖溶地下水流是可行的[7-9]，但在對北方巖溶地下水源地保護(hù)區(qū)劃分的研究中僅有少數(shù)采用了數(shù)值模型方法[10-12]。

肥城盆地是北方巖溶地區(qū)典型的向斜-盆地型巖溶水系統(tǒng)[13]，盆地內(nèi)肥城市城區(qū)以地下水作為唯一的供水水源地，之前未針對飲用地下水水源地開展過數(shù)值模擬研究[14-15]，只有有關(guān)部門曾采用經(jīng)驗(yàn)法對水源地保護(hù)區(qū)范圍進(jìn)行了劃分[16]。采用數(shù)值模型方法，嘗試建立肥城市城區(qū)水源地三維非穩(wěn)定流地下水?dāng)?shù)值模型，通過質(zhì)點(diǎn)追蹤技術(shù)劃定水源地開采井保護(hù)區(qū)范圍，對有效防止巖溶水源污染，保障城市供水安全具有借鑒意義，可以促進(jìn)水源地保護(hù)和城市建設(shè)共同發(fā)展。

1 水源地概況

1.1 地質(zhì)及水文地質(zhì)條件

肥城盆地位于山東省中部偏西，泰山西麓，肥城市北部，是以肥城市北部平原及其周邊山體為主體形成的獨(dú)立地質(zhì)構(gòu)造單元，在區(qū)域構(gòu)造上隸屬華北型沉降。區(qū)域內(nèi)在原有的變質(zhì)巖基底上沉積了淺海相、陸相和海陸交互相的地層。肥城盆地內(nèi)斷裂構(gòu)造發(fā)育，肥城斷裂、馬山斷裂、桃園斷裂及其伴生的網(wǎng)狀斷層，使得大地構(gòu)造更加復(fù)雜。

肥城市城區(qū)水源地位于肥城斷陷盆地東部(圖1)。屬暖溫帶大陸性半濕潤季風(fēng)氣候區(qū)，四季分明，多年平均降水量633mm，汛期降水量占全年降水量的70%以上，多年平均氣溫為12.9℃。

1—行政界；2—肥城盆地邊界；3—斷層；4—模擬第一含水層范圍：5—模擬第二含水層范圍；6—城區(qū)位置；7—水源地開采井圖1 模擬區(qū)位置圖

肥城市城區(qū)水源地主要取水層為碳酸鹽巖類裂隙巖溶水含水巖組，包括碳酸鹽巖裂隙巖溶水含水層和碳酸鹽巖夾碎屑巖巖溶裂隙水含水層。碳酸鹽巖類裂隙巖溶水分布在水源地大部分地區(qū)，在城區(qū)一帶隱伏于第四系以下，在南部地區(qū)則出露于地表，巖性為奧陶系馬家溝群、寒武系三山子組灰?guī)r、白云質(zhì)灰?guī)r，單井涌水量一般1000～3000m3/d；碳酸鹽巖類夾碎屑巖巖溶裂隙水分布在水源地南部丘陵區(qū)以南，巖性為寒武系炒米店組、崮山組頁巖夾薄板狀灰?guī)r、竹葉狀灰?guī)r，單井涌水量一般500～1000m3/d。含水層主要接受大氣降水入滲補(bǔ)給和側(cè)向滲流補(bǔ)給；水源地地下水由南向北徑流至北部受石炭二疊系弱透水層的阻擋后折向西南；人工開采已成為水源地范圍內(nèi)主要的巖溶水排泄方式。水源地范圍內(nèi)發(fā)育有馬山斷裂、桃園斷裂等傾向E—SE的高角度正斷層，多具導(dǎo)水性，對地下水含水系統(tǒng)起一定的控制作用。

1.2 地下水開發(fā)利用情況

肥城市城區(qū)內(nèi)居民密集，工業(yè)發(fā)達(dá)。肥城市自來水公司在城區(qū)水源地范圍內(nèi)現(xiàn)有生產(chǎn)井和備用井32眼，企業(yè)自備井?dāng)?shù)以百計(jì)。根據(jù)肥城市環(huán)境保護(hù)局統(tǒng)計(jì)結(jié)果，目前肥城市城區(qū)水源地工業(yè)和生活用自來水集中供水開采量約為3.0萬m3/d，而2012年工業(yè)自備井取水量則達(dá)到了1.67萬m3/d。區(qū)內(nèi)地下水動(dòng)態(tài)類型為降水入滲-開采型，近年來年際水位變化總體呈下降趨勢，在城區(qū)范圍內(nèi)產(chǎn)生了大范圍的降落漏斗。2015年肥城市城區(qū)水源地地下水總體質(zhì)量均可達(dá)到《地下水質(zhì)量標(biāo)準(zhǔn)》(GB/T14848-1993)的Ⅲ類水以上標(biāo)準(zhǔn)，但與往年相比總體水質(zhì)呈下降趨勢，特別是硝酸鹽、硫酸鹽等離子含量一直呈上升趨勢。水源地主要污染源包括工業(yè)、農(nóng)業(yè)、生活、固體廢棄物、加油站及油庫等，漏斗區(qū)的形成加重了地下水污染[17]。肥城市曾于2000年開展了城區(qū)水源地保護(hù)區(qū)劃，并于2001年發(fā)布實(shí)施了《肥城市城區(qū)水源地保護(hù)區(qū)管理辦法》，但城市化進(jìn)程的加快，經(jīng)濟(jì)社會(huì)的持續(xù)發(fā)展，必然導(dǎo)致水資源需求量的增長，進(jìn)一步加劇了水資源供需矛盾。

2 數(shù)值模型方法的應(yīng)用

2.1 水文地質(zhì)概念模型與數(shù)學(xué)模型

該次研究劃定的模擬區(qū)范圍包括肥城市城區(qū)水源地補(bǔ)給和徑流區(qū)，東北和東部至肥城-王莊弧形斷裂，南至地表分水嶺以南寒武系灰?guī)r與泰山巖群地層界線，北至肥城煤田邊界，西南至孫牛公路，西部以行政區(qū)劃為界。

垂向上將模擬區(qū)內(nèi)含水層概化為兩個(gè)含水層，第一含水層為灰?guī)r裸露區(qū)及第四系，概化為潛水，第二含水層為隱伏灰?guī)r，概化為承壓水。一層北部山前沖洪積平原概化為定流量流入邊界；西部沿河道方向?yàn)槎髁苛鞒鲞吔?；東部和南部為天然隔水邊界；西南部根據(jù)流場，判斷為隔水邊界。二層?xùn)|部邊界定為肥城盆地的東部零流量邊界；南部邊界處主要為寒武系灰?guī)r，賦水性較差，巖溶連通性一般，概化為零流量邊界；中北部邊界為老城鎮(zhèn)的煤田南部邊界，概化為定流量流出邊界，模型中分為3段；西部邊界本身為行政區(qū)邊界及其延伸，且為徑流帶的流出邊界，概化為定流量邊界。含水層之間存在一定的水量交換。

整個(gè)地下水系統(tǒng)接受大氣降水入滲、農(nóng)業(yè)回灌、側(cè)向流入的補(bǔ)給，主要排泄方式是人工開采，以及少量的側(cè)向流出和排泄到地表河流。將模擬區(qū)含水系統(tǒng)概化為非均質(zhì)、垂向各項(xiàng)異性、三維非穩(wěn)定地下水流系統(tǒng)。其控制方程為：

潛水：

式中：Kx，Ky，Kz分別為潛水含水層水平x，y方向和垂向滲透系數(shù)(m/d)，μ為給水度，h為含水層的水位標(biāo)高(m)，ε為含水層的源匯項(xiàng)(L/d)。

承壓水：

式中：Kx，Ky，Kz—分別為承壓含水層水平x，y方向和垂向滲透系數(shù)(m/d)，S為貯水系數(shù)(L/m)，H為含水層的水位標(biāo)高(m)，其他參數(shù)意義同前。

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

模擬區(qū)面積357.8km2，剖分為100m×100m的單元格，選取2011年6月1日至2015年3月31日為該次模型研究的識(shí)別階段，以2015年4月1日至5月31日為驗(yàn)證階段。以模擬區(qū)2011年6月的地下水位實(shí)測值為二層承壓水模擬水位的初始數(shù)據(jù)，以2012年6月的地下水位實(shí)測值近似替代一層潛水模擬水位的初始數(shù)據(jù)。

采用2014年4月19日的Landsat8遙感數(shù)據(jù)，使用ENVI5.1軟件進(jìn)行解譯，采用基于CART的自動(dòng)決策樹分類方法，確定了工作區(qū)范圍內(nèi)的地物類型，結(jié)合巖性特征，綜合確定了一層的參數(shù)分區(qū)，分為12個(gè)區(qū)。二層參數(shù)分區(qū)主要考慮不同巖性區(qū)域的賦水性差異、區(qū)內(nèi)馬山斷裂構(gòu)造的影響等，分為5個(gè)區(qū)。

該次模擬的源匯項(xiàng)主要包括降水入滲、農(nóng)業(yè)灌溉開采、企業(yè)自備井開采和水源地開采。將月降雨量轉(zhuǎn)化為平均降雨入滲強(qiáng)度，根據(jù)前述遙感解譯結(jié)果和地形資料，將降水入滲補(bǔ)給的范圍，根據(jù)補(bǔ)給強(qiáng)度的差異分為4個(gè)區(qū)。區(qū)內(nèi)農(nóng)灌井的數(shù)量眾多，無詳細(xì)井位和開采量數(shù)據(jù)，因此確定4個(gè)主要的井灌區(qū)，按照200m×200m間距布置開采井，根據(jù)《山東省主要農(nóng)作物灌溉定額》(DB37/T1640-2010)中的規(guī)定，核算不同降水保證率下畝均灌溉定額，根據(jù)灌溉制度分解為每月的地下水開采量，核算出模擬期內(nèi)每個(gè)月的農(nóng)灌開采量用于計(jì)算。將收集到的工業(yè)和水源地開采量進(jìn)行處理后，分配到逐月用于模型計(jì)算。

根據(jù)模型最終調(diào)試結(jié)果分析，識(shí)別后的模型概化和參數(shù)符合實(shí)際的水文地質(zhì)條件。模擬的地下水位變化過程與趨勢同實(shí)際的地下水動(dòng)態(tài)過程基本一致，模擬的地下水流場基本可以反映區(qū)域內(nèi)地下水開采形成的降落漏斗(圖2、圖3)。

圖2 觀測孔BH25水位識(shí)別期過程線擬合情況

圖3 模型模擬的識(shí)別期末二層地下水流場擬合情況

2.3 模型的預(yù)測

以數(shù)值模型識(shí)別階段取得的末流場為初始流場，進(jìn)行一個(gè)水文年的地下水流場的預(yù)測。預(yù)測過程中，邊界條件不發(fā)生變化；在源匯項(xiàng)中，選擇保證率為50%的年份的降雨量作為降雨入滲和農(nóng)業(yè)開采的計(jì)算依據(jù)，工業(yè)開采量保持不變，水源地開采量除備用井外，選擇2014年的開采量數(shù)據(jù)作為水井開采量，目前尚未完全使用的華美熱電1-3號井、西水廠1號井、象山井、桃園街小學(xué)井等，類比附近開采井的開采量設(shè)定。取得的模擬結(jié)果作為模擬水源井跡線的依據(jù)。

3 水源地保護(hù)區(qū)劃分

3.1 劃分技術(shù)方法

對三維穩(wěn)定流，modpath的質(zhì)量平衡方程可表示為：

(3)

式中：Vx，Vy，Vz—位線形流動(dòng)流速矢量在各坐標(biāo)軸方向的分量，N為含水層有效孔隙率，W—由含水層單位體積源、匯產(chǎn)生的水量，由GMS依據(jù)公式(1)(2)計(jì)算得出。對于三維非穩(wěn)定流,可以將其視為由一系列的穩(wěn)定流時(shí)間步長組成，在每一個(gè)時(shí)間步長內(nèi)按穩(wěn)定流方法計(jì)算質(zhì)點(diǎn)運(yùn)移矢量，通過矢量累加得到非穩(wěn)定流一定時(shí)間段內(nèi)的質(zhì)點(diǎn)運(yùn)移軌跡[12]。

以modflow中各個(gè)水井所在網(wǎng)格的中心點(diǎn)作為水質(zhì)點(diǎn)的起始位置。利用modpath的反向示蹤功能按100d，1000d確定水質(zhì)點(diǎn)位置。連接質(zhì)點(diǎn)位置構(gòu)成地下水流截獲區(qū)。截獲區(qū)主要沿地下水流向展布，質(zhì)點(diǎn)在水源地上游運(yùn)移跡線較長，呈條帶狀展布；在水源地中游運(yùn)移跡線較短，呈扇形展布。

由于肥城城區(qū)水源地屬于群井水源地，根據(jù)技術(shù)規(guī)范對集中式供水水源地保護(hù)區(qū)范圍的規(guī)定，井群內(nèi)井間距小于或等于保護(hù)區(qū)半徑的2倍時(shí)，以外圍井的外界多邊形為邊界，向外徑向距離為保護(hù)區(qū)半徑的多邊形區(qū)域。對modpath確定的地下水流截獲區(qū)進(jìn)行修正。并結(jié)合水源地保護(hù)區(qū)附近的地標(biāo)、地物特點(diǎn)，充分利用具有永久性的明顯標(biāo)志，最終確定各級保護(hù)區(qū)的界線。

3.2 保護(hù)區(qū)劃分結(jié)果

保護(hù)區(qū)劃定結(jié)果為：一級保護(hù)區(qū)面積合計(jì)1.53km2，二級保護(hù)區(qū)面積合計(jì)28.4km2，準(zhǔn)保護(hù)區(qū)面積約272.9km2。以二級保護(hù)區(qū)為例，計(jì)算得到的1000d地下水流截獲區(qū)面積為3.36km2?？紤]到井群間距，向外徑向出保護(hù)區(qū)半徑，并結(jié)合附近地標(biāo)地物，最終確定以北部康匯河河道和南部的道路為保護(hù)區(qū)邊界，圈定保護(hù)區(qū)面積28.4km2。

圖4 1000d截獲區(qū)范圍及劃定的二級保護(hù)區(qū)示意圖

3.3 分區(qū)可靠性評價(jià)

與前人成果相比[16]，該次研究得到的保護(hù)區(qū)類型從6個(gè)減少為3個(gè)，面積變化明顯，一級保護(hù)區(qū)面積減少為1.53km2，二級保護(hù)區(qū)面積從99.2km2減少為28.43km2，準(zhǔn)保護(hù)區(qū)面積從44.6km2擴(kuò)大為272.9km2。一、二級保護(hù)區(qū)的范圍主要位于建城區(qū)，通過合理計(jì)算得到較小的保護(hù)區(qū)面積，有利于地方政府和管理部門安排開展有關(guān)建設(shè)和保護(hù)工作，減少城市發(fā)展和水源地保護(hù)的沖突。

4 結(jié)論

(1)該次研究重點(diǎn)針對肥城市城區(qū)水源地開展數(shù)值模擬工作，為了彌補(bǔ)資料欠缺，參考了遙感影像、灌溉定額等數(shù)據(jù)，對緣匯項(xiàng)及水文地質(zhì)參數(shù)進(jìn)行厘定，對水文地質(zhì)概念模型進(jìn)行了較為詳細(xì)準(zhǔn)確的刻畫，具有一定的先進(jìn)性。

(2)數(shù)值模型建立的關(guān)鍵在于巖溶含水介質(zhì)的合理概化以及資料的詳細(xì)程度。該次研究較為準(zhǔn)確的概化了城區(qū)水源地的水文地質(zhì)情況，但在水文地質(zhì)參數(shù)、開采量、水位觀測數(shù)據(jù)等資料獲取方面存在欠缺，影響了模型評價(jià)的精度。另外現(xiàn)有的等效多孔介質(zhì)模型，在巖溶地區(qū)使用時(shí)忽略了巖體的非均質(zhì)性，本身極易造成較大的計(jì)算誤差。以上因素降低了分區(qū)的可靠性。未來巖溶地區(qū)數(shù)值模型應(yīng)用中引入不確定性分析，可在一定程度上提高評價(jià)結(jié)果的可靠性。

(3)建議在保護(hù)區(qū)范圍內(nèi)建立完善的地下水監(jiān)測網(wǎng)絡(luò)，結(jié)合當(dāng)前開展的國家地下水監(jiān)測工程和省級地下水監(jiān)測工程等項(xiàng)目，對水源地水位及水質(zhì)實(shí)現(xiàn)動(dòng)態(tài)觀測。按照《飲用水水源保護(hù)區(qū)污染防治管理規(guī)定》相關(guān)要求，結(jié)合肥城市城區(qū)水源地環(huán)境現(xiàn)狀，針對不同級別的保護(hù)區(qū)，采取不同的防范措施和管理措施。大力控制地下水開采量，多方尋找替代水源，提高工業(yè)的單位用水效率，提倡居民生活節(jié)約用水，避免因過度開采引起地下水動(dòng)力條件發(fā)生進(jìn)一步的變化，對水源地水質(zhì)造成污染。

[1] 郝永艷.三給地壘巖溶地下水系統(tǒng)及水源地保護(hù)區(qū)劃分研究[D].太原：太原理工大學(xué),2011.

[2] HJ/T338-2007.飲用水水源地保護(hù)區(qū)劃分技術(shù)規(guī)范[S].

[3] 劉志濤,周群道,楊建華.地下水溶質(zhì)運(yùn)移數(shù)值法和解析法預(yù)測結(jié)果對比分析——以沾化電廠為例[J].山東國土資源,2016,32(7):78-82.

[4] 趙紅梅,肖杰.公式法與數(shù)值模擬法在地下水飲用水源保護(hù)區(qū)劃分中的應(yīng)用——成都平原某水源地為例[J].四川環(huán)境,2013(S1):60-64.

[5] 鄧媛媛,胡立堂,高童,等.吳忠市金積地下水飲用水源地保護(hù)區(qū)劃分[J].南水北調(diào)與水利科技,2013(01):127-131.

[6] 湯新梅.北京市水源八廠水源地保護(hù)區(qū)劃分研究[D].長春：吉林大學(xué),2012.

[7] 付延玲.邯邢水文地質(zhì)南單元巖溶地下水系統(tǒng)數(shù)值預(yù)報(bào)[J].中國巖溶,2002,21(4):269-275.

[8] 錢家忠,吳劍鋒,董洪信.徐州市張集水源地裂隙巖溶水三維等參有限元數(shù)值模擬[J].水力學(xué)報(bào),2003(3):37-41.

[9] 劉猛,束龍倉,劉波.地下水?dāng)?shù)值模擬中的參數(shù)隨機(jī)模擬[J].水利水電科技進(jìn)展,2005,25(6):25-27.

[10] 卜華,陳占成,張良鵬.山東羊莊巖溶水系統(tǒng)飲用水水源地保護(hù)區(qū)劃分探討[J].科技創(chuàng)新導(dǎo)報(bào),2008(20):110-112.

[11] 翟立娟.巖溶水飲用水水源保護(hù)區(qū)劃分技術(shù)方法——以邯鄲市羊角鋪水源地為例[J].中國巖溶,2011(01):47-52.

[12] 李星宇,南天,王新娟,等.數(shù)值模擬方法在隱伏巖溶水源地保護(hù)區(qū)劃分及污染治理中的應(yīng)用[J].中國巖溶,2014(03):280-287.

[13] 梁永平,王維泰.中國北方巖溶水系統(tǒng)劃分與系統(tǒng)特征[J].地球?qū)W報(bào),2010,31(06):860-868.

[14] 魏曉燕.肥城盆地地下水動(dòng)態(tài)數(shù)值模擬研究[D].濟(jì)南：濟(jì)南大學(xué),2014.

[15] 王瑋,李云峰,侯東輝,等.肥城盆地石橫電廠水源地巖溶水系統(tǒng)地下水管理模型[J].地球科學(xué)與環(huán)境學(xué)報(bào),2004,26(3):32-39.

[16] 王萬喜,王啟田,劉福臣,等.水源地保護(hù)區(qū)劃與防護(hù)措施——以肥城盆地為例[J].地下水,2008,30(04):42-44.

[17] 崔素芳,張保祥,范明元,等.肥城盆地地下水水化學(xué)演變規(guī)律研究[J].人民黃河,2015,37(3):75-79.

Partition of Groundwater Source Protection Zones of Feicheng City Based on Numerical Model

WANG Xiaowei1，WANG Yanling2，ZHAO Zhiwei2，CHEN Weiqing2

(1. School of Water Resource and Environment, China University of Geosciences (Beijing), Beijing 100083，China；2. No.5 Exploration Institute of Geology and Mineral Resources, Shandong Tai'an 271000, China)

Karst groundwater had been developed as the drinking water source in many cities in the areas of covered karst areasof northern China. There remains lots of work to do in demarcating karst water drinking source protection zones with the technique of numerical model. This research conducted the partition of karst groundwater drinking source protection zones of Feicheng City in Shandong Province by building a groundwater numerical model and executing particle tracing. The square of the first level protection zone for the drinking water source is 1.53km2, the second level is 28.4km2, and the square beyond protection zone is about 272.9km2. It is showed that protection zones divided by numerical model are reliable and rational compared with the former division. It is suitable for management and protection of groundwater resources by local government.

Karst water; protection of drinking water source; numerical model; particle tracing; Feicheng city

2016-11-21；

2017-01-20；編輯：曹麗麗 作者簡介：王曉瑋(1982—)，男，山東泰安人，博士研究生，主要研究方向?yàn)榈叵滤Y源評價(jià)和管理；E-mail：wangxw@cugb.edu.cn

P641.8

B

王曉瑋，王延嶺，趙志偉，等.基于數(shù)值模型的肥城市地下水源地保護(hù)區(qū)劃分[J].山東國土資源，2017,33(4)：29-33.WANG Xiaowei，WANG Yanling，ZHAO Zhiwei，etc. Partition of Groundwater Source Protection Zones of Feicheng City Based on Numerical Model[J].Shandong Land and Resources, 2017,33(4)：29-33.