流域非點源氮磷污染負荷分布模擬

張 婷，高 雅，李建柱，孫博聞，康愛卿

(1.天津大學水利工程仿真與安全國家重點實驗室，天津 300072；2.中國水利水電科學研究院流域水循環模擬與調控國家重點實驗室，北京 100038)

水資源是人類生活、生產的重要保障，中國的水資源僅為世界平均水平的1/4，水資源短缺問題比較嚴重[1]。近年來，隨著社會經濟的發展，點源污染和非點源污染導致的水環境問題也越來越凸顯，進一步加劇了我國水資源的供需矛盾。

農業非點源污染是區域水體富營養化的重要原因，針對非點源污染時空分布特征等問題，目前常采用SWAT模型進行流域非點源污染負荷計算與水環境對流域管理措施的響應研究[2]。劉博等[3]采用SWAT模型對沙河水庫流域內的非點源污染進行模擬，通過設置不同情景分析了汛期污染與全年污染負荷間的關系，得到了不同土地利用條件下的污染負荷情況。宋蘭蘭等[4]采用SWAT(soil and water assessment tool)模型對復新河氮磷營養物流失及入河負荷與用地類型的關系進行了研究，計算得出了流域內氮磷污染排放量與入河量間的關系，解決了污染物實地監測數據不足的問題。但SWAT模型本身存在缺陷，Pohlert等[5]采用能反映模型模擬值與實測值擬合程度的Nash效率系數(Ens)對SWAT模型在不同時間尺度下的精度進行了評價，發現模型針對月尺度下氮素負荷的Ens可達0.60，但日尺度下Ens則會降低至0.15以下，無法滿足模擬精度需求。Yang等[6]針對SWAT模型自帶的經驗公式和數據庫不適用于中國的問題，改進了SWAT模型中的營養物循環算法，提升了SWAT模型適用性，并應用于海南省松濤水庫，分析了不同化肥施用情景下總氮總磷的負荷情況。

潘家口水庫是引灤入津的源頭，擔負著向天津供水的重任。20世紀90年代以來，潘家口水庫流域年均降水量明顯較少，年均氣溫顯著升高[7]，入庫徑流量明顯減少，除水量短缺的問題外，流域的水質情況也逐漸惡化，嚴重影響了天津市用水安全問題。2016年，潘家口水庫甚至因水質問題而未能向天津市供水，作為流域內的重要水源地，潘家口水庫以上流域水污染問題也應予以重視。

本文擬對SWAT模型農業及土壤數據庫進行調整補充，同時通過多站點率定，提高模型在潘家口水庫流域的適用性。另外，通過非點源污染調查、監測和模擬分析，估算非點源污染負荷，并借助SWAT模型分析非點源污染時空分布規律，揭示潘家口水庫流域氮磷污染的主要驅動因子，研究不同水文條件、土地利用下流域氮磷污染特征。

1 流域概況與SWAT模型構建

1.1 流域概況

灤河發源于河北省承德市境內的巴彥古爾圖山北麓，流經河北省、內蒙古自治區及遼寧省等地區，最終于河北省樂亭縣匯入渤海[8]，流域基本概況如圖1所示，范圍為39°10′N～42°30′N，115°30′E～119°15′E[9]。潘家口水庫為引灤入津的主體工程，由一座攔河大壩和2座副壩組成，水庫控制流域面積33 700 km2，水庫總容量29.3億m3。

圖1 研究區域地理位置Fig.1 Geographic location of study area

潘家口水庫流域內水體污染分為點源污染和非點源污染2種。點源污染排放主要來源于工業直排、生活直排與污水處理廠。據統計，流域內COD、氨氮、總氮與總磷排放總量分別為8 667.5 t、771.3 t、3 652.7 t和183.3 t。

非點源污染劃分為農村生活污染源、化肥農藥污染源、禽畜養殖污染物3方面。(a)農村生活污染采用排放系數法進行估算，依據張大弟等[10]、王桂玲等[11]的研究成果確定各類污染物排放系數，然后乘以各行政區人數得到污染物產生量。(b)化肥農藥污染采用《全國水資源綜合規劃地表水水質評價及污染物排放量調查估算工作補充技術細則》計算公式進行計算，總磷、總氮產生量計算公式：總氮=(氮肥+復合肥×0.3+磷肥×0.185)，總磷=(磷肥+復合肥×0.3)×0.436，氨氮=總氮×10%，COD=總氮×30%。(c)畜禽養殖污染物產生量采用式(1)進行計算。

(1)

式中：W——畜禽養殖污染物產生量，kg；Qi——某種畜禽的飼養數量,頭(只)；Pi——該種畜禽的飼養周期，d；Ri——該種畜禽的排泄系數，kg/(頭(只)·d)；Ci——該種畜禽糞便中的污染物平均含量，%。

3種非點源污染計算中涉及的人口數量、化肥使用量及畜禽飼養量數據來自《2018年河北農村統計年鑒》，最終計算結果見表1。

表1 2017年潘家口水庫流域非點源污染排放量Table 1 Emission amount of non-point source pollution in Panjiakou Reservoir Basin in 2017 t

1.2 模型建立

SWAT模型是美國農業部開發的流域尺度分布式模型，可用于模擬單一流域或多個具有水文聯系流域的水文過程[12]。建模需要的數據包括：分辨率為30 m×30 m的數字高程模型DEM數據；三道河子、寬城、下河南等10個雨量站1998—2018年逐日降水數據；多倫、圍場、承德及青龍4個氣象站同時段逐日數據；來自韓家營、三道河子等6個水文站的實測水文資料。模型將潘家口水庫以上流域劃分為34個子流域(圖2(a))；由中國科學院資源科學與數據中心(http://www.resdc.cn/)獲得的精度1 km的土地利用數據(圖2(b))；由HWSD土壤數據庫(http://www.fao.org)獲得的土壤數據(圖2(c))。

圖2 潘家口水庫流域概況Fig.2 General situation of Panjiakou Reservoir Basin

由于潘家口水庫流域內主要涉及的行政區域為下游的承德市，且流域上游地處山區，人口數量較少，并分散在不同的行政區內，因此以承德市的非點源污染排放量作為流域的非點源污染排放總量用于模型輸入。農村生活污染與畜禽污染年內污染負荷產生量較為固定，在SWAT模型中將2種污染分別根據其各自TN、TP、NH3-N的比例，定義新的化肥類型，以連續施肥的形式進行添加。施肥持續時間設置為1 a，施肥頻次設置為每日1次；對于化肥污染，經查閱年鑒發現，承德市主要種植作物為玉米、大豆、馬鈴薯，生長季節約在每年3—10月，設置施肥時間為每年同時段，施肥頻次設置為每月1次；對于點源污染，由于收集到的資料為年排放總量，因此將其概化為每月輸出量相等，即年內12個月均勻排放至流域。另外，根據流域內10個雨量站1959—2017年共59 a的數據計算得出多年平均降雨量并應用P-Ⅲ曲線進行擬合，設置豐水年、平水年、枯水年設計頻率分別為25%、50%、75%，選取模擬年份中最接近該頻率的年份，分別為2001年、2011年與2014年，作為典型年進行非點源污染規律分析。

1.3 模型率定

選用SWAT-CUP軟件中的SUFI-2算法對模型進行敏感性分析與率定驗證，各參數取值見表2～3。采用2000—2010年和2011—2015年承德、韓家營、潘家口等6個水文站實測徑流資料分別用于模型的率定與驗證；水質采用2017—2018年大杖子(一)等3個國控斷面數據進行率定驗證。

表2 潘家口水庫流域徑流參數率定結果Table 2 Calibration results of runoff parameters in Panjiakou Reservoir Basin

各站點徑流實測值與模擬值的吻合程度較高，R2、Ens均高于0.50，徑流模擬效果較好。因為未獲取到泥沙實測數據，水質率定效果與徑流相比較差，但R2均高于0.57，說明總氮、總磷模擬精度也能滿足要求。模擬值與實測值評價結果見表4～5。

表3 潘家口水庫流域水質參數率定結果Table 3 Calibration results of water quality parameters in Panjiakou Reservoir Basin

表4 各水文站徑流模擬結果Table 4 Simulated runoff results at each hydrological station

表5 國控斷面水質模擬效果評價Table 5 Effect evaluation of water quality simulation at national control sections

2 結 果 分 析

2.1 時間變化結果分析

2.1.1 年際污染負荷特征

在流域徑流形成過程中，降雨、徑流均會受到時間周期的影響，并表現出一定的特性。因此，有必要從年際和年內2個不同的時間尺度分析降雨、徑流與污染負荷間的關系。采用SWAT模型對2011—2018年各年非點源總氮、總磷與徑流量年際變化情況進行模擬分析，結果如圖3所示。

圖3 潘家口水庫流域模擬污染負荷年際變化Fig.3 Interannual variation of simulated pollution load in Panjiakou Reservoir Basin

從圖3可以看出，TN、TP污染負荷各年之間雖存在差異，但與流域出口徑流量關系基本一致。通過計算發現，TN、TP與流域出口徑流量相關系數分別為0.65和0.96。各年泥沙輸出量與TN、TP間的相關系數分別為0.50和0.94。

年際尺度上，TN污染負荷與徑流相關性比泥沙相關性更強，通過對流域內多年TN污染負荷成分進行模擬分析可知，硝態氮、有機氮、氨氮、亞硝酸氮的成分占比分別為53.93%、41.92%、3.92%和0.21%，TN污染中可溶性氮與有機氮負荷較為平均，這是導致上述相關性產生差異的原因之一。而TP與徑流和泥沙的相關性則均為密切相關，說明TP污染物隨徑流沖刷和土壤侵蝕而流失。對流域內總磷成分進行模擬分析發現，顆粒態無機磷的占比高達73.67%。顆粒態無機磷易受沖刷影響，往往隨泥沙流失匯入河道[13]。泥沙產量與污染物總量在豐水年份會有明顯提升，因此，在豐水年控制好流域內的水土流失，采用更加有效的污染控制手段，是減小流域污染負荷的重要途徑。

以2017年的經濟、人口、污染設置為基準，分析不同水平年來水條件下的非點源污染年際變化規律，由表6可知，隨著降水量的減少，地表產水量、TN和TP負荷也都隨之減小，其變化程度與降水量變化密切相關。豐水年的TN負荷是平水年的1.68倍，是枯水年的2.12倍；豐水年TP負荷是平水年的1.52倍，是枯水年的2.21倍。而豐水年、平水年、枯水年3個典型年的氮磷污染入河形態并未有明顯差別，其中TN污染以硝態氮為主，通過淋溶作用進入水體，并在擴散轉移后形成污染[14]，對河道水質情況造成了顯著影響。TP污染則以顆粒態無機磷為主，一般隨泥沙流失進入河道[15]。

表6 不同水平年非點源污染負荷Table 6 Non-point source pollution load in different hydrological years

分析不同水平年下的污染負荷情況與產水量的相關性?？梢园l現平水年與豐水年氮磷污染負荷與產水量相關性較好，TN與產水量間相關系數分別為0.64、0.61；TP與產水量間相關系數則分別為0.69、0.71，呈顯著正相關，說明其明顯受到產水量影響；而枯水年產水量與污染負荷間的相關性則有明顯下降，說明當產水量降低到一定閾值后，污染負荷產生量不是僅由產水量決定，而是受到污染產生方式及沖刷強度等多種因素共同影響。

2.1.2 年內污染負荷特征

根據2011—2018年連續8 a的各月模擬結果分析潘家口水庫逐月入庫污染負荷情況，結果如圖4所示。

圖4 潘家口水庫流域年內降雨徑流與TN、TP負荷分布情況Fig.4 Distribution of annual rainfall,runoff,TN load and TP load in Panjiakou Reservoir Basin

潘家口水庫流域屬溫帶季風氣候，年內降水分布不均，6—8月降水可占全年降水量的70%以上。受降雨徑流因素影響，TN、TP污染負荷也主要集中于汛期，并隨流量與降水量波動呈現出“非汛期低、汛期高”的變化趨勢。6—8月TN污染負荷輸出約占全年的72%；TP污染時間分布較TN污染更為集中，6—8月輸出96%。

從表7可以看出，TN、TP與降雨徑流均有極顯著的相關關系，由此也驗證了對于非點源污染負荷，在時間尺度上降雨是最根本的驅動因子，其形成徑流攜帶污染物到河道；而徑流則是直接促成因子，可直接將河道中的污染物攜帶至流域出口[16]?？傮w來看，汛期是潘家口水庫一年中污染負荷最嚴重的時期，若想有效控制流域內的非點源污染，對汛期的污染控制尤為關鍵。

表7 年內降雨徑流與污染負荷相關系數Table 7 Correlation coefficient between annual rainfall,runoff,and pollution load

為進一步分析不同水文條件下的年內污染負荷情況，明確非點源污染變化規律，還選取了豐、平、枯典型水平年進行污染負荷變化趨勢分析。

從圖5可以看出，三期典型年TN和TP的污染負荷變化趨勢與降雨徑流變化趨勢不盡相同。TN污染負荷變化與峰現時間和徑流的符合趨勢較好，經計算發現各水平年下二者相關系數均高于0.87；TP污染負荷變化與峰現時間則和降雨的符合趨勢較好，除枯水年外，相關系數均高于0.89。分析原因，對TN來說，降水僅能反應剝蝕土壤的能量，無法反應輸出氮的能量，徑流是土壤前期含水量和降水的綜合產物，是氮素遷移的主要途徑[17]，因此相比于降雨，TN污染負荷與徑流間的相關性更為顯著。另外，從污染角度看，枯水年相對于平水年的負荷峰值基本沒有下降，而流量則下降接近50%，可以看出枯水年汛期非點源TN濃度更高，污染情況更為嚴重。

TP污染負荷峰值一般出現在流量峰值前一月，分析豐水年、平水年TP與降水的相關性系數分別為0.89、0.91，均高于當年與徑流的相關系數(0.80、0.84)，說明相較于徑流，TP污染受降水因素影響更大，更易受沖刷影響，這也是由總磷污染物組分所決定的。

2.2 空間變化結果分析

潘家口水庫流域多年平均非點源污染各種土地利用類型TN、TP輸出負荷量見表8?？傮w來看流域多年平均TN負荷量約為2 600 t；多年平均TP負荷量約為17 t。來自耕地的TN負荷與TP負荷分別占到了總量的49.35%和69.76%，產生量遠大于其他類別土地利用類型，證明了農業生產與農村生活污染是重要的非點源污染源。

除對多年平均污染負荷進行分析外，本文對典型水平年下的TN和TP污染負荷情況也進行了分析，以探尋不同水文條件下污染負荷空間分布變化情況，從圖6可以看出，污染負荷總量從豐水年到枯水年呈現下降趨勢。

圖6 典型水平年非點源污染負荷分布Fig.6 Distribution of non-point source pollution load in typical hydrological years

由表9可發現，人口數量、降水量、坡度與產水量與非點源污染負荷的相關系數均在0.25～0.45之間，說明在空間尺度上，TN、TP的污染輸出綜合受到以上多種因素的綜合影響。相比其他因素，總氮、總磷的污染負荷均與子流域泥沙產量關系密切，相關系數分別為0.80和0.95。說明泥沙產出是非點源污染負荷的重要驅動因子，總磷負荷與泥沙的相關性更強，這也與潘家口水庫流域的總磷主要以吸附態磷流失有關。

表9 非點源污染TN和TP污染負荷空間分布與各影響因素的相關系數Table 9 Correlation coefficient of total nitrogen and total phosphorus loads of non-point source pollution and influencing factors

3 結 論

a.從非點源污染排放總量角度進行分析，2017年研究流域TN非點源排放量為94 769.73 t，TP非點源排放量為16 591.29 t，其中最主要的污染源為農村生活污染和化肥農藥污染。進一步對潘家口水庫流域內污染物形態分析發現，TN主要以硝態氮與有機氮形式流失，TP主要以顆粒態無機磷形式流失。

b.分析污染物各水平年時空分布特征發現，相較于平水年，枯水年年內的TN、TP負荷峰值基本沒有下降，而流量則下降接近50%，從而可以說明枯水年汛期非點源總氮濃度更高，污染情況更為嚴重。在對潘家口水庫流域內污染進行治理時，更應關注枯水年下的污染削減情況。在空間尺度上，TN、TP負荷量從大到小排序依次為耕地、林地、草地、荒地，說明農業生產與農村生活污染是潘家口水庫流域重要的非點源污染源。

c.對污染影響因素進行分析發現，豐水年、平水年TP與降水的相關性系數分別為0.89、0.91，均高于當年與徑流的相關系數(0.80、0.84)，說明相較于徑流，TP污染受降水因素影響更大，更易受沖刷影響。TN、TP的污染負荷均與子流域泥沙產量關系密切，相關系數分別為0.80和0.95，說明泥沙是流域非點源污染負荷的重要驅動因子。