張 婷，高 雅，李建柱，孫博聞，康愛(ài)卿

(1.天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2.中國(guó)水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100038)

水資源是人類生活、生產(chǎn)的重要保障，中國(guó)的水資源僅為世界平均水平的1/4，水資源短缺問(wèn)題比較嚴(yán)重[1]。近年來(lái)，隨著社會(huì)經(jīng)濟(jì)的發(fā)展，點(diǎn)源污染和非點(diǎn)源污染導(dǎo)致的水環(huán)境問(wèn)題也越來(lái)越凸顯，進(jìn)一步加劇了我國(guó)水資源的供需矛盾。

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

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

本文擬對(duì)SWAT模型農(nóng)業(yè)及土壤數(shù)據(jù)庫(kù)進(jìn)行調(diào)整補(bǔ)充，同時(shí)通過(guò)多站點(diǎn)率定，提高模型在潘家口水庫(kù)流域的適用性。另外，通過(guò)非點(diǎn)源污染調(diào)查、監(jiān)測(cè)和模擬分析，估算非點(diǎn)源污染負(fù)荷，并借助SWAT模型分析非點(diǎn)源污染時(shí)空分布規(guī)律，揭示潘家口水庫(kù)流域氮磷污染的主要驅(qū)動(dòng)因子，研究不同水文條件、土地利用下流域氮磷污染特征。

1 流域概況與SWAT模型構(gòu)建

1.1 流域概況

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

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

潘家口水庫(kù)流域內(nèi)水體污染分為點(diǎn)源污染和非點(diǎn)源污染2種。點(diǎn)源污染排放主要來(lái)源于工業(yè)直排、生活直排與污水處理廠。據(jù)統(tǒng)計(jì)，流域內(nèi)COD、氨氮、總氮與總磷排放總量分別為8 667.5 t、771.3 t、3 652.7 t和183.3 t。

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

(1)

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

3種非點(diǎn)源污染計(jì)算中涉及的人口數(shù)量、化肥使用量及畜禽飼養(yǎng)量數(shù)據(jù)來(lái)自《2018年河北農(nóng)村統(tǒng)計(jì)年鑒》，最終計(jì)算結(jié)果見(jiàn)表1。

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

1.2 模型建立

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

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

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

1.3 模型率定

選用SWAT-CUP軟件中的SUFI-2算法對(duì)模型進(jìn)行敏感性分析與率定驗(yàn)證，各參數(shù)取值見(jiàn)表2～3。采用2000—2010年和2011—2015年承德、韓家營(yíng)、潘家口等6個(gè)水文站實(shí)測(cè)徑流資料分別用于模型的率定與驗(yàn)證；水質(zhì)采用2017—2018年大杖子(一)等3個(gè)國(guó)控?cái)嗝鏀?shù)據(jù)進(jìn)行率定驗(yàn)證。

表2 潘家口水庫(kù)流域徑流參數(shù)率定結(jié)果Table 2 Calibration results of runoff parameters in Panjiakou Reservoir Basin

各站點(diǎn)徑流實(shí)測(cè)值與模擬值的吻合程度較高，R2、Ens均高于0.50，徑流模擬效果較好。因?yàn)槲传@取到泥沙實(shí)測(cè)數(shù)據(jù)，水質(zhì)率定效果與徑流相比較差，但R2均高于0.57，說(shuō)明總氮、總磷模擬精度也能滿足要求。模擬值與實(shí)測(cè)值評(píng)價(jià)結(jié)果見(jiàn)表4～5。

表3 潘家口水庫(kù)流域水質(zhì)參數(shù)率定結(jié)果Table 3 Calibration results of water quality parameters in Panjiakou Reservoir Basin

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

表5 國(guó)控?cái)嗝嫠|(zhì)模擬效果評(píng)價(jià)Table 5 Effect evaluation of water quality simulation at national control sections

2 結(jié) 果 分 析

2.1 時(shí)間變化結(jié)果分析

2.1.1 年際污染負(fù)荷特征

在流域徑流形成過(guò)程中，降雨、徑流均會(huì)受到時(shí)間周期的影響，并表現(xiàn)出一定的特性。因此，有必要從年際和年內(nèi)2個(gè)不同的時(shí)間尺度分析降雨、徑流與污染負(fù)荷間的關(guān)系。采用SWAT模型對(duì)2011—2018年各年非點(diǎn)源總氮、總磷與徑流量年際變化情況進(jìn)行模擬分析，結(jié)果如圖3所示。

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

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

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

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

表6 不同水平年非點(diǎn)源污染負(fù)荷Table 6 Non-point source pollution load in different hydrological years

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

2.1.2 年內(nèi)污染負(fù)荷特征

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

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

潘家口水庫(kù)流域?qū)贉貛Ъ撅L(fēng)氣候，年內(nèi)降水分布不均，6—8月降水可占全年降水量的70%以上。受降雨徑流因素影響，TN、TP污染負(fù)荷也主要集中于汛期，并隨流量與降水量波動(dòng)呈現(xiàn)出“非汛期低、汛期高”的變化趨勢(shì)。6—8月TN污染負(fù)荷輸出約占全年的72%；TP污染時(shí)間分布較TN污染更為集中，6—8月輸出96%。

從表7可以看出，TN、TP與降雨徑流均有極顯著的相關(guān)關(guān)系，由此也驗(yàn)證了對(duì)于非點(diǎn)源污染負(fù)荷，在時(shí)間尺度上降雨是最根本的驅(qū)動(dòng)因子，其形成徑流攜帶污染物到河道；而徑流則是直接促成因子，可直接將河道中的污染物攜帶至流域出口[16]。總體來(lái)看，汛期是潘家口水庫(kù)一年中污染負(fù)荷最嚴(yán)重的時(shí)期，若想有效控制流域內(nèi)的非點(diǎn)源污染，對(duì)汛期的污染控制尤為關(guān)鍵。

表7 年內(nèi)降雨徑流與污染負(fù)荷相關(guān)系數(shù)Table 7 Correlation coefficient between annual rainfall,runoff,and pollution load

為進(jìn)一步分析不同水文條件下的年內(nèi)污染負(fù)荷情況，明確非點(diǎn)源污染變化規(guī)律，還選取了豐、平、枯典型水平年進(jìn)行污染負(fù)荷變化趨勢(shì)分析。

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

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

2.2 空間變化結(jié)果分析

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

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

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

由表9可發(fā)現(xiàn)，人口數(shù)量、降水量、坡度與產(chǎn)水量與非點(diǎn)源污染負(fù)荷的相關(guān)系數(shù)均在0.25～0.45之間，說(shuō)明在空間尺度上，TN、TP的污染輸出綜合受到以上多種因素的綜合影響。相比其他因素，總氮、總磷的污染負(fù)荷均與子流域泥沙產(chǎn)量關(guān)系密切，相關(guān)系數(shù)分別為0.80和0.95。說(shuō)明泥沙產(chǎn)出是非點(diǎn)源污染負(fù)荷的重要驅(qū)動(dòng)因子，總磷負(fù)荷與泥沙的相關(guān)性更強(qiáng)，這也與潘家口水庫(kù)流域的總磷主要以吸附態(tài)磷流失有關(guān)。

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

3 結(jié) 論

a.從非點(diǎn)源污染排放總量角度進(jìn)行分析，2017年研究流域TN非點(diǎn)源排放量為94 769.73 t，TP非點(diǎn)源排放量為16 591.29 t，其中最主要的污染源為農(nóng)村生活污染和化肥農(nóng)藥污染。進(jìn)一步對(duì)潘家口水庫(kù)流域內(nèi)污染物形態(tài)分析發(fā)現(xiàn)，TN主要以硝態(tài)氮與有機(jī)氮形式流失，TP主要以顆粒態(tài)無(wú)機(jī)磷形式流失。

b.分析污染物各水平年時(shí)空分布特征發(fā)現(xiàn)，相較于平水年，枯水年年內(nèi)的TN、TP負(fù)荷峰值基本沒(méi)有下降，而流量則下降接近50%，從而可以說(shuō)明枯水年汛期非點(diǎn)源總氮濃度更高，污染情況更為嚴(yán)重。在對(duì)潘家口水庫(kù)流域內(nèi)污染進(jìn)行治理時(shí)，更應(yīng)關(guān)注枯水年下的污染削減情況。在空間尺度上，TN、TP負(fù)荷量從大到小排序依次為耕地、林地、草地、荒地，說(shuō)明農(nóng)業(yè)生產(chǎn)與農(nóng)村生活污染是潘家口水庫(kù)流域重要的非點(diǎn)源污染源。

c.對(duì)污染影響因素進(jìn)行分析發(fā)現(xiàn)，豐水年、平水年TP與降水的相關(guān)性系數(shù)分別為0.89、0.91，均高于當(dāng)年與徑流的相關(guān)系數(shù)(0.80、0.84)，說(shuō)明相較于徑流，TP污染受降水因素影響更大，更易受沖刷影響。TN、TP的污染負(fù)荷均與子流域泥沙產(chǎn)量關(guān)系密切，相關(guān)系數(shù)分別為0.80和0.95，說(shuō)明泥沙是流域非點(diǎn)源污染負(fù)荷的重要驅(qū)動(dòng)因子。