雷湘齡 葉長(zhǎng)青 朱麗蓉 薛楊
(1.海南省林業(yè)科學(xué)研究院 海南???571100;2.海南大學(xué)生態(tài)與環(huán)境學(xué)院 海南???570228;3.海南文昌森林生態(tài)系統(tǒng)國(guó)家定位觀測(cè)研究站 海南文昌 571300;4.海南大學(xué)旅游學(xué)院 海南海口 570228)
隨著各國(guó)對(duì)水污染問(wèn)題的重視,點(diǎn)源污染已逐步得到有效的控制和治理,非點(diǎn)源污染成為水污染的主要來(lái)源[1],北京密云水庫(kù)、安徽巢湖、天津于橋水庫(kù)、上海淀山湖、云南洱海等水域,非點(diǎn)源污染比率超過(guò)了點(diǎn)源污染[2]。土地利用類型的變化使土地下墊面發(fā)生改變,不適當(dāng)?shù)耐恋乩梅绞胶娃r(nóng)業(yè)管理措施在降雨-徑流驅(qū)動(dòng)因子的作用下加劇了土壤侵蝕,化肥中過(guò)量的氮、磷污染物附著于泥沙顆粒及可溶性氮磷污染物進(jìn)入水體,從而形成非點(diǎn)源污染[3-5]。
水文模型模擬是定量描述流域非點(diǎn)源污染過(guò)程的重要手段[6-7]。SWAT(Soil and Water Assessment Tool)模型被稱為是以農(nóng)業(yè)和森林為主的流域中較有前途的具有連續(xù)模擬能力的非點(diǎn)源污染模型。Anna等[8]使用SWAT模型對(duì)歐洲南部的多瑙河流域水質(zhì)進(jìn)行評(píng)估,發(fā)現(xiàn)流域內(nèi)主河道的總氮和總磷主要來(lái)源于農(nóng)作物種植施用的化肥;陳健等[9]和耿潤(rùn)哲等[10]的研究結(jié)果也表明,農(nóng)業(yè)活動(dòng)所使用的化肥、農(nóng)藥等帶來(lái)的營(yíng)養(yǎng)物流失是流域內(nèi)發(fā)生非點(diǎn)源污染的主要原因。大量的研究表明,土地利用變化會(huì)對(duì)非點(diǎn)源污染產(chǎn)生影響,Sushil等[11]在印度莫羅爾流域的研究說(shuō)明,采用等高耕種和種植植被過(guò)濾帶能夠減少泥沙和土壤中營(yíng)養(yǎng)物質(zhì)的流失;劉晶晶等[12]模擬湖北省清江流域宜都段2020年的土地利用情景,結(jié)果顯示,在耕地減少、林地增加、城鎮(zhèn)面積擴(kuò)大的情景下,總氮和總磷負(fù)荷分別減少23.77%和29.53%;秦耀民等[13]對(duì)陜西黑河流域不同土地利用情景進(jìn)行模擬,發(fā)現(xiàn)流域內(nèi)非點(diǎn)源污染負(fù)荷隨著林地面積的增加、耕地和未利用地的減少而減少;Lin等[14]在伊利湖西部流域的研究結(jié)果顯示,采用合理的農(nóng)業(yè)措施能夠降低流域內(nèi)地表水體中氮、磷等營(yíng)養(yǎng)物的負(fù)荷。
近年來(lái)非點(diǎn)源污染成為當(dāng)下研究的熱點(diǎn)之一,國(guó)內(nèi)已對(duì)我國(guó)西北部、南部、東部典型流域的非點(diǎn)源污染開展研究,但鮮見(jiàn)關(guān)于熱帶地區(qū)河流的研究報(bào)道。海南省位于我國(guó)的最南端,屬于熱帶季風(fēng)氣候,松濤水庫(kù)是海南省最大的水庫(kù)。本研究以松濤水庫(kù)上游流域 1975—2014年的氣象數(shù)據(jù)和 1988、2002、2016年土地利用圖構(gòu)建SWAT模型,分析土地利用變化(砍伐熱帶雨林種植橡膠等經(jīng)濟(jì)林)對(duì)非點(diǎn)源污染氮負(fù)荷的影響,并識(shí)別造成非點(diǎn)源污染的“源-匯”區(qū)。
松濤水庫(kù)上游流域位于海南省白沙縣內(nèi)(圖1),流域內(nèi)山地多,森林資源豐富,平均氣溫高,日照長(zhǎng)、雨量豐沛,以臺(tái)風(fēng)帶來(lái)的降雨為主;多年平均降雨量為 1 859 mm,年最大降雨量可達(dá)2 637.8 mm,存在明顯的干濕兩季,11月至次年4月為干季,降雨量占全年的15.9%,5—10月為濕季,占全年降雨量的84.1%。
圖1 研究區(qū)域圖
以流域內(nèi)福才水文站點(diǎn)的水質(zhì)數(shù)據(jù)進(jìn)行模擬研究,控制的集水面積為 492.7 km2。根據(jù) 2012年海南省統(tǒng)計(jì)年鑒,白沙縣一年施用氮肥8 688 t,耕地和水旱田面積共有268.38 km2,平均氮肥施用320 kg/hm2(折純氮)。流域內(nèi)有2個(gè)排污口,分別用于南開鄉(xiāng)衛(wèi)生院排污和生活用水排污,氨氮負(fù)荷分別為0.03和0.04 t一年(表1);而根據(jù)福才站的實(shí)測(cè)數(shù)據(jù),氨氮負(fù)荷每年約有5 t;由此可見(jiàn),在該流域內(nèi)水質(zhì)污染以非點(diǎn)源為主。
表1 研究區(qū)內(nèi)點(diǎn)源排放統(tǒng)計(jì)
1.2.1 模型簡(jiǎn)介SWAT模型中,水文循環(huán)陸地階段基于水量平衡的方程為:
式中:SWt表示土壤最終含水量(mm);SWo表示第i天的土壤初始含水量(mm);t表示時(shí)間(d);Rday表示第i天的降水量(mm);Qsurf表示地表徑流總量,單位為mm/hm2;Ea表示第i天的蒸散發(fā)量(mm);Wseep表示第i天離開土壤剖面底部的滲透水流和旁通水流的水量(mm);Qgw表示第i天回歸流的水量(mm)。
修正的通用土壤流失方程(USLE)如下:
式中:sed表示某天的產(chǎn)沙量(t);Qsurf表示地表徑流總量(mm/hm2);qpeak表示洪峰流量(m3/s);areahru表示水文響應(yīng)單元HRU的面積(hm2);KUSLE表示方程中的土壤可侵蝕因子0.013 t·m2·h/(m2·t·cm);CUSLE表示土壤覆蓋與農(nóng)業(yè)管理措施因子;PUSLE表示水土措施保持因子;LSUSLE表示地形因子;CFRG表示粗糙度因子。
1.2.2 數(shù)據(jù)來(lái)源構(gòu)建SWAT模型需要建立流域的空間數(shù)據(jù)庫(kù)和屬性數(shù)據(jù)庫(kù)以保證模型模擬的精確性,構(gòu)建SWAT模型所需數(shù)據(jù)來(lái)源見(jiàn)表2。
表2 構(gòu)建模型所需的數(shù)據(jù)來(lái)源
1.2.3 模型的構(gòu)建
1.2.3 .1 數(shù)字高程模型數(shù)據(jù)(DEM)采用SWAT模型對(duì) DEM 進(jìn)行填洼,計(jì)算河流流向,對(duì)流累積柵格進(jìn)行預(yù)處理,通過(guò)設(shè)置上游匯水面積的閾值劃分子流域,共劃分16個(gè)子流域。
1.2.3 .2 土地利用圖使用 ENVI軟件的監(jiān)督分類方法,結(jié)合谷歌地球影像圖,采用人機(jī)交互的方式對(duì)30 m×30 m遙感影像進(jìn)行解譯,將研究區(qū)的土地利用類型分為七類:漿紙林、水域、橡膠林、耕地、居民區(qū)、果園、天然林。
1.2.3 .3 子流域和水文響應(yīng)單元?jiǎng)澐諷WAT模型屬于分布式的水文模型,通過(guò)離散化來(lái)體現(xiàn)流域的空間分布[15]。在子流域上疊加土地利用圖、土壤類型圖和坡度圖,進(jìn)一步劃分水文響應(yīng)單元。
1.2.4 參數(shù)的率定與驗(yàn)證使用 SWAT-CUP的SUFI-2算法對(duì)SWAT模擬輸出的數(shù)據(jù)進(jìn)行率定和驗(yàn)證。SWAT 模型中涉及徑流與營(yíng)養(yǎng)物負(fù)荷的參數(shù)眾多,本研究通過(guò)查閱相關(guān)文獻(xiàn)和 SWAT-CUP的多次迭代確定參數(shù),利用 SUFI-2 的全局敏感性分析方法評(píng)價(jià)參數(shù)的敏感性,敏感性越高表示該參數(shù)對(duì)徑流或營(yíng)養(yǎng)物負(fù)荷的貢獻(xiàn)率越大,選用的徑流參數(shù)及其率定與驗(yàn)證見(jiàn)表3、4。
表3 SWAT模型徑流敏感性參數(shù)及其率定值
以Nash-Sutcliffe效率系數(shù)Ens和確定性系數(shù)R2為評(píng)價(jià)指標(biāo),R2越趨近于 1,表明實(shí)測(cè)值與模擬值的擬合度越高[16]。當(dāng)Ens>0.75時(shí),模擬結(jié)果良好;當(dāng)0.36≤Ens≤0.75時(shí),結(jié)果令人滿意,說(shuō)明模型適用于該流域;當(dāng)Ens<0.36,模擬效果較差,表明模型并不適用[17]。
表4 徑流參數(shù)的率定與驗(yàn)證
1.2.4 .1 徑流的率定和驗(yàn)證根據(jù) Mann-Kendall法對(duì)年徑流系數(shù)突變檢驗(yàn)的結(jié)果,選取 1977—1996年為模型率定期(1975—1976年設(shè)置為模型預(yù)熱期),1997—2014年為驗(yàn)證期。徑流參數(shù)率定與驗(yàn)證的R2和Ens均達(dá)到0.90以上,說(shuō)明SWAT模型在該流域內(nèi)對(duì)徑流具有較好的模擬效果。
1.2.4 .2 氨氮的率定和驗(yàn)證福才水質(zhì)監(jiān)測(cè)數(shù)據(jù)為2012—2014年,設(shè)置2012年為率定期,2013—2014年為驗(yàn)證期。氨氮率定與驗(yàn)證的Ens均在0.54以上,說(shuō)明SWAT模型適用于該流域營(yíng)養(yǎng)物負(fù)荷的模擬,見(jiàn)表5、6。
表5 SWAT模型氨氮敏感性參數(shù)及其率定值
1.2.4 .3 總氮的率定和驗(yàn)證根據(jù)福才水文站點(diǎn)的相關(guān)監(jiān)測(cè)數(shù)據(jù)得到 2013年的總氮數(shù)據(jù),于2014年分別對(duì)下游的南豐水文站點(diǎn)氨氮和總氮濃度數(shù)據(jù)進(jìn)行監(jiān)測(cè),計(jì)算南豐站點(diǎn)氨氮占總氮的百分比,約為25%;以此比例估算福才站點(diǎn)2014年的總氮濃度,并對(duì)估算出的總氮數(shù)據(jù)進(jìn)行率定與驗(yàn)證(表7、圖2~4),以2013年為率定期,2014年為驗(yàn)證期。
圖2 1988、2002、2016年土地利用下徑流的率定與驗(yàn)證
表6 氨氮的率定與驗(yàn)證
表7 總氮的率定與驗(yàn)證
圖3 1988、2002、2016年土地利用下氨氮的率定與驗(yàn)證
圖4 1988、2002、2016年土地利用下總氮(TN)的率定與驗(yàn)證
研究區(qū)內(nèi)土地利用類型以天然林和經(jīng)濟(jì)林為主,1988—2016年漿紙林、橡膠林和耕地的面積有所增加,橡膠林面積占比由1988年的9.74%增加至 2016年的 11.30%,漿紙林由 5.19%增加至9.49%;而天然林面積逐漸減少,由 80.15%降至64.63%(表8),全年總氨空間分布情況見(jiàn)圖5。
圖5 1988、2002、2016年土地利用下的全年總氮空間分布圖
表8 土地利用類型及其面積比
降雨是導(dǎo)致地表土壤養(yǎng)分流失的外部動(dòng)力條件[18]。雨滴落下沖刷地表,使得土壤中的氮、磷等營(yíng)養(yǎng)元素流失[19],在降雨-徑流驅(qū)動(dòng)因子的作用下進(jìn)入水體,從而產(chǎn)生了非點(diǎn)源污染。松濤水庫(kù)上游流域 5—10月為濕季,降雨量占全年的84.1%,其總氮輸出約占全年 85%;干季降雨量占全年的 15.9%,總氮輸出占比約為 15%。有研究表明,山東省小清河流域6—9月的降雨量占全年的70%,在該時(shí)間段內(nèi)非點(diǎn)源污染負(fù)荷的輸出占比為全年的一半以上[20];清水河非點(diǎn)源污染物總氮、總磷與月均降雨量、產(chǎn)水量的變化趨勢(shì)基本一致,當(dāng)降雨量達(dá)到峰值時(shí),總氮、總磷輸出量也到達(dá)最大值[21]。上述研究均與本文的研究結(jié)論相符,因此在該流域內(nèi)非點(diǎn)源污染負(fù)荷以濕季輸出為主。
降雨產(chǎn)生徑流,而徑流在遷移的過(guò)程同時(shí)還受到土地利用方式的影響,因此不同的土地利用方式下非點(diǎn)源污染物的輸出存在較大的差異[22]。采用SWAT模型分別模擬1988、2002、2016年三期土地利用下總氮負(fù)荷的變化,結(jié)果表明,1988—2016年流域全年氮負(fù)荷量表現(xiàn)為顯著的增加趨勢(shì);與 1988年相比,2002年流域內(nèi)全年總氮負(fù)荷增加了114.6%,2016年增加了171.65%,在此期間流域內(nèi)天然林面積減少,耕地、漿紙林、橡膠林的面積增加。經(jīng)濟(jì)林面積的增加意味著施肥量的增加,但作物對(duì)氮肥當(dāng)季的利用率僅為30%~35%[23],剩余的營(yíng)養(yǎng)元素附著于土壤顆粒,當(dāng)發(fā)生降雨時(shí),隨徑流遷移入水體,從而造成了流域內(nèi)非點(diǎn)源污染氮負(fù)荷量的增加。
由總氮負(fù)荷的空間分布(圖5)可知,1988年土地利用下的總氮負(fù)荷量最大的是7和16號(hào)子流域,2002年是 3、5、7、8、16號(hào)子流域的總氮輸出量最大;相比于 1988年,2002年屬于總氮負(fù)荷最大梯度等級(jí)的子流域增加了3個(gè),其中3和8號(hào)子流域的耕地面積增加,5號(hào)子流域的橡膠面積增加;16個(gè)子流域中總氮輸出增長(zhǎng)幅度最大的是1號(hào)子流域,漲幅為 223.58%,其耕地面積增加。2016年3、5、6、7、8、16號(hào)子流域的總氮輸出量最大,相比于1988年,總氮負(fù)荷最大梯度等級(jí)的子流域增加了4個(gè),增長(zhǎng)幅度最大的是8號(hào)子流域,漲幅為407.66%,其漿紙林和耕地面積大幅度增加??偟敵隽扛叩淖恿饔騼?nèi)均以耕地和經(jīng)濟(jì)林為主,1988—2016年間天然林面積減少,由人工開墾種植的耕地、經(jīng)濟(jì)林面積增加;土地下墊面的改變,施肥量的增加,造成氮負(fù)荷量的增加。
土地利用類型的單位面積氮負(fù)荷量從高到低依次為:耕地>果園>橡膠林>漿紙林>天然林。不同土地利用類型的氮負(fù)荷量與植被的覆蓋度相關(guān),土壤中氮元素的流失是地表徑流與土壤中氮素相互作用的結(jié)果,當(dāng)植被覆蓋度增加,土壤顆粒與地表徑流充分作用,使水流速度減慢,同時(shí)使得隨著徑流遷移的粗顆粒沉淀;反之,泥沙中細(xì)顆粒含量增加,所以植被覆蓋度的增加有利于泥沙中細(xì)顆粒的富集,從而導(dǎo)致泥沙中全氮含量增加[24]。在研究區(qū)內(nèi),耕地的單位面積氮負(fù)荷量最大,達(dá)到7.70 kg/hm2,在同等施肥條件下,相比于果園、橡膠林和漿紙林,耕地的植被覆蓋率最低;其次是果園,而天然林在當(dāng)?shù)厣L(zhǎng)多年,植被覆蓋率高,無(wú)人工干預(yù)不施用化肥,因此單位面積氮負(fù)荷最低,為1.76 kg/hm2,其土壤中的氮營(yíng)養(yǎng)物質(zhì)主要來(lái)源于生物固氮和微生物[25]。
土地利用類型對(duì)全年總氮輸出的貢獻(xiàn)率由大到小依次為:天然林>橡膠林>漿紙林>耕地>果園。天然林的單位面積氮負(fù)荷量最小,但其面積占比最大,因此輸出氮負(fù)荷量的貢獻(xiàn)率最大。1988—2016年,不同土地利用類型的氮負(fù)荷貢獻(xiàn)率與其面積變化相關(guān),耕地、橡膠林、漿紙林的面積呈增加的趨勢(shì),氮負(fù)荷貢獻(xiàn)率隨之增加,如1988年的耕地面積為25.56 km2,2016年為46.75 km2,其貢獻(xiàn)率從3.37%漲到了15.25%(表9);1988—2016年漿紙林面積增加了21.19 km2,氮負(fù)荷貢獻(xiàn)率增加了一倍;相反,天然林面積減少了76.43 km2,貢獻(xiàn)率從1988年的58.4%下降至2016年的34.33%。
對(duì)1988、2002、2016年研究區(qū)內(nèi)各子流域不同土地利用類型的面積與其輸出的總氮量進(jìn)行Pearson相關(guān)性分析,見(jiàn)表8。1988年天然林與總氮負(fù)荷呈極顯著相關(guān);2002年果園與總氮負(fù)荷顯著相關(guān),漿紙林、天然林與總氮負(fù)荷呈極顯著相關(guān);2016年耕地、橡膠林與總氮負(fù)荷顯著相關(guān),漿紙林、天然林與總氮負(fù)荷呈極顯著相關(guān)。天然林在流域內(nèi)面積占比最大,因此在三期土地利用下均與總氮負(fù)荷呈極顯著相關(guān),但隨著面積的減小,其相關(guān)系數(shù)逐年降低;相反,耕地、橡膠林、漿紙林的面積隨年限有所增加,氮負(fù)荷量隨之增加,其相關(guān)系數(shù)逐漸增大,甚至在2016年達(dá)到顯著或極顯著相關(guān)。
將景觀生態(tài)學(xué)的“源-匯”理論應(yīng)用于非點(diǎn)源污染中,“源”是土壤養(yǎng)分流失的源頭,“匯”則是能夠接納遷移物質(zhì),并對(duì)非點(diǎn)源污染具有滯緩作用的場(chǎng)所[26-27]。結(jié)合不同土地利用類型氮負(fù)荷的貢獻(xiàn)率(表9)和Pearson相關(guān)分析(表10),研究區(qū)內(nèi)耕地在1988和2002年與總氮負(fù)荷呈負(fù)相關(guān),但隨著面積的增加,2016年達(dá)到顯著相關(guān),其貢獻(xiàn)率也逐年增加;橡膠林、漿紙林也隨著面積的增加,相關(guān)系數(shù)和貢獻(xiàn)率增加,因此耕地、橡膠林和漿紙林對(duì)非點(diǎn)源污染氮負(fù)荷起到“源”的作用。天然林由于在流域內(nèi)面積占比最大,其貢獻(xiàn)率最高且與總氮負(fù)荷呈極顯著相關(guān),而天然林具有涵養(yǎng)水源、保持水土等生態(tài)功能,因此起到“匯”的作用,滯緩非點(diǎn)源污染的發(fā)生。
表9 不同土地利用類型的總氮負(fù)荷
表10 土地利用類型的面積與總氮負(fù)荷的Pearson相關(guān)系數(shù)
降雨是非點(diǎn)源污染產(chǎn)生的直接驅(qū)動(dòng)力[28],研究區(qū)內(nèi)存在明顯的干濕兩季,濕季降雨量占全年的 84.1%,其總氮輸出量占全年的 85%,總氮負(fù)荷與月降雨量的變化趨勢(shì)基本一致,隨著降雨量的增加,流域內(nèi)總氮輸出升高。
在降雨條件下,非點(diǎn)源污染的負(fù)荷強(qiáng)度還受到土地利用方式的影響。不同土地利用類型的單位面積氮負(fù)荷量從高到低依次為:耕地>果園>橡膠林>漿紙林>天然林,其與植被覆蓋度和化肥施用量相關(guān),因此耕地的單位面積氮負(fù)荷量最大,相反,天然林最小。同理,流域內(nèi)的氮負(fù)荷變化和空間分布與土地利用的分布相關(guān),1988—2016年天然林面積減少,耕地、橡膠林、漿紙林和果園種植面積增加,在此期間總氮負(fù)荷呈增加的趨勢(shì);相比于 1988年,在相同的劃分條件下 2002年總氮負(fù)荷最大梯度等級(jí)的子流域增加了3個(gè),2016年增加了4個(gè),總氮增加的子流域均是由于耕地或其他經(jīng)濟(jì)林的種植面積增加了。不同土地利用總氮輸出的貢獻(xiàn)率從高到低依次為:天然林>橡膠林>漿紙林>耕地>果園,貢獻(xiàn)率的大小與面積有關(guān);Pearson相關(guān)性分析與貢獻(xiàn)率分析的結(jié)果基本一致,天然林在流域內(nèi)所占面積最大,氮負(fù)荷貢獻(xiàn)率最高,與總氮負(fù)荷的關(guān)系均表現(xiàn)為極顯著相關(guān),但隨著面積的減小,其相關(guān)系數(shù)和貢獻(xiàn)率逐漸降低;耕地、橡膠林、漿紙林隨著面積的增加,貢獻(xiàn)率升高,與總氮負(fù)荷的相關(guān)系數(shù)逐漸增大,到了2016年達(dá)到顯著或極顯著相關(guān)。
當(dāng)前非點(diǎn)源污染已成為水體污染的一大主要污染來(lái)源,而由農(nóng)業(yè)造成的非點(diǎn)源貢獻(xiàn)率已超過(guò)50%。松濤水庫(kù)是海南省重點(diǎn)飲用水源保護(hù)區(qū),其上游流域山多地廣,森林資源豐富,以農(nóng)業(yè)為主導(dǎo)產(chǎn)業(yè),因此如何有效地防控非點(diǎn)源的發(fā)生已成為不容小覷的問(wèn)題。本研究成果可為當(dāng)?shù)胤乐无r(nóng)業(yè)非點(diǎn)源污染提供參考。
研究使用SWAT模型對(duì)松濤水庫(kù)上游流域進(jìn)行了非點(diǎn)源污染負(fù)荷的模擬,分析了土地利用變化對(duì)氮負(fù)荷的影響,結(jié)果如下:
(1)對(duì)以1988、2002和2016年土地利用數(shù)據(jù)為基礎(chǔ)構(gòu)建的SWAT模型進(jìn)行率定與驗(yàn)證,對(duì)月徑流率定和驗(yàn)證的納什效率系數(shù)Ens均為0.90以上,氨氮和總氮的Ens分別為0.54和0.46以上,說(shuō)明SWAT模型適用于該流域的非點(diǎn)源污染模擬。
(2)降雨是造成非點(diǎn)源污染的外部動(dòng)力條件,濕季總氮負(fù)荷全年輸出約占85%,干季總氮負(fù)荷全年約占15%。1988—2016年流域內(nèi)氮負(fù)荷量呈增加的趨勢(shì),氮負(fù)荷的空間分布與土地利用類型相關(guān),總氮輸出量高的子流域內(nèi)均以耕地和經(jīng)濟(jì)林為主。
(3)單位面積氮負(fù)荷量從高到低依次為:耕地>果園>橡膠林>漿紙林>天然林;土地利用類型對(duì)全年總氮輸出的貢獻(xiàn)率由大到小依次為:天然林>橡膠林>漿紙林>耕地>果園;天然林與氮負(fù)荷的關(guān)系為極顯著相關(guān),隨著面積的減小,其Pearson相關(guān)系數(shù)逐漸減?。浑S著面積的增加,耕地、橡膠林、漿紙林與總氮負(fù)荷的Pearson相關(guān)系數(shù)逐漸增大,在2016年達(dá)到顯著或極顯著相關(guān)。