蒲河流域氮污染负荷模拟及时空分布

韩博，周丽丽，范昊明，黄东浩

(沈阳农业大学水利学院，110866，沈阳)

非点源污染(面源污染)是指溶解性或固体污染物在降水和径流冲刷作用下汇入受纳水体而引起的水体污染，其主要来源包括土壤侵蚀、化肥与农药的使用过量、城市和公路径流、畜禽养殖和农业与农村废弃物等。近年来，随着农业活动的增加和生态植被的破坏，由农药、化肥、耕作等产生的农业非点源污染问题日益突出，尤其表现在水资源的污染与破坏。土壤养分吸附于泥沙进入水体或溶解在径流中流失，都会造成严重的水体非点源污染［1］。氮是作物生长必要的营养因子，同时也是引起水体富营养化和水质变差的根源因子，这已引起全世界的高度关注。资料显示，欧洲由于农业活动所造成的大量氮污染达到入海通量的60%［2］，在瑞典最南端的谢夫灵厄流域，由于农业活动而造成的河流运输氮量占流域总氮量的84% ～87%［3］，荷兰由农业非点源污染所提供的总氮量也占到水环境污染总量的60%之多［4］，在我国南方的三湖(巢湖、滇池和太湖)地区湖泊富营养化的现象日趋严重，其由农业活动所造成的总氮污染贡献率达到60% ～70%［5］。郭红岩等［6］对太湖流域非点源氮污染下水质影响进行定量化研究，研究结果显示:农田氮排放量占总量的72.7%，同样在我国北方的吉林省双阳水库，每年仅氮的负荷量也达到389.4t［7］。由此可见，对流域氮素污染负荷的研究，可为日后治理流域非点源污染提供规划依据。

SWAT模型是由美国农业部农业研究中心(USDA，Agricultural Research Service)开发的分布式流域水文模型，目前在农业非点源污染研究中得到了广泛的应用，模型已被应用于流域径流模拟、流量预测和非点源污染评价等方面，并取得了显著研究成果。万超等［8］应用SWAT模型对潘家口水库上游农业非点源污染负荷进行计算，研究以不同化肥施用水平进行模拟，模拟结果认为合理施用化肥对减少非点源污染有重大意义。王中根等［9］应用SWAT模型建立海河流域分布式水文模型，研究中对土地利用方式和土壤类型进行了重分类处理，并应用多个站点的实测数据进行结果校正与验证，其研究成果为SWAT模型在国内大型流域的研究应用提供范例。王林等［10］应用SWAT模型对东南沿海的晋江流域建立水文模型，通过年、月径流的模拟分析可知，SWAT模型适用于降水量丰沛的东南沿海湿润地区并且模拟效率受降水量及径流调节影响较大。谢元博［11］在泾河流域建立非点源污染模型，通过模拟分析可知污染负荷的产输出与当年的降雨径流量有密切关系。应用SWAT模型进行流域模拟及分析已在世界范围内得到了有效证实，但在我国东北三省地区，针对流域面积不足3 000 km2的中尺度流域农业非点源氮素污染研究实属罕见;因此，笔者以蒲河流域为研究区域，应用基于GIS的SWAT模型，进行蒲河流域氮素模拟并对其进行时空分布研究，为蒲河流域的生态恢复研究提供依据。

1 研究区概况

蒲河位于辽宁省沈阳市中部的产粮区(E 122°40′48″～123°56′32″，N 41°21′53″～ 42°4′15″)，发源于铁岭县横道河子乡，在棋盘山风景区上游望滨乡石砬子村进入沈阳市辖区，流经沈北新区、东陵、于洪、新民市和辽中5个县(市)区，于辽中县老观坨乡黑鱼沟村入浑河。河道全长205 km，流域面积2 248 km2。蒲河较大支流有九龙河和小浑河，主要排干有沈新辽、于家台、乌伯牛等。蒲河流域属温带大陆性季风气候，多年平均气温8.2℃。多年平均降水量647.4 mm，降水主要集中在7、8月，2个月降水量占年降水量的50%左右。蒲河流域上游为低山区，占流域面积的5.92%，区内山峦重叠，山高坡陡，但区内植被丰富，生长较好。中上游为丘陵区，占流域面积的13.48%，区内丘陵起伏，植被覆盖度较低;中下游为平原区，占流域面积的80.6%。该区地势平坦，土壤肥沃，是沈阳市的主要粮食产区。土壤类型主要有棕壤、草甸土、水稻土、风沙土和沼泽土等5类。

2 流域模型建立材料与方法

2.1 SWAT模型原理

SWAT模型是在20世纪后期由美国农业部农业研究中心推出的分布式流域水文模型，能够对具有复杂下垫面条件、多种土地利用方式和不同管理措施下的流域径流、泥沙及化学物质的产输出进行模拟计算。SWAT模型内主要含有水文过程子模型、土壤侵蚀子模型和污染负荷子模型3部分。

SWAT模型中水文过程子模型采用的水量平衡表达式为

式中:St为土壤最终含水量，mm;S0为土壤初始含水量，mm;t为时间步长，d;Rday为第 i天降水量，mm;Qsurf为第 i天地表径流量，mm;Ea为是第 i天蒸发量，mm;Wseep为第i天存在于土壤剖面地层的渗透量和侧流量，mm;Qgw为第i天地下水含量，mm。

SWAT模型中土壤侵蚀子模型采用的修正后通用土壤流失方程为:

式中:msed为土壤流失量，t;qpeak为洪峰径流，m3/s;Ahru为水文响应单元(HRU)的面积，hm2;KUSLE为土壤侵蚀因子;CUSLE为植被覆盖和管理因子;PUSLE为保持措施因子;LUSLE为坡长因子;SUSLE为坡度因子;WCFRG为粗碎屑因子。

SWAT模型能够追踪流域内氮磷随径流和泥沙的迁移和转化，土壤中氮磷形态转化是由氮循环和磷循环来控制实现的。

2.2 流域模型基础数据库构建

模型数据库构建按应用要求分为空间数据库和属性数据库2大类。空间数据库依托于GIS的支持，建立包括数字高程模型(DEM)、土壤类型、土地利用等空间数字图层。属性数据库则涵盖研究区水文、气象、水库、土壤属性、农业管理措施等文本数据，其中:研究区DEM格式为“.GeoTIFF”，来源于ASTER GEDM数据;土壤类型图格式为“.Grid”，来源于中国科学院南京土壤研究所与中国农业部共同推出1∶100万数字化土壤图;土地利用类型图格式为“.Grid”，来源于地球系统科学数据共享平台;水文、气象和水库数据格式为“.xls”，来源于国家气象信息中心和辽河流域水文年鉴;研究区土壤属性和农业管理措施数据格式为“.xls”，来源于辽宁土种志、辽宁土壤和辽宁省统计年鉴。根据SWAT模型需求，将经纬度坐标转化为大地坐标系，为便于对研究区域水文循环特征的描述，研究选用保留面积性质的Albers等积圆锥投影作为投影方法，并将空间数字图层坐标系进行统一设置。

2.3 流域模型应用

2.3.1 模型应用过程 应用蒲河流域数字高程模型(DEM)在SWAT模型中进行河网提取并利用实际水系对其进行校准修正，将蒲河流域内老观陀乡黑鱼沟村所在位置作为总流域出口进行模拟，同时将蒲河流域内的中型水库——棋盘山水库和各雨量站所在的位置作为子流域出口点，以方便水库定位和校准。经过模拟运行，将蒲河流域划分为23个子流域，划分子流域如图1所示。将土地利用类型图和土壤类型图及其模型所需的土地类型查询表和土壤类型查询表载入模型中。模型采用优势土壤(或优势土地利用)法进行HRU划分，通过设定阈值在子流域里划分出多个水文响应单元，模型结合蒲河流域实际情况，确定土地利用和土壤面积阈值均为4%，将蒲河流域划分为215个HRU。将各气象信息测站表导入模型中，模型则根据测站位置读取流域各气象信息(包括日降水数据、日最高气温、日最低气温、风速、相对湿度、太阳辐射等)。在模型中加入水库位置及其水库属性数据。根据流域土层实际状况及查询辽宁土种志和辽宁土壤，在模型中填入土层氮磷初始浓度。结合流域农业管理实际状况及辽宁省统计年鉴，在模型中输入农业管理措施数据。运行模型并利用实测数据对模型参数校准修正，提高模型精度，最终达到模型适宜性评价效果并应用校准模型对研究区氮污染负荷进行时空特征分析。

图1 蒲河流域子流域划分Fig．1 Sub-basin plans in Puhe river basin

2.3.2 模型参数校准 模型参数校准采用自动与手动调参相结合的方式，通过多次模拟运行，使模拟值与实测值拟合较好，实现参数最优化取值。模型校准按照先率定径流量再率定水质的顺序操作。利用多年实测径流资料对模型径流参数进行敏感性分析，率定分析结果显示:对模型输出影响较大的包括影响地下水与地表径流之间交换关系的基流α系数(ALPHA_BF)，影响径流过程的SCS径流曲线系数(CN2)和地表径流滞后时间(SURLAG)，影响蒸发过程的土壤蒸发补偿系数(Esco)，影响土壤水分过程的土壤深度(SOL_Z)，影响水层渗透的深蓄水层渗透系数(Rchrg_Dp)，影响河道汇流过程的河道有效水力传导系数(CH_K2)和主河道曼宁系数值(CH_N)，影响降雪和融雪过程的6月21日融雪系数(SMFMX)和12月21日融雪系数(SMFMN)。通过查阅相关文献和资料结合研究区实际状况，确定影响水质的参数为:氮渗透系数(NPERCO)，磷渗透系数(PPERCO)［12］。模型率定参数见表1。

表1 模型的率定参数Tab．1 Calibration parameters of SWAT model

应用实测数据对模拟数值进行校准拟合:校准期径流相对误差Re=-0.117 5，在±20%范围以内;R2=0.719 0＞0.6;Nash-Sutcliffe效率系数Ens=0.658 0＞0.5。验证期Re=-0.068 9，在 ±20%范围内;R2=0.882 0＞0.6;Ens=0.856 8＞0.5。模拟效果较好，总体反应了流域的实际情况。SWAT模型能够追踪流域内氮磷随径流和泥沙的迁移和转化，因此，在径流量拟合较好的前提下，氮的污染模拟也具有一定的可信性和参考价值。

3 结果与分析

3.1 流域氮素污染时间分布

表2 蒲河流域监测断面模拟结果输出Tab．2 Simulation results output of Puhe river basin monitoring sections

3.1.2 月际分布规律 通过对多年径流量与模拟氮污染负荷量分析可知，径流与TN有较好的相关性，R2=0.813 6。选取枯水年(2000年)、平水年(2004年)、丰水年(2010年)为代表年份，对不同水文年各月TN量进行模拟分析，TN模拟结果月输出值如表3所示，不同水文年径流与TN月拟合图如图2所示。分析图表可知:TN量在年内各月呈波动变化，1—4月的TN量较低，5—6月开始波动，7—9月达到最大，10—12月有所回落。丰水年(2010年)，降水量多，径流量大，径流峰值出现在8月份，TN量也因降水丰沛而呈增加趋势，并在9月份达到含量高峰;平水年(2004年)，径流峰值出现在汛期，在7月份达到最大值，TN量自8月份增长显著，9月份出现增长高峰;枯水年(2000年)，由于受1999年枯水年降水影响，加之该年降水量也极少，年内径流峰值出现在前半年，而TN却在8月份达到含量峰值。由此可见，TN流失具有明显的滞后性。

表3 TN模拟结果月输出Tab．3 Month results output of TN simulation103t/month

依研究区分布特点，自流域上游到下游选取5个监测断面，分别为望滨、道义、大河泡、马家窝棚和流域出口处(位置见图1)，结合1996—2010年降水量与模拟污染数值对流域氮素污染月际贡献率进行分析，监测断面月际平均径流、TN贡献率如图3所示，由图可知，径流与TN月际贡献率分配趋势大体一致，自流域上游望滨到流域出口，各监测断面径流与TN贡献率均在月际间出现高峰值:望滨在5—6月;道义在6—7月;大河泡在7—11月;马家窝棚在11—12月;流域出口在1—4月，自流域上游至下游各监测断面处贡献率高峰值从5—6月开始呈递推式发展，且望滨以上TN月均贡献率仅占全流域的1.93%;道义至望滨占全流域的11.09%;大河泡至道义占全流域的17.91%;马家窝棚至大河泡占全流域的25.56%;流域出口至马家窝棚占全流域的43.51%。流域TN贡献区集中在流域出口至马家窝棚水域。

3.2 流域氮素污染空间分布

图2 径流与TN月拟合图Fig．2 Runoff and TN month fitting figure

图3 监测断面月际径流、TN贡献率Fig．3 Month runoff and TN contribution ratio of monitoring section

图4 流域有机氮空间分布图Fig．4 Basin organic nitrogen space distribution

图5 流域NO3--N空间分布图Fig．5 Basin NO3--N space distribution

流域內TN平均质量浓度达到0.70 mg/L，子流域1的TN平均质量浓度为1.93 mg/L，根据国家环境保护总局颁布的地表水环境质量标准［15］判别，TN(以N计)质量浓度≤2 mg/L属于Ⅴ类水质，流域内子流域1属于氮素污染关键区，应特别注意统筹规划，与此同时，3、9、10号子流域TN浓度也达到了Ⅳ类水质，也是氮素污染的次要关键区域，对于此类地区应加强管理，防治流域氮素污染加剧，恶化流域水质。

4 结论与讨论

2)根据TN月模拟结果分析可知，TN量在年内各月呈波动变化，受汛期影响，TN量在5—6月开始波动，7—9月达到最大，10—12月有所回落，1—4月的TN量较低，对比流域径流月分布情况知，TN流失具有明显的滞后性。原因在于5—6月进入汛期，随径流量增加会加速吸附在土壤颗粒上氮物质迁移和转化，促使水体中淋溶作用和矿化作用发生发展，进而增加水体中TN含量。由此可见，径流量增加是导致TN量增加的直接原因。

3)自流域上游到下游选取5个监测断面，结合1996—2010年降水量与模拟污染数值对流域氮素污染月际贡献率进行分析知，自流域上游至下游各监测断面处贡献率高峰值从5—6月开始呈递推式发展，且流域TN贡献区集中在流域出口至马家窝棚水域。这是由于辽宁地区自5—6月降水量开始逐渐增多，随降水径流产生，流域中TN量增加，在地表径流、侧向流、渗流作用下氮素自流域上游至下游逐渐扩散增多，且流域出口至马家窝棚水域附近农耕地面积较大。由于农耕地产生的泥沙是化肥和农药的重要携带者［16］，这也为TN量积累提供有利条件，使该区成为TN主要贡献区域。

4)应用模型模拟数据，以子流域为研究单元分析流域氮素空间分布知，流域有机氮和NO-3-N整体呈现由东北向西南扩散，上游及中游地区含氮量高于下游地区，且随降雨量的增加而增加，子流域1的有机氮和NO-

3-N由于受地形坡度和降水丰沛影响而含量较多，子流域4和5则因林地茂盛而含量较少。根据流域內TN平均浓度和国家环境保护总局颁布的地表水环境质量判别标准知，流域内子流域1为氮素污染关键区，主要是由于地形坡度较大，随径流冲刷的泥沙中含有大量氮素，沉积在1号子流域，导致TN平均浓度过高。与此同时3、9、10号子流域TN平均浓度小于1号子流域，也是氮素污染的次要关键区域，同样应加以规划管理重视。

［1］ 叶芝菡，刘宝元，符素华，等.土壤侵蚀过程中的养分富集率研究综述［J］.中国水土保持科学，2009，7(1):124-130

［2］ Edwin D O.Control of water pollution from agriculture［M］.Rome:Food and Agriculture Organization of the United Nations Rome，1996:40

［3］ 苑韶峰，吕军，俞劲炎.氮、磷的农业非点源污染防治方法［J］.水土保持学报，2004，18(1):122-125

［4］ Paul C M B.Nutrient Emissions from Agriculture in the Netherlands，Causes and Remedies［J］.Water Science and Technology，1996，33(4/5):183-189

［5］ 蒋鸿昆，高海鹰，张奇.农业面源污染最佳管理措施(BMPs)在我国的应用［J］.农业环境与发展，2006，23(4):64-67

［6］ 郭红岩，王晓蓉，朱建国，等.太湖流域非点源氮污染对水质影响的定量化研究［J］.农业环境科学学报，2003，22(2):150-153

［7］ 李海杰.吉林省双阳水库汇水区农业非点源污染研究［D］.长春:吉林农业大学，2007:65-74

［8］ 万超，张思聪.基于GIS的潘家口水库面源污染负荷计算［J］.水力发电学报，2003(2):62-68

［9］ 王中根，朱新军，夏军，等.海河流域分布式SWAT模型的构建［J］.地理科学进展，2008，27(4):1-6

［10］王林，陈兴伟.基于3个站点校准与验证的晋江流域径流模拟［J］.中国水土保持科学，2007，5(6):21-26

［11］谢元博.泾河流域非点源污染特性分析与SWAT模型模拟［D］.西安:西安理工大学，2010:60-67

［12］陈媛，郭秀锐，程水源，等.基于SWAT模型的三峡库区大流域不同土地利用情景对非点源污染的影响研究［J］.农业环境科学学报，2012，31(4):798-806

［13］罗巧，王克林，王勤学.基于SWAT模型的湘江流域土地利用变化情景的径流模拟研究［J］.中国生态农业学报，2011，19(6):1431-1436

［14］孟晓云，于兴修，泮雪芹.云梦湖流域土地利用变化对非点源氮污染负荷的影响［J］.环境科学，2012，33(6):1789-1794

［15］易雯.《地表水环境质量标准》中氮、磷指标体系及运用中有关问题的探讨［J］.环境保护，2004(8):10-11

［16］杨巍，汤洁，李昭阳，等.基于SWAT模型的大伙房水库汇水区径流与泥沙模拟［J］.水土保持研究，2012，19(2):77-81