甬江流域土地利用方式對面源磷污染的影響：基于SWAT模型研究

張招招，程軍蕊，畢軍鵬，徐宇婕，李 昀，王 侃*

（1.寧波大學地理與空間信息技術(shù)系，寧波 315211；2.寧波大學生態(tài)環(huán)境研究所，寧波 315211）

土地利用方式?jīng)Q定了下墊面性質(zhì)和地表景觀結(jié)構(gòu)，影響著流域水循環(huán)與物質(zhì)遷移等生態(tài)過程[1]。農(nóng)作物化肥施用、畜禽養(yǎng)殖和農(nóng)村居民生活排放大量的磷素負載于不同的土地利用類型[2]，通過徑流、土壤侵蝕和農(nóng)田排水等不同的地表流失機制以面源污染的方式進入地表水環(huán)境[3]，造成受納水體的水質(zhì)質(zhì)量下降[4]。而磷素不僅是生物體生長的必要元素之一，也是水體富營養(yǎng)化及藻類生長的關鍵因子[5]，因此如何定量分析土地利用方式對面源污染磷負荷的影響，已成為地表水環(huán)境保護和水體富營養(yǎng)化研究的難點和熱點之一。

目前，土地利用方式對面源磷污染的影響研究主要有小區(qū)試驗法[6]、輸出系數(shù)法[7]和機理模型法等方法。小區(qū)實驗法與輸出系數(shù)法難以開展長時間序列、大尺度的研究，促使機理模型迅速發(fā)展。Borah和Be?ra[8-9]從數(shù)學基礎和模型應用兩大方面綜述了AG?NPS、HSPF、SWAT等11種流域面源污染模型，認為SWAT最適用于以農(nóng)業(yè)占主導地位的流域連續(xù)模擬。國外對SWAT模型的應用研究主要有水質(zhì)評估與預測、面源污染負荷估算及情景分析、農(nóng)業(yè)管理措施對水文水質(zhì)的影響等諸多領域[10]，而在國內(nèi)SAWT模型研究主要為各中小尺度流域徑流、泥沙和水質(zhì)的模擬[11]，農(nóng)業(yè)措施管理[12]、氣候變化[13]、土地利用情景等對水文水質(zhì)[14]的影響等研究。流域水質(zhì)模擬和土地利用的集成應用研究相對較少[15]，相關研究更多關注土地利用歷史變化和未來情景模擬下的水文水質(zhì)變化[16-17]，且從子流域尺度去探究污染負荷問題[18]，忽略了子流域內(nèi)部土地利用類型、土地利用結(jié)構(gòu)、下墊面性質(zhì)等差異的影響，同時研究更注重污染物時空特征，缺少不同土地利用類型單位面積污染負荷差異比較。甬江是寧波地區(qū)的母親河，控制水體污染持續(xù)惡化既緊迫又繁復。本文利用SWAT水文水質(zhì)模型，從攜帶多種下墊面信息的水文響應單元（HRU）尺度出發(fā)，探究甬江流域不同土地利用方式之間以及相同土地利用方式內(nèi)部面源磷污染負荷的差異及其影響，以期為甬江流域的水環(huán)境污染治理與規(guī)劃管理提供基礎理論支撐，并為其他流域的水環(huán)境問題治理提供科學參考。

1 材料與方法

1.1 研究區(qū)概況

甬江流域位于浙江省東北部的東海之濱，是浙江省八大水系之一，由甬江干流、奉化江和姚江共同構(gòu)成，流至寧波鎮(zhèn)海區(qū)入東海。甬江流域匯水區(qū)面積約為3 787.89 km2，地理范圍見圖1。流域匯水區(qū)域包含16個飲用水水源保護區(qū)，涉及56個水環(huán)境功能區(qū)。地貌以山丘和平原兩種類型為主，南部為天臺山脈的余脈，西南部為天臺山支脈四明山脈。流域多年平均降水量1 592.3 mm，主要集中于3—6月的春雨連梅雨和8—9月的臺風雨和秋雨。流域內(nèi)共有簡育低活性強酸土（ACh）、聚鐵網(wǎng)紋低活性強酸土（ACp）、腐殖質(zhì)低活性強酸土（ACu）、簡育高活性強酸土（ALh）、聚鐵網(wǎng)紋高活性強酸土（ALp）、人為土（ATc）、不飽和雛形土（CMd）、鐵鋁性雛形土（CMo）、腐殖質(zhì)雛形土（CMu）、石灰性沖積土（FLc）、飽和沖積土（FLe）、飽和潛育土（GLe）、不飽和疏松巖性土（RGd）、鈉質(zhì)鹽土（SCh）14種土壤類型。

圖1 甬江流域地理位置圖Figure 1 Geographic position of Yong River Basin

1.2 數(shù)據(jù)來源

構(gòu)建研究區(qū)SWAT模型，需要收集研究區(qū)域的土地利用、土壤、氣象等大量基礎數(shù)據(jù)，本文主要數(shù)據(jù)及其來源見表1。

1.3 主要數(shù)據(jù)的選擇與處理

土地利用數(shù)據(jù)：在第二次全國土地調(diào)查數(shù)據(jù)中甬江流域包括林地、水田、城市、村莊等32種二級用地類型，本文參照《土地利用現(xiàn)狀分類》（GB/T 21010—2017）國家標準重分類為7種一級用地類型：林地、耕地、園地、建筑用地、水域、草地、未利用地，分別占總流域面積的41.16%、27.07%、4.15%、21.18%、5.97%、0.35%、0.12%。

土壤數(shù)據(jù)：我國廣泛應用的土壤數(shù)據(jù)是以發(fā)生學分類劃分的土壤類型，比如全國第二次土壤普查數(shù)據(jù)，而模型接受的土壤數(shù)據(jù)是以土壤診斷層和診斷特性劃分的土壤類型，所以這兩類土壤數(shù)據(jù)的土壤質(zhì)地粒徑分級標準不同，因此本文選擇在國際上通用且與模型土壤質(zhì)地粒徑分級標準一致的來自于南京土壤所的中國土壤系統(tǒng)分類數(shù)據(jù)。土壤屬性數(shù)據(jù)庫中參數(shù)通過土壤屬性數(shù)據(jù)庫直接獲取與轉(zhuǎn)換、默認值以及土壤水特性計算軟件SPAW等建立。

氣象數(shù)據(jù)：目前研究者們主要利用實測站點和衛(wèi)星反演數(shù)據(jù)作為氣象模塊的輸入，以衛(wèi)星反演的再分析數(shù)據(jù)在數(shù)據(jù)的選擇、降尺度的方法的處理過程中不可避免地降低模擬數(shù)據(jù)精度。因此本文選擇傳統(tǒng)氣象監(jiān)測站的數(shù)據(jù)作為輸入，采用1990—2010年鄞州和慈溪氣象站逐日氣象數(shù)據(jù)構(gòu)建了研究區(qū)天氣發(fā)生器，并統(tǒng)計處理模擬期2008—2014年的氣象要素數(shù)據(jù)。

污染源數(shù)據(jù)：污染源數(shù)據(jù)主要包括點源數(shù)據(jù)和面源數(shù)據(jù)。點源主要有來自于污水處理廠及工業(yè)企業(yè)的廢水排放直接統(tǒng)計，面源主要有農(nóng)村生活、畜禽養(yǎng)殖、農(nóng)業(yè)種植等污染源，此類型數(shù)據(jù)中以行政單元進行統(tǒng)計的數(shù)據(jù)，根據(jù)平均密度法進行子流域數(shù)據(jù)計算，即將行政區(qū)統(tǒng)計數(shù)據(jù)理想化為均勻分布。①農(nóng)村生活污染：對于納管的村莊污染物負荷按0計；而未納管處理的污染物排放量按排放系數(shù)法估算。②畜禽養(yǎng)殖：豬、牛、羊、家禽等畜禽養(yǎng)殖數(shù)量統(tǒng)計，根據(jù)畜禽養(yǎng)殖主要污染物排放系數(shù)進行計算。③農(nóng)業(yè)種植：寧波市地區(qū)歷史上以雙季稻種植為主，近年來種植模式逐步轉(zhuǎn)變?yōu)殡p季稻、單季稻與小麥輪作、單季稻與油料作物輪作等更加多元的種植方式[19]。復雜的種植模式和區(qū)域之間種植差異給建模帶來極大困難，本文將流域內(nèi)的種植模式概化為較為廣泛的單季稻與小麥輪作，分別在種植和生長期按比例施肥。

1.4 SWAT模型的構(gòu)建與模擬

模型設置統(tǒng)一的地球橢球體為大地2000，投影為高斯克呂格投影，柵格精度為30 m。建模過程按如下4部分完成：①子流域劃分：該模型流域劃分的原理是基于D8算法確定水系與流域信息，而甬江的平原河網(wǎng)區(qū)地勢平緩，此算法難以正確劃分，因此將部分縱橫交錯的真實河網(wǎng)水系修改為樹枝狀并用burn-in來解決。最終將子流域面積閾值設為5 km2，劃分了65個子流域。②水文響應單元（HRU）劃分：HRU考慮了子流域內(nèi)土地利用類型、土壤類型、地形的多樣性及其相互組合形成復雜的下墊面狀況，從而極大增加了模型模擬精度，本文最后劃分了既能保證模擬精度又能有效運行的676個HRU。③模型參數(shù)寫入與編輯：將按模型需求處理好的氣象站點索引表和各站各氣象要素實際監(jiān)測表導入，并寫入12類文件初步確定模型參數(shù)。通過編輯界面，輸入各子流域點源形式的排污量，編輯各子流域農(nóng)田管理措施時間，作物類型、施肥量、施肥種類等。④模型運行：模型初期所有的參數(shù)都是從零增加至模擬值，勢必會對模擬結(jié)果的精準性造成影響，因此本文設置了2年預熱期，模擬時間為2008年1月1日至2014年12月31日。

表1 建模數(shù)據(jù)及其來源Table 1 The required data and source for modeling

2 結(jié)果與分析

2.1 參數(shù)率定結(jié)果及模型驗證

本文通過SWAT-CUP 2012進行敏感性分析，從眾多參數(shù)中篩選出敏感程度較大的20個參數(shù)，借助SUFI-2參數(shù)估計最優(yōu)化方法對姚江大閘站2010—2012年徑流和總磷進行逐月率定，2013—2014年進行逐月驗證，率定期與驗證期的徑流和總磷模擬值與實測值見圖2、圖3。選取決定系數(shù)R2和納什效率系數(shù)Ens進行模擬效果適宜性評價，徑流在率定期與驗證期R2分別達0.83、0.80，Ens分別達到0.80、0.82；總磷在率定期和驗證期R2分別達到0.85、0.71，Ens分別達到0.78、0.67。一般認為，Ens大于 0.5，R2大于0.6時，模擬效果令人滿意，所以本研究構(gòu)建的甬江流域SWAT水文水質(zhì)模型可達到模擬要求，將SWAT-CUP中的參數(shù)移植到SWAT模型后進行新的模擬，其結(jié)果為最佳模擬結(jié)果。

2.2 不同土地利用方式下面源污染輸出差異分析

2.2.1 不同土地利用方式下的產(chǎn)流、產(chǎn)沙、總磷輸出

如表2所示，各土地利用方式的產(chǎn)流深度無論在何種坡度下，建設用地的不透水面導致的產(chǎn)流深度最大且至少高于其他土地利用方式100 mm以上，其次為人為干預地表程度最強的耕地較高，再其次為人為干預強度較低的園地，產(chǎn)流最小的為植被覆蓋度最高的林地。但由于產(chǎn)流強度不僅與覆蓋方式和坡度相關性強，還受到土壤類型及其屬性、降雨要素等的強烈干擾，造成各土地利用方式的產(chǎn)流能力變化趨勢各不相同。

圖2 徑流月校準與驗證結(jié)果Figure 2 Calibration and validation result for runoff

圖3 總磷月校準與驗證結(jié)果Figure 3 Calibration and validation result for total phosphorus

各土地利用方式的產(chǎn)沙負荷隨著坡度的增加而增大，且在坡度級別為3～4級時產(chǎn)沙速度變化最快，耕地的產(chǎn)沙單位負荷遠高于其他土地利用方式，分別為建設用地、林地、園地的42.10、8.96、7.27倍，園地略高于林地，而建設用地的產(chǎn)沙能力遠低于其他土地利用方式。

建設用地總磷負荷最低范圍為0.32～0.38 kg·hm-2，耕地負荷最高范圍高達2.15～11.5 kg·hm-2，而園地和林地輸出較一致范圍均在0.34～1.73 kg·hm-2，各土地利用方式的總磷輸出強度均隨坡度的增加而增強，在坡度級別為3～4級時總磷輸出強度變化最大。耕地的磷單位負荷最大，單位面積磷流失約為林地的28倍，由于居民的日常生活產(chǎn)生的磷素較少而致建設用地磷流失強度最低；林地雖人為活動干涉最少，具有最佳保持水土、削減洪峰以及減少泥沙輸送量等重要生態(tài)水文功能，但輸出強度卻略低于園地。

2.2.2 不同土地利用方式下的產(chǎn)流、產(chǎn)沙與面源磷的關系分析

為探究產(chǎn)流和產(chǎn)沙對磷輸出的影響關系，本文以土地利用方式、土壤類型和坡度均不同的676個水文響應單元（HRU）作為樣區(qū)，篩選出建設用地在坡度等級為2級而土壤類型為ACp的17個HRU、林地坡度等級為3級而土壤類型為RGd的11個HRU、園地坡度等級為3級而土壤類型為ACu的6個HRU、耕地坡度等級為2級而土壤類型為ACp的14個HRU，以每個HRU各年的產(chǎn)流深度、產(chǎn)沙負荷與總磷的輸出建立一元線性回歸模型，如圖4所示。林地、園地、耕地的產(chǎn)流深度與磷單位負荷的相關性較好（R2=0.63～0.68，P＜0.001），表明產(chǎn)流深度與磷單位負荷這二者之間存在明顯的線性關系，而建設用地的線性擬合精度卻不佳（R2=0.21，P＜0.001）主要是由于地表被人工重新塑造產(chǎn)生獨特的產(chǎn)流機制所造成的；建設用地、林地、園地的產(chǎn)沙單位負荷與磷單位負荷顯著相關（R2=0.83～0.88，P＜0.001），表明產(chǎn)沙單位負荷與磷單位負荷之間也存在明顯的線性關系，而耕地的擬合度較差（R2=0.31，P＜0.001），這是由于本研究的HRU樣本位于不同行政區(qū)內(nèi)，而不同的行政區(qū)內(nèi)施肥管理水平不同，因此導致耕地出現(xiàn)擬合度較低的結(jié)果。總的來說，產(chǎn)沙和面源磷之間的回歸模型的預測能力均高于產(chǎn)流和磷的預測能力，表明了流域內(nèi)面源磷流失的主要載體為泥沙。

2.3 同類土地利用方式下的面源污染負荷的空間負荷差異分析

相同土地利用方式下磷負荷空間差異性十分明顯，圖5為各土地利用方式下各HRU的年均磷流失負荷圖?？臻g分布差異一方面來源于化肥施用水平、作物管理方式等人為干預程度和方式的不同，另一方面來源于坡度與土壤類型等地表性質(zhì)的不同。耕地磷流失一方面受制于坡度、土壤等自然條件，另一方面受制于不同區(qū)域作物管理、施肥水平不同等人類活動強度的影響，其高強度流失區(qū)主要分布于坡度較大且化肥施用強度最高的奉化區(qū)；建設用地磷流失規(guī)律與土壤類型無關，對坡度的敏感性較低，與人類的活動相關性最強，因此磷的高強度流失區(qū)主要集中分布于人口高密度區(qū)；而林地和園地人為干擾力最弱，磷流失空間輸出主要由地表性質(zhì)不均引起，高強度流失區(qū)集中于坡度較大的區(qū)域。

建設用地和耕地在流域尺度內(nèi)各區(qū)域人為活動的影響極其復雜且具有很大的隨機性，其基礎數(shù)據(jù)難以全面且精準統(tǒng)計，因此只從定性角度去探討了面源污染輸出強度的主導因子。而林地和園地磷流失空間差異主要由地表性質(zhì)不均引起，因此以廣泛分布于流域山地、丘陵區(qū)的林地為研究對象，定量探討地表性質(zhì)中最為敏感的地表坡度和土壤類型兩因子對磷流失的影響。圖6為各土壤類型在不同坡度級別下面源磷負荷圖。各土壤類型均隨坡度的增加磷流失強度不斷遞增，坡度級別為6級時磷流失平均負荷達到最高值3.48 kg·hm-2，是坡度級別為2級時的6.90倍；其次坡度級別越高增速越快，當坡度在2～3級之間時，磷流失增長率僅為33.68%；而坡度級別為5～6級時，增長率已達85.38%。其次，不同的土壤類型在相同坡度級別下磷流失負荷也各不相同，且坡度越大因土壤類型的不同帶來的磷輸出強度差異越大，比如土壤類型RGd在坡度2級下磷流失是ACu的1.15倍，而在坡度級別為6級下RGd是ACu的1.42倍。

表2 不同坡度級別下不同土地利用方式的年均產(chǎn)流深度、產(chǎn)沙單位負荷和總磷單位負荷Table 2 Annual runoff depth，sediment loading and total phosphorus loading under different land uses for different levels of slope

圖4 各土地利用類型年均產(chǎn)流深度與總磷單位負荷、產(chǎn)沙單位負荷與磷單位負荷的關系圖Figure 4 Relationship between runoff depth and sediment loading and total phosphorus loading under different land uses

圖5 各土地利用類型面源總磷空間負荷強度（kg·hm-2）Figure 5 Land use types of non-point source TP pollution spatial loading intensity（kg·hm-2）

圖6 各土壤類型在不同坡度級別下面源磷負荷Figure 6 Non-point source pollution TP load of various soil types in different slope grades

3 討論

前人的研究表明：土壤亞表層具有較強的固磷能力，致使土壤滲流中的磷濃度較低[20]，但土壤表層中粒徑小的礦物質(zhì)和土壤有機質(zhì)可以吸附大量的磷，這些磷素隨著降水過程通過徑流溶解與泥沙搬運進入水體，隨之入河造成水體污染和肥料資源的浪費[21]。而各土地利用方式磷流失量除土壤富磷水平之外，還取決于他們的產(chǎn)流、產(chǎn)沙過程，這給土地利用對面源污染的影響的定量研究帶來困難。與傳統(tǒng)小區(qū)試驗法和輸出系數(shù)法相比，SWAT機理模型在模擬尺度、模擬效率、模擬精度方面具有較大模擬優(yōu)勢。SWAT模型以流域為研究范圍，以日為時間步長，與人為控制觀測小區(qū)條件來避免干擾因素的傳統(tǒng)小區(qū)試驗法相比，在保證模擬精度的前提下極大地拓展了研究的空間尺度與時間尺度；同時SWAT模型在大量資料觀測的基礎上通過計算機模擬，與逐月進行野外水文觀測的小區(qū)試驗法相比，模型具有極高的模擬效率；此外，SWAT模型以水文響應單元為基本模擬單元，充分考慮了研究區(qū)的下墊面性質(zhì)，并輸入研究期氣象、水文、水質(zhì)等實測資料，與研究者以經(jīng)驗確定輸出系數(shù)參數(shù)的系數(shù)輸出法相比，極大地提高了模型模擬精度。

由于傳統(tǒng)小區(qū)試驗法的模擬精度較高，因此本文利用小區(qū)試驗法所取得的成果來佐證本研究中通過SWAT機理模型所取得的結(jié)果。但同時研究區(qū)的位置、研究的時期，土地利用類型劃分方式等的不同，均會給研究結(jié)果帶來略微差異。林地和園地的人為干擾力較弱且植被覆蓋度高，其殘體覆于土壤表層有助于降水的下滲、能夠減緩降水對地表的沖擊力，因此地表產(chǎn)流速度，產(chǎn)流和產(chǎn)沙均較小，磷流失也弱。園地磷素在各坡度級別下年流失范圍為0.68～3.78 kg·hm-2之間，均值為 1.81 kg·hm-2。而 Liu等[22]通過對丹江口流域進行小區(qū)觀測取得了與本研究較一致的結(jié)果，得到橘子園的磷素損失范圍為1.3～1.7 kg·hm-2，并利用指數(shù)模型很好地描述了泥沙和磷流失的關系；而林地磷年流失范圍在0.55～2.80 kg·hm-2之間，均值為 1.42 kg·hm-2，緩坡地為 0.55 kg·hm-2。高超等[23]研究了太湖流域的竹園、板栗的年均磷流失量在平坡地上分別為0.28、0.66 kg·hm-2，竹園低于本研究的結(jié)果而板栗園略高于本研究結(jié)果；耕地由于作物的種植，對土壤頻繁的翻耕施肥等人為擾動，土壤表層有機質(zhì)含量低，其次耕地土壤容易結(jié)塊，不利于降水下滲而導致水土的保持性差，其產(chǎn)水產(chǎn)沙水平高，磷流失遠高于其他方式用地，在各坡度級別下磷年流失范圍在2.63～18.29 kg·hm-2之間，平坡耕地為2.63 kg·hm-2、而緩坡為 4.60 kg·hm-2，分別為林地、園地的 8.96、7.27倍。沈連峰等[24]表明淮河流域花卉的年均磷流失量為3.98 kg·hm-2，且產(chǎn)流、產(chǎn)沙與總磷均具有顯著的線性關系，其結(jié)果高于本研究結(jié)果主要是由于淮河流域為旱地耕作區(qū)而本研究區(qū)耕地以水田為主，而水田的磷流失較弱[25]。建設用地分布于緩坡及以下平坦區(qū)域，其地表被人工重新塑造，降水下滲極弱容易造成洪峰等完全不同于其他地類的性質(zhì)，其產(chǎn)流水平雖最高，但產(chǎn)沙水平最低，因此磷流失規(guī)律與土壤類型無關，對坡度的敏感性也較低，與人類的活動相關性最強，由于居民生活釋放較少而磷的年流失最低僅為0.32～0.38 kg·hm-2，均值為0.34 kg·hm-2。

此外，同類型土地利用方式的磷流失負荷的空間負荷差異也不容忽視。比如耕地坡度級別每增加一級磷流失平均增加1.67倍，而林地不管在什么地表狀況下相比于其他土地利用方式均具有最佳的保持水土、涵養(yǎng)水源等重要生態(tài)水文功能，但在坡度為6級土壤類型為ACh的易侵蝕區(qū)其磷流失強度達4.01 kg·hm-2，接近耕地的緩坡流失強度。因此在土地資源緊缺的社會大發(fā)展背景下，不僅要合理規(guī)劃管理不同土地利用方式之間的轉(zhuǎn)換，也要注重相同土地利用方式下不同區(qū)域轉(zhuǎn)換的適宜度。

4 結(jié)論

（1）本文以決定系數(shù)與納什效率系數(shù)驗證SWAT水文水質(zhì)模型，其徑流、泥沙、面源磷在率定期與驗證期的R2均大于0.78，Ens均大于0.67，表明模型模擬精度具有良好的可靠性。其次將模型模擬與傳統(tǒng)小區(qū)試驗法結(jié)果進行比較分析，認為模型在開展大中流域尺度，長時間序列的土地利用方式對面源磷污染影響的研究中具有明顯優(yōu)勢。

（2）建設用地的不透水面導致的年均產(chǎn)流深度最大達729.52 mm，高于其他土地利用方式100 mm以上；耕地的產(chǎn)沙單位負荷遠高于其他土地利用方式達44.68 t·hm-2，分別為建設用地、林地、園地的 42.10、8.96、7.27倍；建設用地總磷負荷最低范圍為0.32～0.38 kg·hm-2，耕地最高范圍高達2.15～11.5 kg·hm-2，園地和林地輸出較低范圍均在0.34～1.73 kg·hm-2。

（3）同類土地利用方式下磷輸出空間差異性十分顯著，耕地磷流失主要分布于施肥水平高的奉化區(qū)；建設用地高強度流失區(qū)主要集中分布于人口高密度區(qū)的建成區(qū)中心；林地和園地高強度流失區(qū)主要集中在坡度較大的土壤易侵蝕區(qū)。