桃林口水库流域氮磷污染负荷对土地利用的响应分析

杨妮娟, 王晓云, 李建柱, 张 婷, 冯 平

(1.天津大学 水利工程仿真与安全国家重点实验室, 天津 300350; 2.中水北方勘测设计研究有限责任公司, 天津 300222)

1 研究背景

非点源污染是目前我国水环境治理的主要对象[1-3],由于其过程具有随机性和分散性,污染的形成与流失机理较为复杂,难以对其进行有效的管理与控制[4-5]。非点源污染受众多因素影响,其中土地利用方式直接影响着流域的水文循环与物质迁移过程,成为影响非点源污染的关键因子[6-7],由于近年来极端气候事件频发和人类活动加剧,土地利用方式的时空变化显著,因此明晰不同土地利用方式对污染流失的影响以及流域环境对土地利用时空变化的响应具有重要意义。

国内外学者采用不同的方法对非点源污染进行了研究。邓华等[8]采用田间试验法研究在不同的土地利用方式下石盘丘小流域非点源污染的流失变化,并针对总氮、总磷的流失特征分别提出了相应的治理措施,为流域的污染治理提供了科学依据。但试验法需要花费大量的时间和精力,通常只适合于面积较小的流域。付绍桐等[9]利用“源-汇”理论分析了太原市景观空间特征对非点源污染的影响,在“源-汇”理论的基础上,Huang等[10]结合最小阻力累积(minimum cumulative resistance, MCR)模型对黄河三角洲的非点源污染进行了动态模拟,指出耕地的污染风险等级为中度以上的面积占耕地总面积的47.57%。Srinivas等[11]采用多元回归方法探究了印度Ganges流域的土地利用类型与关键水质参数之间的关系。也有学者利用输出系数模型[12-14]、潜在非点源污染指数(potential non-point pollution index, PNPI)模型[15-16]研究了流域非点源污染对不同土地利用方式的响应,均指出耕地和城镇用地的污染流失较为严重。经验统计模型对数据要求不高,操作简便[17-18],但没有反映真实的水文过程和物质的迁移转化规律,不能及时响应水文要素的动态变化[19-20]。因此,探究土地利用方式的时空动态变化对非点源污染的影响,使用机理模型更为合适[21]。Shi等[22]、李明涛等[23]分别采用SWAT(soil and water assessment tool)模型、HSPF(hydrological simulation program fortran)模型模拟了景观格局变化对非点源污染的影响,发现农田的高度集中在一定程度上会加剧非点源污染负荷。宗敏等[24]以辽宁中部城市群为研究对象,利用SWMM(storm water management model)研究了城市化的快速发展带来的城市非点源污染变化。天目湖流域因茶园开发使土地利用方式发生了改变,赖正清等[25]利用APEX(air pollutants exposure)-SWAT集成模型研究了该变化对流域非点源污染负荷的影响,结果表明茶园用地的污染负荷贡献率逐年上升。

目前,研究土地利用方式对流域非点源污染的影响主要是通过对比分析污染负荷在单一土地利用类型下产出量的差异,而实际中流域的计算单元往往包含多种土地利用类型,并不都以某一类型为主,但考虑土地利用类型的不同组合对流域非点源污染的影响研究还较少。本研究以桃林口水库流域为研究对象,通过构建SWAT非点源污染负荷模型研究不同土地利用结构对各子流域氮、磷流失率的影响,并定量分析研究区在2014和2019年两种土地利用模式下非点源污染产量的变化,以期为桃林口水库流域合理的土地利用规划、水环境保护及秦皇岛市水资源管理提供科学依据。

2 研究区概况

选择位于秦皇岛市青龙满族自治县青龙河上游的桃林口水库流域为研究区。青龙河发源于承德市平泉县境内,是滦河的第二大支流。桃林口水库是国家和河北省“八五”“九五”期间的重点建设工程,其水质情况对当地的生活和发展都有着重大影响。桃林口水库流域控制面积为4 949 km2,占青龙河流域总面积的79%,其地理位置为118°37′～119°37′E,39°51′～41°07′N,海拔在92～1 841 m之间,地形总体上呈东北高西南低。桃林口水库流域地理位置及水系分布如图1所示。

图1 桃林口水库流域地理位置及水系分布

研究区属北温带半湿润大陆性季风气候,四季特征变化明显,昼夜温差大,无霜期较长,多年平均气温为24.7 ℃,年平均降水量为500～700 mm,降水主要集中在7—9月份,约占全年降水量的60%～70%,多年平均径流量为2.55×108m3。流域内的土地利用方式以林地为主,棕壤土和褐土是主要的土壤类型,主要农作物为小麦和玉米。根据实测资料,桃林口水库流域内多年平均总磷(TP)、总氮(TN)浓度分别为0.022和3.934 mg/L,根据《地表水环境质量标准》(GB 3838—2002),水库TP负荷满足Ⅱ类水质标准,TN负荷为Ⅴ类水质标准。桃林口水库是秦皇岛市重要的饮用水源地,流域内点源污染得到了有效控制,污染来源主要以非点源形式存在。

3 数据来源与研究方法

3.1 数据来源与处理

构建桃林口水库流域SWAT模型需要的基础数据包括数字高程模型(digital elevation model,DEM)、土壤类型、土地利用类型等空间数据,以及气象、土壤物理及化学属性、污染源以及水文水质等属性数据,数据来源及时空尺度见表1。

表1 模型数据库所需基础数据及其来源

研究区非点源污染来源主要分为农村生活、农药化肥和禽畜养殖3种类型。农村生活污染与畜禽养殖污染负荷产量较为固定,因此根据两种污染源的TN、TP和NH3—N的比例定义新的化肥类型,并以连续施肥的形式输入SWAT模型[26-27],施肥持续时间设置为1 a,施肥频次设置为每日1次[28],同时综合考虑各子流域的人口与面积来确定子流域的污染排放量。对于农药化肥污染,查阅年鉴发现桃林口水库流域主要农作物为小麦和玉米,小麦在3—4月份播种,7月份收获;玉米在6—7月份播种,10月份收获,因此结合两种作物的生长时间及所需施肥量,将施肥时间设置为每年的4—9月,施肥频次设置为每月1次。

3.2 研究方法

3.2.1 SWAT模型基本原理 SWAT模型可以模拟流域在不同情况下的水循环过程,能较好地预测人类活动对流域内水、泥沙以及化学物质的长期影响。水量平衡方程是SWAT模型的基础[29-30],其表达式如下:

(1)

式中:SWt为第i天的土壤最终含水量,mm;SW0为土壤初始含水量,mm;t为时间,d;Rday为第i天的降水量,mm;Qsurf为第i天的地表径流量,mm;Ea为第i天的蒸发量,mm;Wseep为第i天存在于土壤剖面地层的渗透量和侧流量,mm;Qgw为第i天进入主河道的地下水量,mm。

SWAT模型的水质模块来源于QUAL2E模型,能对流域中多种形式的氮和磷进行迁移转化模拟,氮和磷的转化分别由氮循环和磷循环控制。在SWAT模型模拟污染物的过程中首先对子流域进行离散,将其划分为属性一致的水文响应单元(hydrologic response unit, HRU),基于HRU对污染物的流失与迁移进行分析,再汇总到各个子流域,然后模拟污染物在河道及水库等水体中的迁移转化。

3.2.2 模型参数率定与验证 SWAT模型的率定与验证采用“先径流后营养物质”的顺序,选取桃林口水库上游双山子水文站2004—2020年实测逐月数据对径流模拟参数进行校准,其中2004—2005年为预热期,2006—2016年为率定期,2017—2020年为验证期。对水质参数的率定和验证采用桃林口水库坝上逐月实测水质数据,基于率定好的SWAT水量模块模拟2019—2020年TN和TP负荷,2019年为率定期,2020年为验证期。使用SWAT-CUP软件的SUFI-2算法对模型参数进行率定和验证,该算法计算效率高[31],并可以同时分析多个参数的敏感性。

采用决定系数(R2)和纳什效率系数(Ens)作为模型适用性评价指标,其表达式分别为[32]:

(2)

(3)

3.2.3 土地利用方式对非点源污染影响分析 考虑土地利用结构以及不同研究年份的土地利用变化研究土地利用方式对非点源污染的影响。基于流域2019年的土地利用数据将土地利用类型重分类为林地、耕地、草地、住宅用地和其他5种类型,并利用率定好的SWAT模型对流域的非点源污染进行模拟,分别统计各子流域2019年TN、TP负荷输出量和5种土地利用类型面积占比,根据子流域内5种土地利用类型的面积占比情况,将子流域划分为8种不同的土地利用结构,见表2。

表2 土地利用结构分类

以2019年为研究年份分析不同土地利用结构对桃林口水库流域TN、TP负荷的影响。考虑到降雨径流对非点源污染负荷的影响较大[35-37],选取降雨径流特征相似的2014和2019年作为研究年份,分析不同年份的土地利用变化对桃林口水库流域TN、TP负荷的影响。

4 结果与分析

4.1 子流域划分

基于桃林口水库流域DEM提取的河道水系,将研究流域划分为29个子流域,如图2所示。将模型中土地利用类型面积、土壤类型面积以及坡度阈值均设置为10%,按照相同的土地利用、土壤属性和坡度组合将研究区进一步划分为2 124个水文响应单元(HRUs)。

图2 研究区子流域划分

4.2 模型率定和验证

首先对模型参数进行敏感性分析,选取对径流和水质有重要影响的31个参数进行参数率定,率定参数敏感性及最优值见表3。SWAT模型模拟的率定期和验证期桃林口水库流域双山子水文站径流量和水库TN、TP负荷与实测值对比见图3,双山子水文站径流量模拟效果及桃林口水库TN、TP负荷模拟评价结果见表4、5。

表3 模型参数敏感性及最优值率定结果

表4 研究区双山子水文站径流量模拟评价结果

图3 研究区径流量及库内TN、TP负荷模拟值与实测值对比

由表4可知,率定期桃林口水库流域径流量模拟结果的R2=0.93,Ens=0.92;验证期R2=0.82,Ens=0.80,说明构建的SWAT模型能较好地模拟该流域的水文特征,从图3(a)也可以看出,模拟的各月平均径流量能够反映出该流域径流的年内丰枯变化,并且峰值的模拟效果较好。表5显示,率定期桃林口水库TN、TP负荷模拟结果的R2分别为0.95、0.94,验证期R2也均在0.85及以上,表明校准后的SWAT模型能较好地反映研究流域TN和TP的迁移转化情况。因此,构建的SWAT非点源污染模型对研究区的水文特征及TN和TP负荷的模拟具有较强的适用性。

表5 桃林口水库TN、TP负荷模拟评价结果

4.3 土地利用方式对非点源污染的影响

4.3.1 不同土地利用结构下子流域的非点源污染特征差异 利用SWAT模型对桃林口水库流域非点源污染进行模拟,得到研究区2019年各子流域的TN、TP负荷量,由于各子流域在空间上的污染物施用量存在较大差异,本研究以TN、TP负荷流失率(流失量/施用量)作为污染指标,以此反映土地利用结构对子流域污染流失的影响,研究区不同土地利用结构下各子流域TN、TP负荷流失率见图4。

图4 研究区不同土地利用结构下各子流域TN、TP负荷流失率(2019年)

从图4中可以看出,研究区土地利用结构为林地(Ⅰ)的子流域数量最多,TN、TP负荷在不同土地利用结构下的变化特征基本一致。在单一土地利用类型(Ⅰ、Ⅱ、Ⅲ、Ⅷ)中,子流域的TN、TP负荷流失率按从大到小排序为耕地(Ⅱ)>草地(Ⅲ)>住宅用地(Ⅷ)>林地(Ⅰ);在多种土地利用类型的组合(Ⅳ、Ⅴ、Ⅵ、Ⅶ)中,土地利用结构为林地-耕地-草地(Ⅶ)的子流域TN、TP负荷流失率最低,土地利用结构为耕地-草地(Ⅴ)的子流域负荷流失率最高;与相对应的单一土地利用类型的土地利用结构相比,林地的组合能有效降低子流域的污染流失率。

4.3.2 土地利用类型的年际变化对流域非点源污染的影响 2014、2019年研究区的土地利用转移矩阵见表6。分析表6可知,与2014年相比,2019年的林地和住宅用地面积有所增加,耕地和草地面积相应减小,2014年的土地利用类型中分别有25.26 km2的耕地、153.10 km2的草地转化成2019年的林地。表7和图5为2014、2019年研究区各土地利用类型面积及TN、TP负荷量。由表7可以看出,2019年桃林口水库流域TN、TP总负荷量较2014年明显减少,分别减少了34.43%、23.75%,说明林地面积的增加以及耕地和草地面积的减少对研究流域的TN、TP污染负荷有一定的削减作用。图5显示,研究区不同土地利用类型的单位面积产污量差异较大,林地的产污量最小,耕地的产污量最大;在研究年份(2014、2019年)中林地的平均单位面积TN、TP负荷量分别为10.73和0.40 kg/km2,而耕地的平均单位面积TN、TP负荷量分别高达87.69和4.28 kg/km2。

表6 2014、2019年研究区土地利用转移矩阵 km2

表7 2014、2019年研究区各土地利用类型面积与TN、TP负荷量

图5 2014、2019年研究区各土地利用类型面积及单位面积TN、TP负荷量

5 讨 论

桃林口水库流域SWAT模型对非点源污染的模拟具有较强的适用性,但从水质模拟结果中发现模拟值较实测值偏小,其主要原因可能为:(1)SWAT模型中需要对子流域土地利用和土壤类型的面积阈值进行设定,因此在模拟过程中会忽视对某些小面积区域的分析计算,但这些区域污染输出量可能会较大;(2)未考虑点源污染排放、大气沉氮等情况,因而造成了误差。

随着近几年环境保护力度的加大,退耕还林措施效果明显,2019年桃林口水库流域TN、TP输出量与2014年输出量相比有所减小。在单一的土地利用类型下,流域内林地的污染物流失量最小,草地次之,耕地最大[38-39],林地中树木的根系能固持土壤,具有较强的土壤保持能力,且能较好地过滤上方径流中的淤泥,通过改变径流对TN、TP负荷的运输过程,可以在一定程度上抑制污染物的迁移与转化[40-42]。桃林口水库流域的草地多为牧草地,所涉及的行政区养殖业较为发达,存在较高的牲畜养殖污染,因此以草地为主的区域污染流失也较严重,这一结论与Zeiger等[43]在James河流域的研究结果一致。与相对应的单一土地利用类型的土地利用结构相比,林地的组合能有效降低子流域的负荷流失率,说明林地对非点源污染具有较强的调控作用。

6 结 论

根据桃林口水库流域水文气象、水质、土地利用、土壤类型等资料,构建了流域SWAT非点源污染模型,对TN和TP负荷进行了模拟,分析了不同土地利用方式对流域非点源污染负荷的影响,主要结论如下:

(1)采用双山子水文站2004—2020年实测径流量和桃林口水库2019—2020年实测TN、TP负荷量数据,对构建的桃林口水库流域SWAT非点源污染模型进行参数率定和验证。率定期径流决定系数R2=0.93,纳什效率系数Ens=0.92,TN、TP负荷决定系数R2分别为0.95、0.94,纳什系数Ens分别为0.84、0.89;验证期径流决定系数R2=0.82,纳什系数Ens=0.80,TN、TP负荷决定系数R2也均在0.85及以上,说明构建的SWAT模型适合桃林口水库流域氮、磷负荷特征的模拟。

(2)不同土地利用结构下的子流域非点源污染特征差异较大,TN、TP负荷的分布规律相似,土地利用结构为林地的子流域污染流失率最低,耕地的污染流失率最高;与相对应的单一土地利用类型相比,林地的组合能有效降低污染流失率。

(3)不同土地利用方式的产污能力差异显著,林地的单位面积产污量最低,耕地最高;土地利用方式的年际变化对流域非点源污染也有较大影响,与2014年相比,2019年桃林口水库流域林地面积增加而耕地与草地面积减少,TN、TP总负荷量分别减少了34.43%、23.75%。