王喜峰，賈仰文，牛存穩(wěn)

(1.中國社會科學(xué)院數(shù)量經(jīng)濟與技術(shù)經(jīng)濟研究所，北京 100732;2.中國水利水電科學(xué)研究院，北京 100038)

隨著我國工業(yè)化和城鎮(zhèn)化的不斷發(fā)展，水資源利用量也在不斷增加，甚至超過了區(qū)域水資源承載能力，造成了嚴(yán)重的水危機。同時，大量的污染物隨著廢水排入河道、湖泊以及地下含水層中，導(dǎo)致了水環(huán)境惡化，危害到可持續(xù)發(fā)展、人類健康和生活環(huán)境。

開采地下水是許多流域和地區(qū)最為重要的應(yīng)對水危機的方式之一。然而，由于過度的農(nóng)業(yè)和工業(yè)活動，許多地區(qū)的地下水環(huán)境已不容樂觀。氮是地下水常見污染物，主要來自于農(nóng)業(yè)的面源污染和工業(yè)污水排放到河流形成的線源污染。因此，對地下含水層氮污染負(fù)荷的研究對地下水環(huán)境管理有著重要意義。

為了應(yīng)對嚴(yán)重的水危機，大量的研究關(guān)注于對地下含水層氮污染物來源的辨識和測算方面［1－2］，許多學(xué)者還進(jìn)行了地下水氮污染的研究［3－6］，但是對地下含水層氮污染負(fù)荷來源的研究還不是很多。許多研究只對不能被作物吸收的氮進(jìn)行研究，以提高肥料的利用率［5］。然而，對于進(jìn)入地下含水層中的氮肥比例，考慮到污染物運移和反應(yīng)的機制等因素，導(dǎo)致相關(guān)方面的研究不足以對地下水的水質(zhì)做出判斷［7－10］。

此外，研究工具也有相對局限性。由于尺度效應(yīng)的存在，對相關(guān)區(qū)域的實驗不能滿足對污染物來源的研究。然而，對地下水污染負(fù)荷的影響因素主要還是由水流和污染物的量決定的。在對地下含水層污染物來源的研究中，污染物的量由每個城市(區(qū)域)的污染物使用量獲得，水流情況由水循環(huán)模擬得到。因此，水循環(huán)模型應(yīng)該給出相對準(zhǔn)確的水流數(shù)據(jù)，包括地表、土壤層和含水層的相互作用數(shù)據(jù)。因為污染物在土壤層的復(fù)雜反應(yīng)，因此檢測數(shù)據(jù)應(yīng)該取代土壤層中污染物反應(yīng)后的數(shù)據(jù)。本文建立了一個地下水環(huán)境綜合模擬框架(IGESF)。該框架包括地下水?dāng)?shù)值模擬、污染物運移模擬、分布式水文模型、水環(huán)境模擬以及檢測數(shù)據(jù)。

1 地下水環(huán)境綜合模擬框架開發(fā)與驗證

1.1 地下水綜合模擬框架簡介

如圖1所示，大尺度流域地下水環(huán)境綜合模擬框架是由4個模塊綜合模擬而成。在來自河道線源地下水污染負(fù)荷模擬時，第1模塊“二元”水循環(huán)模型(WEP-L)為第2模塊地表水系環(huán)境模型提供坡面匯流和河道匯流動力條件;第1模塊WEP-L向第3模塊MODFLOW提供河道水位模擬結(jié)果模擬河道進(jìn)入到地下水的水量;第2模塊地表水系環(huán)境模型則模擬各河段的污染物濃度;根據(jù)模擬的濃度和河道補給地下水的水量計算污染物線源污染負(fù)荷。在研究地下水污染物面源負(fù)荷時，通過分析第1模塊WEP-L模擬的面狀地下水分布式補給量及埋深處實測污染物濃度場得到地下水分布式面源污染負(fù)荷。將地下水線源和面源負(fù)荷輸入第4模塊地下水污染物運移模型中，模擬污染物在地下水含水層通過對流、彌散稀釋、吸附和消去等作用后的污染物濃度情況。

1.2 IGESF各個模塊介紹

1)大尺度空間插值方法。大尺度空間插值方法是模型間數(shù)據(jù)交換的主要工具，本文應(yīng)用王喜峰等［11］開發(fā)的改進(jìn)的RDS方法。

2)“二元”水循環(huán)模型是基于物理機制的分布式水文模型WEP-L。WEP-L(water and energy transfer processes in large river basins)模型是建立在WEP模型基礎(chǔ)之上，以適應(yīng)水循環(huán)模型在大尺度流域上的應(yīng)用。

3)地表水系環(huán)境模型(WEQ)是建立在分布式水文模型(WEP-L)的基礎(chǔ)上，對流域面源和點源產(chǎn)生的氮的產(chǎn)生量和入河量進(jìn)行了模擬。該模型由牛存穩(wěn)［1］開發(fā)。

4)MODFLOW是modular three-dimensional finite-difference ground-water flow model(三維有限差分地下水流數(shù)值模擬模型)的簡稱，由美國地質(zhì)調(diào)查局(USGS)McDonald和Harbaugh于20世紀(jì)80年代開發(fā)的。

5)MT3DMS是應(yīng)用最廣泛的地下水溶質(zhì)運移模型，由鄭春苗教授等于1999年開發(fā)，是基于地下水污染物運移模型MT3D的第二代地下含水層溶質(zhì)運移模型。

圖1 地下水環(huán)境綜合模擬框架示意圖

1.3 IGESF的驗證

首先，進(jìn)行分布式水循環(huán)模型(WEP-L)的開發(fā)和驗證。華北平原的分布式水循環(huán)模型是海河流域分布式水循環(huán)模型的一部分，由國家基礎(chǔ)研究發(fā)展規(guī)劃項目“海河流域水循環(huán)演變機理與水資源高效利用”項目組開發(fā)驗證。模擬結(jié)果 (圖2和表1)表明，模型模擬效果較好，基本反映了華北平原水循環(huán)狀況，為地下水水流模擬提供了合理的地下水循環(huán)情況。

其次，河道水環(huán)境模擬結(jié)果較為成功，這在牛存穩(wěn)［1］的研究中很好地反映出來。河道水環(huán)境模擬對河道與地下含水層中的相互作用有很好的支持作用。

再次，在陸垂裕博士研究的基礎(chǔ)上，建立了華北平原地下水水流模型MODFLOW。與之前的研究不同，本次研究將河道模擬以及河道與地下水的水流相互作用反映出來。這個主要模擬活動包括降雨、從河道以及地下水取用灌溉過程、渠道滲漏、蒸散發(fā)、地下水利用、河道滲漏等?？紤]到華北平原面積大約為130000 km2，MODFLOW 模型采用了4 km×4 km的格子進(jìn)行模擬，設(shè)置了上下兩層(淺層地下水、深層地下水)共16100個格子。邊界條件和初始條件見圖3和圖4。

表1 華北平原水循環(huán)模擬開發(fā)與驗證

圖2 水循環(huán)模擬狀況

圖3 華北平原網(wǎng)格和邊界條件

模擬期限為1996—2000年。模擬為月度模擬，所以共有60個模擬時段，采用的步長為天，水流的初始條件為1995年12月31日的水場情況。由于模擬框架需要河道和地下含水層的相互作用情況，因此河道也設(shè)置在MODFLOW中，如果模擬網(wǎng)格有河道，將該網(wǎng)格設(shè)置為河道單元，并且計算河道與地下含水層的相互作用。河道的參數(shù)化情況見圖5。

圖6和圖7為華北平原MODFLOW的模擬結(jié)果，可見淺層地下水和深層地下水水場模擬結(jié)果與實測結(jié)果相似。

圖4 華北平原初始條件

圖5 華北平原河道參數(shù)化情況

圖6 淺層地下水的實測與模擬對比

圖7 深層地下水模擬與實測情況對比

最后，運用MT3DMS模擬了氮在地下含水層的運移情況。地下水水流模擬由MODFLOW提供，污染物運移的主要參數(shù)參考文獻(xiàn)［6，12，14］。采用的主要反應(yīng)模型是一階不可逆反應(yīng)形式，同時模型綜合考慮了吸附項和彌散項。表2給出了主要的參數(shù)，圖8給出了主要的氮污染物濃度的模擬和實測情況。

表2 華北平原MT3DMS主要的參數(shù)情況

圖8 地下含水層氮污染濃度模擬與實測對比

2 結(jié)果分析

2.1 地下含水層氮污染負(fù)荷的來源分析

根據(jù)1995—2004年IGESF和華北平原的實測數(shù)據(jù)，山前平原地區(qū)(565 t/16 km2)和唐山海濱區(qū)域(739 t/16 km2)的地下含水層氮污染物負(fù)荷比平原中部地區(qū)大得多(見圖9)，這與工農(nóng)業(yè)的實際情況有很大的相關(guān)性。如圖10所示，來自面源的地下含水層氮污染負(fù)荷大約是38400 t/年。特別是在1996年，來自面源的地下水氮污染負(fù)荷達(dá)到98400 t。1996年進(jìn)入到地下含水層中的來自面源污染的氮污染負(fù)荷大于面源進(jìn)入到河道的氮污染負(fù)荷(81200 t)。但是在2002年，只有9900 t來自面源的氮污染負(fù)荷進(jìn)入地下含水層中。

總的來說，1995—2004年每年大約有38400 t來自面源的氮污染負(fù)荷進(jìn)入到地下含水層中。在干旱年份，只有氮肥施用量1%的氮隨著較少的入滲水流進(jìn)入到地下含水層中;在雨量充沛年份，大約有氮肥施用量5%的氮進(jìn)入到地下含水層中。

圖9 1995—2004年來自面源的地下含水層氮污染負(fù)荷

圖10 歷年來自面源的地下含水層氮污染負(fù)荷

根據(jù)圖11可以看出，每年由河道進(jìn)入到地下含水層中的氮大約為26000 t(1995—2004年)，是每年進(jìn)入到河道中的氮污染物負(fù)荷的30%。根據(jù)模擬結(jié)果，城市的下游河道是氮污染物進(jìn)入地下含水層的集中區(qū)域，特別是在大城市的下游地區(qū)，例如北京、天津的下游地區(qū)都有大量的氮污染負(fù)荷經(jīng)河道進(jìn)入地下含水層中，這與實際情況是相符的。根據(jù)圖12，在2000年之后，經(jīng)河道進(jìn)入地下含水層中的氮污染負(fù)荷顯著上升。另外北京、天津和濮陽應(yīng)降低河道的氮污染情況，從現(xiàn)有情況來看，這些區(qū)域含水層氮污染情況不容樂觀。

圖11 1995—2004年經(jīng)河道進(jìn)入地下含水層的氮污染負(fù)荷

圖12 經(jīng)河道進(jìn)入地下含水層中的氮污染負(fù)荷情況

2.2 氮污染負(fù)荷的影響因素

本節(jié)通過研究氮肥使用量與地下含水層氮污染負(fù)荷的關(guān)系發(fā)現(xiàn):盡管土壤結(jié)構(gòu)、種植結(jié)構(gòu)以及包氣帶情況差異巨大，但是，地下含水層氮污染負(fù)荷與氮肥使用量有很大的線性相關(guān)性(相關(guān)系數(shù)達(dá)到0.61)，如圖13所示。因此，氮肥的使用是地下含水層氮污染負(fù)荷的主要影響因素。

圖13 氮肥使用與地下含水層氮污染負(fù)荷相關(guān)性

此外，本部分也測算了地下含水層氮污染負(fù)荷與降雨和灌溉用水之和的相關(guān)性。這主要考慮到降雨和灌溉是氮污染進(jìn)入到地下含水層中的主要且必要的動力條件。根據(jù)IGESF的模擬結(jié)果，本文發(fā)現(xiàn)降雨和灌溉用水的和與地下含水層氮污染負(fù)荷有很強的相關(guān)性，且呈現(xiàn)指數(shù)曲線關(guān)系，如表3和圖14所示。

表3 各地市地下含水層氮污染負(fù)荷與降雨和灌溉量和的相關(guān)性

圖14 各地市地下含水層氮污染負(fù)荷與降雨和灌溉和的關(guān)系

3 結(jié)束語

考慮到現(xiàn)有研究的不足以及地下水環(huán)境管理的實際需要，本文以華北平原為例，提出了測度地下含水層氮污染負(fù)荷來源的方法。建立的地下水環(huán)境綜合模擬框架(IGESF)不僅可以模擬氮污染負(fù)荷，對其他污染負(fù)荷的研究也有一定的借鑒作用。地下水環(huán)境綜合模擬框架綜合了分布式水文模型、水環(huán)境模型、地下水?dāng)?shù)值模擬、污染物運移模型以及實測數(shù)據(jù)。根據(jù)該方法，1995—2004年大約有38400 t/年的氮污染負(fù)荷經(jīng)面源進(jìn)入地下含水層中，大約26000 t/年的氮污染負(fù)荷經(jīng)河道進(jìn)入到地下含水層中。從空間來看，唐山濱海地區(qū)和北京、石家莊等山前地區(qū)比中部地區(qū)的氮污染負(fù)荷大。從影響地下水氮污染負(fù)荷的影響因素來看，氮肥使用、降水和灌溉是影響氮污染負(fù)荷的主要因素，應(yīng)當(dāng)加強節(jié)水和節(jié)肥等相關(guān)對策措施以降低地下含水層的氮污染負(fù)荷。

［1］Niu Cunwen.Integrated Modeling and Application of Water Resources and Water Environmental in Basin Scale［M］.Beijing:China Institute of Water Resources and Hydropower Research，2008.

［2］Nyborg M，Malhi S S，Solberg E D.Carbon Storage and Light Fraction C in a Grassland Dark Grray Chemical Soil as Influenced by N and S fertilization［J］.Can J Soil Svi，1990，79:317－320.

［3］Zhang Weiying.Survey about Nitrate pollution of groundwater by usage of nitrogenous fertilizer in North area of China［J］.Plant nutrition and fertilizer sciences，1995(2):80－87.

［4］Bouwer H.Effect of Irrigated Agriculture on Groundwate［J］.Journal of Irrigation and Drainage Engineering，1987，113(1):4－15.

［5］Liu Hongbin.Li Zhihong.Characteristics of Nitrate Distribution and Accumulation in Soil Profiles under Main Agro-land Use Types in Beijing［J］.China Agricultural Sciences，2004，39(5):1－6.

［6］Zheng Chunmiao，Liu Jie，Cao Guoliang.Research on the Groundwater Quality of North China Plain.A Special Research in National Key Basic Research Development Program(973):Water cycle Evolution Mechanism in Haihe River Basin and efficient use of water resources［M］.Beijing:Peking University，2010.

［7］Liu X X，Ju F，Zhang J.Nitrogen dynamics and Bughets in a Winter Wheat-maize Cropping System in the North China Plain［J］.Field Crops Research，2003(2):111－124.

［8］Shen Shanmin.The contribution of Nitrogenous fertilizer to China’s Agricultural Development and Losses of nitrogen in agriculture［J］.SOIL，2002(39):12－24.

［9］Goderya D J，Dahab M F.Woldt W E.Incorporation of Spatial Variability in Modeling Non-point Source Groundwater Nitrate Pollution［J］.Water Science and Technology，1996，4(5):233－240.

［10］Li Shiqing，Li Shengxiu.Nitrate leaching in farmland ecosystem of semi-arid areas［J］.Journal of Applied Ecology，2000，11(2):240－243.

［11］Wang X，Zhou Z，Jia Y，et al.Amendatory Reversed Distance Square method to interpolate precipitation of long series in large-scale basins［C］//Proceedings of the 9th Imternational Conference on Hydroinformatics.2010:322－330.

［12］Lu G.Clement C.Natural attenuation of BTEX compounds，model development and field-scale application［J］.Groundwater，1999，37(5):707－717.

［13］Gelhar L W，Welty K W.A Critical Review of Data on Field-scale Dispersion in Aquifers［J］.Water Resources Research，1992，28(7):1955－1974.

［14］Jia Yangwen，Ni Guangheng.Development of WEP model and its application to an urban watershed［J］.Hydrology Process，2001，15(11):2175－2194.

［15］Wang Xifeng.Impacts of Climate Change on Groundwater and Groundwater Environment of Large-scale Region［M］.Beijing:China Institute of Water Resources and Hydropower Research，2011.