张 婷，高 雅，李建柱，孙博闻，康爱卿

(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

分析不同水平年下的污染负荷情况与产水量的相关性。可以发现平水年与丰水年氮磷污染负荷与产水量相关性较好，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]。总体来看，汛期是潘家口水库一年中污染负荷最严重的时期，若想有效控制流域内的非点源污染，对汛期的污染控制尤为关键。

表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。总体来看流域多年平均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，说明泥沙是流域非点源污染负荷的重要驱动因子。