付曉剛 唐仲華 劉彬濤 藺林林 卜華 閆佰忠 張永禮

摘要：羊莊盆地地下水資源豐富，近年來盆地內(nèi)地下水硝酸鹽污染呈加重趨勢(shì)?；赩isual Modflow軟件建立了盆地地下水水流與溶質(zhì)運(yùn)移數(shù)學(xué)模型，對(duì)羊莊斷塊地下水硝酸鹽的遷移與擴(kuò)散進(jìn)行了數(shù)值模擬，預(yù)測(cè)了5a和10 a后硝酸鹽在地下水中的分布狀況以及濃度變化趨勢(shì)。結(jié)果表明：在預(yù)測(cè)時(shí)段，地下水中硝酸鹽濃度呈明顯上升趨勢(shì)，模擬區(qū)硝酸鹽濃度較高的地區(qū)主要分布在轆轆井區(qū)域，在地下水動(dòng)力場(chǎng)的影響下，硝酸鹽高濃度區(qū)域主要向西南方向擴(kuò)展。氮肥的過量施用和盆地內(nèi)養(yǎng)殖場(chǎng)糞便的隨意堆放是地下水硝酸鹽濃度持續(xù)增高的主要原因。

關(guān)鍵詞：數(shù)值模擬：Visual Modflow：硝酸鹽污染：羊莊盆地

中圖分類號(hào)：P641.8

文獻(xiàn)標(biāo)志碼：A

doi： 10.3969/j.issn.1000-1379.2018.06.020

羊莊盆地是滕州市乃至棗莊市地下水資源最為豐富的優(yōu)質(zhì)飲用水源區(qū)[1]，是滕州市和棗莊市新城區(qū)的城市供水水源地，地下水資源的開發(fā)為當(dāng)?shù)亟?jīng)濟(jì)發(fā)展和人民生活水平提高做出了重要貢獻(xiàn)，但是受人類活動(dòng)影響，羊莊盆地地下水中硝酸鹽污染呈加劇趨勢(shì)。含硝酸鹽的水源對(duì)人體危害較大，可導(dǎo)致高鐵蛋白癥[2]，因此從供水安全角度出發(fā)，對(duì)羊莊盆地硝酸鹽的動(dòng)態(tài)變化規(guī)律及影響因素進(jìn)行全面分析研究顯得十分迫切。

影響地下水硝酸鹽動(dòng)態(tài)變化的因素復(fù)雜，地下水?dāng)?shù)值模擬可綜合水文地質(zhì)和T農(nóng)業(yè)、生活污染等方面的數(shù)據(jù)，直接建立其相互間的聯(lián)系，使得人們能夠定量分析和預(yù)測(cè)含水層中硝酸鹽濃度的動(dòng)態(tài)變化規(guī)律[3]。因此選用地下水?dāng)?shù)值模擬軟件Visual Modflow建立地下水流和溶質(zhì)運(yùn)移耦合模型[4]，對(duì)羊莊盆地地下水中硝酸鹽的遷移與擴(kuò)散進(jìn)行數(shù)值模擬，預(yù)測(cè)含水層中硝酸鹽濃度的變化，并對(duì)其主要影響因素進(jìn)行系統(tǒng)分析，以期為地下水硝酸鹽污染控制提供依據(jù)。

1 研究區(qū)概況

羊莊盆地是大型山間構(gòu)造盆地，地形總體上東北高、西南低，高程為50～ 600 m[1]。盆地四周較封閉，出口狹窄，以巖溶網(wǎng)絡(luò)為介質(zhì)，具有輸水導(dǎo)水能力和一定調(diào)蓄功能，存在基本統(tǒng)一的水動(dòng)力場(chǎng)，屬于典型的封閉式泉排型巖溶水系統(tǒng)。羊莊盆地由辛召斷塊、山亭斷塊和羊莊斷塊3個(gè)水文地質(zhì)子單元組成[1，5-6]，受地形、地貌和地層巖性制約，盆地3個(gè)子單元間的水文地質(zhì)條件差異較大。辛召斷塊、山亭斷塊2個(gè)子單元位于盆地中上游，多為基巖裸露的低山、丘陵、山間盆地等，地勢(shì)高，地形起伏大，地下水資源較為貧乏[1]。位于盆地下游的羊莊斷塊子單元，斷裂構(gòu)造較多，巖溶發(fā)育程度高，富水性強(qiáng)，是整個(gè)盆地巖溶地下水循環(huán)徑流和富集儲(chǔ)存的主要場(chǎng)所，也是羊莊盆地水文地質(zhì)單元的主體部分。

由于羊莊盆地內(nèi)水質(zhì)取樣點(diǎn)集中在羊莊斷塊，且盆地內(nèi)的主要水源地集中在該斷塊，因此根據(jù)盆地地下水監(jiān)測(cè)點(diǎn)位、主要污染源以及主要水源地分布情況，最終確定研究范圍為羊莊斷塊水文地質(zhì)子單元。該子單元具有明確的邊界條件和完善的地下水補(bǔ)、徑、排、蓄條件，是一個(gè)相對(duì)獨(dú)立、較為完整的地下水子系統(tǒng)。大氣降水是該斷塊地下水的主要補(bǔ)給來源，總體徑流方向?yàn)橛蓶|北向西南，最終通過人工開采排泄。新薛河河谷穿越了整個(gè)羊莊斷塊，河水與地下水存在密切的互補(bǔ)互排關(guān)系[5，7]。

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

為了建立地下水流及溶質(zhì)運(yùn)移數(shù)學(xué)模型，對(duì)實(shí)際水文地質(zhì)條件進(jìn)行概化，建立水文地質(zhì)概念模型。根據(jù)實(shí)際水文地質(zhì)條件和研究需要，將含水系統(tǒng)概化為非均質(zhì)各向同性的二維非穩(wěn)定流[8]。模擬范圍為羊莊斷塊水文地質(zhì)子單元，總面積約134 km2，其邊界條件：西南及南部邊界為地表與地下水分水嶺相重合的隔水邊界：西部邊界中段和北段為化石溝隔水?dāng)嗔褞Ъ暗叵滤炙畮X組成的零流量邊界，南段為通過新薛河和第四系孔隙含水層向盆地外排泄的給定流量邊界[9]：東北和北部為以曹王墓斷裂為界的側(cè)向徑流補(bǔ)給邊界[9]；上邊界主要接受大氣降水、灌溉人滲補(bǔ)給以及新薛河補(bǔ)給和排泄：根據(jù)鉆孔資料統(tǒng)計(jì)，含水層大多集中在埋深220 m以上，其下巖溶發(fā)育微弱，因此將220 m埋深的巖層概化為隔水底板，即下邊界。

溶質(zhì)運(yùn)移模型主要考慮研究區(qū)農(nóng)業(yè)施肥和養(yǎng)殖場(chǎng)污水及糞便淋濾水人滲引起的地下水硝酸鹽污染，因此將研究區(qū)的養(yǎng)殖場(chǎng)處理為污染點(diǎn)源。另外，由于農(nóng)業(yè)生產(chǎn)施肥會(huì)導(dǎo)致地下水面源污染，因此在模型中根據(jù)已有農(nóng)田分布情況，處理為凈補(bǔ)給溶質(zhì)濃度。在溶質(zhì)運(yùn)移模型中，研究區(qū)北邊界處理為補(bǔ)給濃度邊界，東、西、南邊界均處理為零通量邊界。此次評(píng)價(jià)中蒸發(fā)濃度按0計(jì)算[10]。

3 數(shù)學(xué)模型始濃度；г1為給定濃度邊界；f（x，y，t）為г1邊界上的硝酸鹽濃度。

4 數(shù)值模擬

4.1 地下水流模型

利用Visual Modflow軟件進(jìn)行地下水流數(shù)值模擬計(jì)算。選擇2013年5月26日至2016年9月26日的水位觀測(cè)資料對(duì)模型進(jìn)行率定。水流模型經(jīng)過多次運(yùn)行、反復(fù)調(diào)整參數(shù)，最終使得觀測(cè)孔處計(jì)算水位和實(shí)測(cè)水位達(dá)到了最佳擬合[11]，見圖1。

最終研究區(qū)降水人滲補(bǔ)給系數(shù)劃分為6個(gè)分區(qū)，巖溶含水層滲透系數(shù)劃分為15個(gè)分區(qū)，分區(qū)情況見圖2、圖3，研究區(qū)1～6分區(qū)的大氣降水人滲補(bǔ)給系數(shù)分別為0.195、0.139、0.188、0.206、0.192、0.100，巖溶含水層水文地質(zhì)參數(shù)見表1。

4.2 溶質(zhì)運(yùn)移模型

4.2.1 彌散系數(shù)的確定

水動(dòng)力彌散系數(shù)是進(jìn)行溶質(zhì)運(yùn)移計(jì)算的最重要參數(shù)。為了準(zhǔn)確確定含水層的彌散系數(shù)，在羊莊盆地進(jìn)行了地下水彌散試驗(yàn)。試驗(yàn)將七鉬酸銨作為示蹤劑，投源孔選擇曹王墓斷裂附近的莊里30號(hào)孔，投入七鉬酸銨75 kg。在下游南塘、侯莊、羊莊水源地布置了3眼觀測(cè)孔進(jìn)行鉬酸根離子監(jiān)測(cè)。彌散試驗(yàn)投源井和觀測(cè)井位置見圖4。

試驗(yàn)采用逐點(diǎn)求參法求取水動(dòng)力彌散系數(shù)。取各觀測(cè)孔鉬酸根離子峰值時(shí)刻t1和各觀測(cè)孔鉬酸根離子檢出時(shí)刻t2及其對(duì)應(yīng)的濃度C1、C2，利用下式計(jì)算縱橫向水動(dòng)力彌散系數(shù)。式中：D_L為縱向彌散系數(shù)：t1為時(shí)刻1，即各觀測(cè)孔鉬酸根離子峰值時(shí)刻；t2為時(shí)刻2，即各觀測(cè)孔鉬酸根離子檢出時(shí)刻：x為觀測(cè)孔到投源孔的距離；C1為時(shí)刻tl時(shí)鉬酸根離子含量：C2為時(shí)刻t2時(shí)鉬酸根離子含量：u為地下水實(shí)際流速。

根據(jù)各觀測(cè)孔中示蹤劑濃度隨時(shí)間變化的監(jiān)測(cè)數(shù)據(jù)，繪制各觀測(cè)孔示蹤劑濃度C與監(jiān)測(cè)時(shí)間t的關(guān)系曲線，見圖5。地下水彌散系數(shù)計(jì)算參數(shù)見表2，地下水彌散度計(jì)算結(jié)果見表3。

4.2.2 模型硝酸鹽污染源強(qiáng)確定

研究區(qū)畜禽養(yǎng)殖場(chǎng)較多，大多分布在盆地周邊灰?guī)r裸露丘陵地區(qū)，主要養(yǎng)殖豬、雞、鴨。這些養(yǎng)殖場(chǎng)污水、糞便就地排放、堆積，均未采取防滲措施，污水及糞便淋濾水人滲導(dǎo)致地下水出現(xiàn)明顯的點(diǎn)狀、帶狀硝酸鹽污染。從水樣分析來看，以往及現(xiàn)在較大的養(yǎng)殖場(chǎng)分布地段或其附近下游一帶，地下水均出現(xiàn)了明顯的硝酸鹽超標(biāo)現(xiàn)象，養(yǎng)殖場(chǎng)規(guī)模越大、就地排泄物越多、持續(xù)時(shí)間越長(zhǎng)，超標(biāo)越嚴(yán)重。根據(jù)多年地下水污染源調(diào)查資料，取滲濾液硝酸鹽濃度的多年平均值60.0mg/L。鑒于畜禽養(yǎng)殖場(chǎng)常年污染地下水，因此將養(yǎng)殖場(chǎng)概化為連續(xù)補(bǔ)給的給定濃度點(diǎn)狀補(bǔ)給源賦值于模型中。

據(jù)調(diào)查，研究區(qū)農(nóng)田主要種植小麥和玉米，化肥主要施用尿素和復(fù)合肥，施用量為1 500～1 800 k～（ hm2.a），土壤中硝酸根離子在30 mg/L左右。模型計(jì)算中將其向地下水輸送的硝酸鹽濃度設(shè)定為30.0 mg/L，按耕地分布范圍以面狀補(bǔ)給的形式賦值于溶質(zhì)運(yùn)移模型中。

4.2.3 溶質(zhì)運(yùn)移模型的識(shí)別

依據(jù)研究區(qū)已有實(shí)測(cè)水質(zhì)數(shù)據(jù)，模擬期選2013年10月到2015年10月，其中以2014年10月硝酸鹽實(shí)測(cè)濃度場(chǎng)識(shí)別模型，以2015年10月硝酸鹽濃度場(chǎng)驗(yàn)證模型。經(jīng)過多次參數(shù)調(diào)整，最終使模型硝酸鹽計(jì)算濃度等值線與實(shí)測(cè)濃度等值線達(dá)到最佳擬合狀態(tài)，見圖6、圖7。最終確定的彌散度分區(qū)見圖8，彌散度共分為3個(gè)區(qū)，1～3區(qū)取值分別為34.2、12.3、5.6 m。

從圖6和圖7可以看出，研究區(qū)東南部硝酸鹽含量較高，并且呈放射狀向四周逐步遷移擴(kuò)散。計(jì)算與實(shí)測(cè)硝酸鹽濃度等值線的整體擬合程度良好，說明溶質(zhì)運(yùn)移模型中的參數(shù)取值是正確的，所建立的模型能較真實(shí)預(yù)測(cè)研究區(qū)地下水中硝酸鹽濃度場(chǎng)的變化規(guī)律[12]。

5 研究區(qū)地下水硝酸鹽演化趨勢(shì)預(yù)測(cè)

用識(shí)別驗(yàn)證后的模型，現(xiàn)狀各源匯項(xiàng)不變的情況下，預(yù)測(cè)5a及10 a后的硝酸鹽濃度等值線，見圖9～圖11。預(yù)測(cè)結(jié)果表明：在預(yù)測(cè)時(shí)段，地下水中硝酸鹽濃度呈明顯上升趨勢(shì)，模擬區(qū)硝酸鹽濃度較高的地區(qū)主要分布在轆轆井區(qū)域，在地下水動(dòng)力場(chǎng)影響下，硝酸鹽高濃度區(qū)域主要向西南方向擴(kuò)展。

6 結(jié)語

基于Visual Modflow軟件建立了地下水流與溶質(zhì)運(yùn)移數(shù)值模型，對(duì)羊莊斷塊地下水硝酸鹽的遷移與擴(kuò)散進(jìn)行了數(shù)值模擬，預(yù)測(cè)了10 a后硝酸鹽在地下水中的分布狀況以及濃度變化趨勢(shì)。

由硝酸鹽濃度預(yù)測(cè)結(jié)果可知，若不采取措施，羊莊盆地地下水中硝酸鹽濃度將呈增大趨勢(shì)，這主要是當(dāng)?shù)剞r(nóng)民對(duì)氮肥的過量施用和盆地內(nèi)養(yǎng)殖場(chǎng)糞便隨意堆放造成的。因此建議采取必要的措施控制農(nóng)業(yè)施肥的面狀污染源，并對(duì)以往養(yǎng)殖場(chǎng)發(fā)生過污染的場(chǎng)地進(jìn)行土壤和地下水污染修復(fù)，以削減殘留污染物對(duì)地下水可能造成的二次污染。根據(jù)預(yù)測(cè)結(jié)果可知，轆轆井區(qū)域地下水中的硝酸鹽長(zhǎng)期存在濃度較高的現(xiàn)象，主要原因是該處以往建有多個(gè)規(guī)模較大的養(yǎng)殖場(chǎng)，雖然這些對(duì)地下水產(chǎn)生污染的養(yǎng)殖場(chǎng)已于2010年前后被清理拆除，但當(dāng)時(shí)未進(jìn)行必要的土壤修復(fù)，而完全靠含水層自身凈化修復(fù)需要一個(gè)很長(zhǎng)的時(shí)間過程。為此，建議在地下水受到明顯污染的地區(qū)，盡快采用人工修復(fù)方式降低污染組分的含量。

參考文獻(xiàn)：

[1] 蘇秋實(shí).山東省羊莊水源地地下水資源預(yù)測(cè)及開發(fā)利用研兗[D].濟(jì)南：山東師范大學(xué)，2009：1-2.

[2] 張新民.低碳農(nóng)業(yè)的外部性與市場(chǎng)失靈[J].天津農(nóng)業(yè)科學(xué)，2012，18（2）：61-64.

[3] 張翠云，馬琳娜，張勝，等.Visual Modflow在石家莊市地下水硝酸鹽污染模擬中的應(yīng)用[J].地球?qū)W報(bào)，2007，28（6）：561-566.

[4] 孫志浩，王開章，孔凡亮.基于FEFLOW的淺層地下水水環(huán)境演化模擬[J].中國(guó)農(nóng)村水利水電，2009（ 10）：5 -7，11.

[5] 卜華，陳占成，張良鵬.飲用水水源地保護(hù)區(qū)劃分研究：以山東羊莊盆地地下水水源地為例[J].地質(zhì)調(diào)查與研究，2008，31（3）：236-241.

[6] 曲英森，郭留福，馮振中，等.山東省棗莊市羊莊盆地水文地質(zhì)初步勘探報(bào)告[R].兗州：山東省地質(zhì)局第三水文地質(zhì)工程地質(zhì)大隊(duì).1981：31-33.

[7] 吳愛民，李傳謨，文唐章，等.山東省羊莊盆地巖溶水均衡試驗(yàn)報(bào)告[R].兗州：山東省地礦局第三水文地質(zhì)工程地質(zhì)大隊(duì).1988：50-51.

[8] 張俊華，于福榮，李志萍.FEFLOW在黃龍T業(yè)園地下水?dāng)?shù)值模擬中的應(yīng)用[J].人民黃河，2011，33（6）：69-70.

[9] 張秋蘭，崔亞莉，邵景力，等.基于數(shù)值模型的山東省羊莊水源地保護(hù)區(qū)劃分淺析[J].水文地質(zhì)T程地質(zhì)，2008，35（增刊1）：423 -427.

[10] 郝蓉.巴彥淖爾市狼山山前沖洪積扇水環(huán)境質(zhì)量評(píng)價(jià)與溶質(zhì)運(yùn)移模擬的研究[D].呼和浩特：內(nèi)蒙古農(nóng)業(yè)大學(xué)，2010：48-50.

[11] 馮雯娟，劉國(guó)東，甘蓉，等.基于GMS的巖溶區(qū)火電廠地下水污染物運(yùn)移模擬[J].人民黃河，2015，37（3）：89-91.

[12] 崔欣.鎮(zhèn)賚縣生態(tài)環(huán)境脆弱區(qū)引水灌溉對(duì)地下水及環(huán)境的影響研究[D].長(zhǎng)春：吉林大學(xué)，2008：68 - 69.