韓博,周麗麗,范昊明,黃東浩
(沈陽農(nóng)業(yè)大學(xué)水利學(xué)院,110866,沈陽)
非點(diǎn)源污染(面源污染)是指溶解性或固體污染物在降水和徑流沖刷作用下匯入受納水體而引起的水體污染,其主要來源包括土壤侵蝕、化肥與農(nóng)藥的使用過量、城市和公路徑流、畜禽養(yǎng)殖和農(nóng)業(yè)與農(nóng)村廢棄物等。近年來,隨著農(nóng)業(yè)活動(dòng)的增加和生態(tài)植被的破壞,由農(nóng)藥、化肥、耕作等產(chǎn)生的農(nóng)業(yè)非點(diǎn)源污染問題日益突出,尤其表現(xiàn)在水資源的污染與破壞。土壤養(yǎng)分吸附于泥沙進(jìn)入水體或溶解在徑流中流失,都會(huì)造成嚴(yán)重的水體非點(diǎn)源污染[1]。氮是作物生長必要的營養(yǎng)因子,同時(shí)也是引起水體富營養(yǎng)化和水質(zhì)變差的根源因子,這已引起全世界的高度關(guān)注。資料顯示,歐洲由于農(nóng)業(yè)活動(dòng)所造成的大量氮污染達(dá)到入海通量的60%[2],在瑞典最南端的謝夫靈厄流域,由于農(nóng)業(yè)活動(dòng)而造成的河流運(yùn)輸?shù)空剂饔蚩偟康?4% ~87%[3],荷蘭由農(nóng)業(yè)非點(diǎn)源污染所提供的總氮量也占到水環(huán)境污染總量的60%之多[4],在我國南方的三湖(巢湖、滇池和太湖)地區(qū)湖泊富營養(yǎng)化的現(xiàn)象日趨嚴(yán)重,其由農(nóng)業(yè)活動(dòng)所造成的總氮污染貢獻(xiàn)率達(dá)到60% ~70%[5]。郭紅巖等[6]對太湖流域非點(diǎn)源氮污染下水質(zhì)影響進(jìn)行定量化研究,研究結(jié)果顯示:農(nóng)田氮排放量占總量的72.7%,同樣在我國北方的吉林省雙陽水庫,每年僅氮的負(fù)荷量也達(dá)到389.4t[7]。由此可見,對流域氮素污染負(fù)荷的研究,可為日后治理流域非點(diǎn)源污染提供規(guī)劃依據(jù)。
SWAT模型是由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心(USDA,Agricultural Research Service)開發(fā)的分布式流域水文模型,目前在農(nóng)業(yè)非點(diǎn)源污染研究中得到了廣泛的應(yīng)用,模型已被應(yīng)用于流域徑流模擬、流量預(yù)測和非點(diǎn)源污染評價(jià)等方面,并取得了顯著研究成果。萬超等[8]應(yīng)用SWAT模型對潘家口水庫上游農(nóng)業(yè)非點(diǎn)源污染負(fù)荷進(jìn)行計(jì)算,研究以不同化肥施用水平進(jìn)行模擬,模擬結(jié)果認(rèn)為合理施用化肥對減少非點(diǎn)源污染有重大意義。王中根等[9]應(yīng)用SWAT模型建立海河流域分布式水文模型,研究中對土地利用方式和土壤類型進(jìn)行了重分類處理,并應(yīng)用多個(gè)站點(diǎn)的實(shí)測數(shù)據(jù)進(jìn)行結(jié)果校正與驗(yàn)證,其研究成果為SWAT模型在國內(nèi)大型流域的研究應(yīng)用提供范例。王林等[10]應(yīng)用SWAT模型對東南沿海的晉江流域建立水文模型,通過年、月徑流的模擬分析可知,SWAT模型適用于降水量豐沛的東南沿海濕潤地區(qū)并且模擬效率受降水量及徑流調(diào)節(jié)影響較大。謝元博[11]在涇河流域建立非點(diǎn)源污染模型,通過模擬分析可知污染負(fù)荷的產(chǎn)輸出與當(dāng)年的降雨徑流量有密切關(guān)系。應(yīng)用SWAT模型進(jìn)行流域模擬及分析已在世界范圍內(nèi)得到了有效證實(shí),但在我國東北三省地區(qū),針對流域面積不足3 000 km2的中尺度流域農(nóng)業(yè)非點(diǎn)源氮素污染研究實(shí)屬罕見;因此,筆者以蒲河流域?yàn)檠芯繀^(qū)域,應(yīng)用基于GIS的SWAT模型,進(jìn)行蒲河流域氮素模擬并對其進(jìn)行時(shí)空分布研究,為蒲河流域的生態(tài)恢復(fù)研究提供依據(jù)。
蒲河位于遼寧省沈陽市中部的產(chǎn)糧區(qū)(E 122°40′48″~123°56′32″,N 41°21′53″~ 42°4′15″),發(fā)源于鐵嶺縣橫道河子鄉(xiāng),在棋盤山風(fēng)景區(qū)上游望濱鄉(xiāng)石砬子村進(jìn)入沈陽市轄區(qū),流經(jīng)沈北新區(qū)、東陵、于洪、新民市和遼中5個(gè)縣(市)區(qū),于遼中縣老觀坨鄉(xiāng)黑魚溝村入渾河。河道全長205 km,流域面積2 248 km2。蒲河較大支流有九龍河和小渾河,主要排干有沈新遼、于家臺、烏伯牛等。蒲河流域?qū)贉貛Т箨懶约撅L(fēng)氣候,多年平均氣溫8.2℃。多年平均降水量647.4 mm,降水主要集中在7、8月,2個(gè)月降水量占年降水量的50%左右。蒲河流域上游為低山區(qū),占流域面積的5.92%,區(qū)內(nèi)山巒重疊,山高坡陡,但區(qū)內(nèi)植被豐富,生長較好。中上游為丘陵區(qū),占流域面積的13.48%,區(qū)內(nèi)丘陵起伏,植被覆蓋度較低;中下游為平原區(qū),占流域面積的80.6%。該區(qū)地勢平坦,土壤肥沃,是沈陽市的主要糧食產(chǎn)區(qū)。土壤類型主要有棕壤、草甸土、水稻土、風(fēng)沙土和沼澤土等5類。
SWAT模型是在20世紀(jì)后期由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心推出的分布式流域水文模型,能夠?qū)哂袕?fù)雜下墊面條件、多種土地利用方式和不同管理措施下的流域徑流、泥沙及化學(xué)物質(zhì)的產(chǎn)輸出進(jìn)行模擬計(jì)算。SWAT模型內(nèi)主要含有水文過程子模型、土壤侵蝕子模型和污染負(fù)荷子模型3部分。
SWAT模型中水文過程子模型采用的水量平衡表達(dá)式為
式中:St為土壤最終含水量,mm;S0為土壤初始含水量,mm;t為時(shí)間步長,d;Rday為第 i天降水量,mm;Qsurf為第 i天地表徑流量,mm;Ea為是第 i天蒸發(fā)量,mm;Wseep為第i天存在于土壤剖面地層的滲透量和側(cè)流量,mm;Qgw為第i天地下水含量,mm。
SWAT模型中土壤侵蝕子模型采用的修正后通用土壤流失方程為:
式中:msed為土壤流失量,t;qpeak為洪峰徑流,m3/s;Ahru為水文響應(yīng)單元(HRU)的面積,hm2;KUSLE為土壤侵蝕因子;CUSLE為植被覆蓋和管理因子;PUSLE為保持措施因子;LUSLE為坡長因子;SUSLE為坡度因子;WCFRG為粗碎屑因子。
SWAT模型能夠追蹤流域內(nèi)氮磷隨徑流和泥沙的遷移和轉(zhuǎn)化,土壤中氮磷形態(tài)轉(zhuǎn)化是由氮循環(huán)和磷循環(huán)來控制實(shí)現(xiàn)的。
模型數(shù)據(jù)庫構(gòu)建按應(yīng)用要求分為空間數(shù)據(jù)庫和屬性數(shù)據(jù)庫2大類??臻g數(shù)據(jù)庫依托于GIS的支持,建立包括數(shù)字高程模型(DEM)、土壤類型、土地利用等空間數(shù)字圖層。屬性數(shù)據(jù)庫則涵蓋研究區(qū)水文、氣象、水庫、土壤屬性、農(nóng)業(yè)管理措施等文本數(shù)據(jù),其中:研究區(qū)DEM格式為“.GeoTIFF”,來源于ASTER GEDM數(shù)據(jù);土壤類型圖格式為“.Grid”,來源于中國科學(xué)院南京土壤研究所與中國農(nóng)業(yè)部共同推出1∶100萬數(shù)字化土壤圖;土地利用類型圖格式為“.Grid”,來源于地球系統(tǒng)科學(xué)數(shù)據(jù)共享平臺;水文、氣象和水庫數(shù)據(jù)格式為“.xls”,來源于國家氣象信息中心和遼河流域水文年鑒;研究區(qū)土壤屬性和農(nóng)業(yè)管理措施數(shù)據(jù)格式為“.xls”,來源于遼寧土種志、遼寧土壤和遼寧省統(tǒng)計(jì)年鑒。根據(jù)SWAT模型需求,將經(jīng)緯度坐標(biāo)轉(zhuǎn)化為大地坐標(biāo)系,為便于對研究區(qū)域水文循環(huán)特征的描述,研究選用保留面積性質(zhì)的Albers等積圓錐投影作為投影方法,并將空間數(shù)字圖層坐標(biāo)系進(jìn)行統(tǒng)一設(shè)置。
2.3.1 模型應(yīng)用過程 應(yīng)用蒲河流域數(shù)字高程模型(DEM)在SWAT模型中進(jìn)行河網(wǎng)提取并利用實(shí)際水系對其進(jìn)行校準(zhǔn)修正,將蒲河流域內(nèi)老觀陀鄉(xiāng)黑魚溝村所在位置作為總流域出口進(jìn)行模擬,同時(shí)將蒲河流域內(nèi)的中型水庫——棋盤山水庫和各雨量站所在的位置作為子流域出口點(diǎn),以方便水庫定位和校準(zhǔn)。經(jīng)過模擬運(yùn)行,將蒲河流域劃分為23個(gè)子流域,劃分子流域如圖1所示。將土地利用類型圖和土壤類型圖及其模型所需的土地類型查詢表和土壤類型查詢表載入模型中。模型采用優(yōu)勢土壤(或優(yōu)勢土地利用)法進(jìn)行HRU劃分,通過設(shè)定閾值在子流域里劃分出多個(gè)水文響應(yīng)單元,模型結(jié)合蒲河流域?qū)嶋H情況,確定土地利用和土壤面積閾值均為4%,將蒲河流域劃分為215個(gè)HRU。將各氣象信息測站表導(dǎo)入模型中,模型則根據(jù)測站位置讀取流域各氣象信息(包括日降水?dāng)?shù)據(jù)、日最高氣溫、日最低氣溫、風(fēng)速、相對濕度、太陽輻射等)。在模型中加入水庫位置及其水庫屬性數(shù)據(jù)。根據(jù)流域土層實(shí)際狀況及查詢遼寧土種志和遼寧土壤,在模型中填入土層氮磷初始濃度。結(jié)合流域農(nóng)業(yè)管理實(shí)際狀況及遼寧省統(tǒng)計(jì)年鑒,在模型中輸入農(nóng)業(yè)管理措施數(shù)據(jù)。運(yùn)行模型并利用實(shí)測數(shù)據(jù)對模型參數(shù)校準(zhǔn)修正,提高模型精度,最終達(dá)到模型適宜性評價(jià)效果并應(yīng)用校準(zhǔn)模型對研究區(qū)氮污染負(fù)荷進(jìn)行時(shí)空特征分析。
圖1 蒲河流域子流域劃分Fig.1 Sub-basin plans in Puhe river basin
2.3.2 模型參數(shù)校準(zhǔn) 模型參數(shù)校準(zhǔn)采用自動(dòng)與手動(dòng)調(diào)參相結(jié)合的方式,通過多次模擬運(yùn)行,使模擬值與實(shí)測值擬合較好,實(shí)現(xiàn)參數(shù)最優(yōu)化取值。模型校準(zhǔn)按照先率定徑流量再率定水質(zhì)的順序操作。利用多年實(shí)測徑流資料對模型徑流參數(shù)進(jìn)行敏感性分析,率定分析結(jié)果顯示:對模型輸出影響較大的包括影響地下水與地表徑流之間交換關(guān)系的基流α系數(shù)(ALPHA_BF),影響徑流過程的SCS徑流曲線系數(shù)(CN2)和地表徑流滯后時(shí)間(SURLAG),影響蒸發(fā)過程的土壤蒸發(fā)補(bǔ)償系數(shù)(Esco),影響土壤水分過程的土壤深度(SOL_Z),影響水層滲透的深蓄水層滲透系數(shù)(Rchrg_Dp),影響河道匯流過程的河道有效水力傳導(dǎo)系數(shù)(CH_K2)和主河道曼寧系數(shù)值(CH_N),影響降雪和融雪過程的6月21日融雪系數(shù)(SMFMX)和12月21日融雪系數(shù)(SMFMN)。通過查閱相關(guān)文獻(xiàn)和資料結(jié)合研究區(qū)實(shí)際狀況,確定影響水質(zhì)的參數(shù)為:氮滲透系數(shù)(NPERCO),磷滲透系數(shù)(PPERCO)[12]。模型率定參數(shù)見表1。
表1 模型的率定參數(shù)Tab.1 Calibration parameters of SWAT model
應(yīng)用實(shí)測數(shù)據(jù)對模擬數(shù)值進(jìn)行校準(zhǔn)擬合:校準(zhǔn)期徑流相對誤差Re=-0.117 5,在±20%范圍以內(nèi);R2=0.719 0>0.6;Nash-Sutcliffe效率系數(shù)Ens=0.658 0>0.5。驗(yàn)證期Re=-0.068 9,在 ±20%范圍內(nèi);R2=0.882 0>0.6;Ens=0.856 8>0.5。模擬效果較好,總體反應(yīng)了流域的實(shí)際情況。SWAT模型能夠追蹤流域內(nèi)氮磷隨徑流和泥沙的遷移和轉(zhuǎn)化,因此,在徑流量擬合較好的前提下,氮的污染模擬也具有一定的可信性和參考價(jià)值。
表2 蒲河流域監(jiān)測斷面模擬結(jié)果輸出Tab.2 Simulation results output of Puhe river basin monitoring sections
3.1.2 月際分布規(guī)律 通過對多年徑流量與模擬氮污染負(fù)荷量分析可知,徑流與TN有較好的相關(guān)性,R2=0.813 6。選取枯水年(2000年)、平水年(2004年)、豐水年(2010年)為代表年份,對不同水文年各月TN量進(jìn)行模擬分析,TN模擬結(jié)果月輸出值如表3所示,不同水文年徑流與TN月擬合圖如圖2所示。分析圖表可知:TN量在年內(nèi)各月呈波動(dòng)變化,1—4月的TN量較低,5—6月開始波動(dòng),7—9月達(dá)到最大,10—12月有所回落。豐水年(2010年),降水量多,徑流量大,徑流峰值出現(xiàn)在8月份,TN量也因降水豐沛而呈增加趨勢,并在9月份達(dá)到含量高峰;平水年(2004年),徑流峰值出現(xiàn)在汛期,在7月份達(dá)到最大值,TN量自8月份增長顯著,9月份出現(xiàn)增長高峰;枯水年(2000年),由于受1999年枯水年降水影響,加之該年降水量也極少,年內(nèi)徑流峰值出現(xiàn)在前半年,而TN卻在8月份達(dá)到含量峰值。由此可見,TN流失具有明顯的滯后性。
表3 TN模擬結(jié)果月輸出Tab.3 Month results output of TN simulation103t/month
依研究區(qū)分布特點(diǎn),自流域上游到下游選取5個(gè)監(jiān)測斷面,分別為望濱、道義、大河泡、馬家窩棚和流域出口處(位置見圖1),結(jié)合1996—2010年降水量與模擬污染數(shù)值對流域氮素污染月際貢獻(xiàn)率進(jìn)行分析,監(jiān)測斷面月際平均徑流、TN貢獻(xiàn)率如圖3所示,由圖可知,徑流與TN月際貢獻(xiàn)率分配趨勢大體一致,自流域上游望濱到流域出口,各監(jiān)測斷面徑流與TN貢獻(xiàn)率均在月際間出現(xiàn)高峰值:望濱在5—6月;道義在6—7月;大河泡在7—11月;馬家窩棚在11—12月;流域出口在1—4月,自流域上游至下游各監(jiān)測斷面處貢獻(xiàn)率高峰值從5—6月開始呈遞推式發(fā)展,且望濱以上TN月均貢獻(xiàn)率僅占全流域的1.93%;道義至望濱占全流域的11.09%;大河泡至道義占全流域的17.91%;馬家窩棚至大河泡占全流域的25.56%;流域出口至馬家窩棚占全流域的43.51%。流域TN貢獻(xiàn)區(qū)集中在流域出口至馬家窩棚水域。
圖2 徑流與TN月擬合圖Fig.2 Runoff and TN month fitting figure
圖3 監(jiān)測斷面月際徑流、TN貢獻(xiàn)率Fig.3 Month runoff and TN contribution ratio of monitoring section
圖4 流域有機(jī)氮空間分布圖Fig.4 Basin organic nitrogen space distribution
圖5 流域NO3--N空間分布圖Fig.5 Basin NO3--N space distribution
流域內(nèi)TN平均質(zhì)量濃度達(dá)到0.70 mg/L,子流域1的TN平均質(zhì)量濃度為1.93 mg/L,根據(jù)國家環(huán)境保護(hù)總局頒布的地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)[15]判別,TN(以N計(jì))質(zhì)量濃度≤2 mg/L屬于Ⅴ類水質(zhì),流域內(nèi)子流域1屬于氮素污染關(guān)鍵區(qū),應(yīng)特別注意統(tǒng)籌規(guī)劃,與此同時(shí),3、9、10號子流域TN濃度也達(dá)到了Ⅳ類水質(zhì),也是氮素污染的次要關(guān)鍵區(qū)域,對于此類地區(qū)應(yīng)加強(qiáng)管理,防治流域氮素污染加劇,惡化流域水質(zhì)。
2)根據(jù)TN月模擬結(jié)果分析可知,TN量在年內(nèi)各月呈波動(dòng)變化,受汛期影響,TN量在5—6月開始波動(dòng),7—9月達(dá)到最大,10—12月有所回落,1—4月的TN量較低,對比流域徑流月分布情況知,TN流失具有明顯的滯后性。原因在于5—6月進(jìn)入汛期,隨徑流量增加會(huì)加速吸附在土壤顆粒上氮物質(zhì)遷移和轉(zhuǎn)化,促使水體中淋溶作用和礦化作用發(fā)生發(fā)展,進(jìn)而增加水體中TN含量。由此可見,徑流量增加是導(dǎo)致TN量增加的直接原因。
3)自流域上游到下游選取5個(gè)監(jiān)測斷面,結(jié)合1996—2010年降水量與模擬污染數(shù)值對流域氮素污染月際貢獻(xiàn)率進(jìn)行分析知,自流域上游至下游各監(jiān)測斷面處貢獻(xiàn)率高峰值從5—6月開始呈遞推式發(fā)展,且流域TN貢獻(xiàn)區(qū)集中在流域出口至馬家窩棚水域。這是由于遼寧地區(qū)自5—6月降水量開始逐漸增多,隨降水徑流產(chǎn)生,流域中TN量增加,在地表徑流、側(cè)向流、滲流作用下氮素自流域上游至下游逐漸擴(kuò)散增多,且流域出口至馬家窩棚水域附近農(nóng)耕地面積較大。由于農(nóng)耕地產(chǎn)生的泥沙是化肥和農(nóng)藥的重要攜帶者[16],這也為TN量積累提供有利條件,使該區(qū)成為TN主要貢獻(xiàn)區(qū)域。
4)應(yīng)用模型模擬數(shù)據(jù),以子流域?yàn)檠芯繂卧治隽饔虻乜臻g分布知,流域有機(jī)氮和NO-3-N整體呈現(xiàn)由東北向西南擴(kuò)散,上游及中游地區(qū)含氮量高于下游地區(qū),且隨降雨量的增加而增加,子流域1的有機(jī)氮和NO-
3-N由于受地形坡度和降水豐沛影響而含量較多,子流域4和5則因林地茂盛而含量較少。根據(jù)流域內(nèi)TN平均濃度和國家環(huán)境保護(hù)總局頒布的地表水環(huán)境質(zhì)量判別標(biāo)準(zhǔn)知,流域內(nèi)子流域1為氮素污染關(guān)鍵區(qū),主要是由于地形坡度較大,隨徑流沖刷的泥沙中含有大量氮素,沉積在1號子流域,導(dǎo)致TN平均濃度過高。與此同時(shí)3、9、10號子流域TN平均濃度小于1號子流域,也是氮素污染的次要關(guān)鍵區(qū)域,同樣應(yīng)加以規(guī)劃管理重視。
[1] 葉芝菡,劉寶元,符素華,等.土壤侵蝕過程中的養(yǎng)分富集率研究綜述[J].中國水土保持科學(xué),2009,7(1):124-130
[2] Edwin D O.Control of water pollution from agriculture[M].Rome:Food and Agriculture Organization of the United Nations Rome,1996:40
[3] 苑韶峰,呂軍,俞勁炎.氮、磷的農(nóng)業(yè)非點(diǎn)源污染防治方法[J].水土保持學(xué)報(bào),2004,18(1):122-125
[4] Paul C M B.Nutrient Emissions from Agriculture in the Netherlands,Causes and Remedies[J].Water Science and Technology,1996,33(4/5):183-189
[5] 蔣鴻昆,高海鷹,張奇.農(nóng)業(yè)面源污染最佳管理措施(BMPs)在我國的應(yīng)用[J].農(nóng)業(yè)環(huán)境與發(fā)展,2006,23(4):64-67
[6] 郭紅巖,王曉蓉,朱建國,等.太湖流域非點(diǎn)源氮污染對水質(zhì)影響的定量化研究[J].農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào),2003,22(2):150-153
[7] 李海杰.吉林省雙陽水庫匯水區(qū)農(nóng)業(yè)非點(diǎn)源污染研究[D].長春:吉林農(nóng)業(yè)大學(xué),2007:65-74
[8] 萬超,張思聰.基于GIS的潘家口水庫面源污染負(fù)荷計(jì)算[J].水力發(fā)電學(xué)報(bào),2003(2):62-68
[9] 王中根,朱新軍,夏軍,等.海河流域分布式SWAT模型的構(gòu)建[J].地理科學(xué)進(jìn)展,2008,27(4):1-6
[10]王林,陳興偉.基于3個(gè)站點(diǎn)校準(zhǔn)與驗(yàn)證的晉江流域徑流模擬[J].中國水土保持科學(xué),2007,5(6):21-26
[11]謝元博.涇河流域非點(diǎn)源污染特性分析與SWAT模型模擬[D].西安:西安理工大學(xué),2010:60-67
[12]陳媛,郭秀銳,程水源,等.基于SWAT模型的三峽庫區(qū)大流域不同土地利用情景對非點(diǎn)源污染的影響研究[J].農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào),2012,31(4):798-806
[13]羅巧,王克林,王勤學(xué).基于SWAT模型的湘江流域土地利用變化情景的徑流模擬研究[J].中國生態(tài)農(nóng)業(yè)學(xué)報(bào),2011,19(6):1431-1436
[14]孟曉云,于興修,泮雪芹.云夢湖流域土地利用變化對非點(diǎn)源氮污染負(fù)荷的影響[J].環(huán)境科學(xué),2012,33(6):1789-1794
[15]易雯.《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》中氮、磷指標(biāo)體系及運(yùn)用中有關(guān)問題的探討[J].環(huán)境保護(hù),2004(8):10-11
[16]楊巍,湯潔,李昭陽,等.基于SWAT模型的大伙房水庫匯水區(qū)徑流與泥沙模擬[J].水土保持研究,2012,19(2):77-81