北运河上游非点源污染负荷模拟与最佳管理措施评估研究
时间：2023-04-20 23:35:02 来源：千叶帆 本文已影响人
宋卓远，李华林，于佩丹，谢晨新，张守红,2,3*
1.北京林业大学水土保持学院，北京 100083
2.山西吉县森林生态系统国家野外科学观测研究站，山西 临汾 042200
3.北京林业大学，北京市水土保持工程技术研究中心，北京 100083
随着社会发展，工业化和城市化进程加快，水污染问题不断加剧[1].非点源污染是水污染的主要来源之一，而氮、磷是非点源污染的重要组成部分，分别占我国水污染的81%和93%[2-4]，治理非点源氮、磷污染是水污染防治的重点.由于非点源污染具有随机性、广泛性、滞后性、潜伏性和隐蔽性等特点[5-6]，使其定量监测、模拟和治理十分困难.准确估算流域氮、磷负荷，定量解析其空间分布特征，并选择适宜的最佳管理措施，对流域水污染治理有至关重要的意义[7-9].
目前，国内外非点源污染模拟研究主要采用经验模型和物理模型模拟的方法[10-11].经验模型基于统计数据或因果分析，建立输入数据和输出数据之间的经验关系，模拟流域非点源污染负荷或风险[12].例如，输出系数模型(ECM)[13]、农业非点源污染指数(APPI)模型[14]和潜在非点源污染风险指数(PNPI)模型[15]等.这类模型数据需求量较少、计算及操作简单，能快速评估区域非点源污染负荷或风险[16].但经验模型无法定量模拟非点源污染产生、迁移和转化的物理过程[17].物理模型是根据污染物产生和迁移转化机理构建的数学模型，如 SWAT模型[18]、AGNPS模型[19]、ANSWERS模型[20]、SWMM模型[21]等.这类模型考虑了污染过程的内在机理，能进行长时间序列模拟，且模拟精度较高[22-23].其中，SWAT模型因能充分考虑水循环及物质循环过程、充分利用遥感数据被广泛应用[24-25].例如：向鑫等[26]在清溪河流域建立SWAT模型，并解析该流域非点源污染负荷时空分布规律；
Briak等[27]在摩洛哥北部卡拉亚流域，使用SWAT模型评估了BMPs对沉积物的影响.基于SWAT模型模拟结果，能够对流域非点源污染关键源区进行识别，有利于实现非点源污染的针对性治理.此外，通过SWAT模型自带的子流域输入编辑模块，还能实现不同最佳管理参数设置及布设，对措施布设前后模型输出结果比较能实现最佳管理措施生态效益的评估.因此，该研究选择SWAT模型对非点源污染负荷进行模拟.
最佳管理措施(BMPs)是指能有效控制点源和非点源污染的措施或措施组合，可分为工程措施和非工程措施[28-31].对BMPs生态环境效益进行评估是判断措施可行性的前提.例如，Uribe等[32]筛选出了对里奥格兰德湖上游流域氮、磷污染削减效果最好的BMPs；
Himanshu等[33]评估了BMPs控制污染物及泥沙流失的效果；
丁洋[34]评估了不同单项和组合BMPs对妫水河流域总氮、总磷的削减效果.以往研究通过对不同管理措施进行生态效益评价，筛选出适合治理流域非点源污染的BMPs，实现了非点源污染高效治理.由于受地形、土壤、土地利用、降水和人类活动等因素影响，流域内的不同区域对流域污染总负荷的贡献不同.同时，最佳管理措施的选择和布设也需要投入大量的人力、财力、物力.因此，不考虑污染负荷的空间分布特征，在流域内全面铺开式的布设非点源污染治理措施，会导致资源的浪费.为进一步提高非点源污染治理效率，使资源调配更合理，就需要对污染严重且范围相对较小的区域即关键源区进行优先治理[35-36].然而，当前对于BMPs的研究多集中在全流域BMPs筛选和效益评价[37-38]，缺乏针对非点源污染关键源区不同BMPs生态效益评价的研究.北运河是北京市重要的排水通道和连接京津冀的重要生态走廊.但近年来流域季节性断流频发，河流自净能力退化，水质不断恶化.因此，加强北运河上游非点源污染关键源区治理对北运河流域的水质改善至关重要，同时也是支撑京津冀生态文明建设的重要工作之一.
该研究基于北运河上游流域土地利用类型、DEM、土壤类型、气象、水文等数据，利用SWAT模型对北运河上游流域总氮、总磷负荷进行模拟，解析2019年北运河上游总氮、总磷负荷空间分布特征，采用单位负荷指数法识别研究区非点源污染关键源区，并评估关键源区布设不同BMPs的总氮、总磷削减率，以期为评估筛选BMPs提供科学参考，为北运河上游流域非点源污染防治提供科学依据.
1.1 研究区概况
北运河是海河的支流之一，发源于北京市海淀及昌平西部山区，向南流经通州区，最后汇入海河.北运河上游由南沙河、北沙河及东沙河三条支流组成，上游流域产流汇集于沙河水库.北运河上游流域地形以山地丘陵为主，地势呈西北部高、东南部低的特征，且起伏较大(见图1).研究区气候类型为温带大陆性季风气候，夏季高温多雨、冬季寒冷干燥.多年平均降雨量为500～600 mm，且年内分配不均，降雨量主要集中在6—9月[39].研究区土壤类型以褐土及潮土为主，土地利用类型主要为林地、耕地、城镇用地(见图2).由于研究区位于京津冀核心区，近年来流域人口快速增加，经济高速发展，这也导致了流域水污染问题加剧，其GB 3838—2002《地表水环境质量标准》劣Ⅴ类水的占比在20%以上[40].
图1 北运河上游的地理位置Fig.1 Location of the upper Beiyun River Basin
图2 北运河上游的土地利用类型Fig.2 Land use types of the upper Beiyun River Basin
1.2 研究方法
1.2.1 SWAT简介
SWAT(soil water and assessment tools，SWAT)模型是由美国农业部农业研究中心开发的分布式水文模型，用于模拟流域的径流、泥沙以及污染物的产生、迁移和转化过程[41-42].该研究根据DEM将北运河上游流域划分为90个子流域，再根据土地利用、土壤和坡度等数据将该流域划分为583个水文响应单元(hydrological response unit，HRU)，最后输入气象、农业管理措施和点源污染等数据构建SWAT模型.构建的SWAT模型对北运河上游流域2017—2019年的总氮、总磷负荷进行模拟.该研究使用的主要数据及其来源见表1.总氮、总磷污染强度的计算公式：
表1 基础数据及其来源Table 1 Sources of basic data
式中：[TN]为总氮的污染强度，kg/hm2；
ORGN为有机氮的污染强度，kg/hm2；
NSURQ为地表径流硝酸盐的污染强度，kg/hm2；
LAT_Q_NO3为侧向流硝酸盐的污染强度，kg/hm2；
GWNO3为地下水硝酸盐的污染强度，kg/hm2；
[TP]为总磷的污染强度，kg/hm2；
ORGP 为有机磷的污染强度，kg/hm2；
SOLP为可溶性磷的污染强度，kg/hm2；
SEDP为无机磷的污染强度，kg/hm2.
1.2.2 SWAT模型率定及验证
该研究采用SWAT-CUP软件对SWAT模型进行率定与验证.首先对参数进行敏感性分析，筛选出与径流相关性强的参数，然后基于沙河闸2018年日入库流量数据对径流进行率定，并采用决定系数(r2)和纳什系数(NSE)来评估径流模拟的精度.当r2＞0.6且NSE＞0.5时，可认为模型模拟精度达到要求[43].
式中，Oi为i时刻的实测值，O为实测值的平均值，Pi为i时刻的模拟值，P为模拟值的平均值.
由于缺少完整的长时间序列水质数据，基于北京市水科学技术研究院提供的北运河2018—2019年每月水质监测数据，选择海淀区温泉镇的东埠头沟(79号子流域)为验证点，对比分析了总氮、总磷负荷实测值与模拟值的差异以验证SWAT模型污染物负荷模拟精度.
1.2.3 关键源区识别
基于SWAT模型2017—2019年总氮、总磷输出结果，采用单位面积负荷指数法计算子流域总氮、总磷负荷强度，并划分为不同总氮、总磷污染强度等级，根据分级结果识别关键源区[44-45].单位面积负荷指数法是关键源区识别的常用方法[46]，其计算公式：
式中：Pi为子流域i在某时段内的污染负荷强度，kg/hm2；
PTi为子流域i在某时段内产生的污染负荷量，kg；
Ai为子流域i的面积，hm2.
1.2.4 管理措施设置及评估
结合北运河上游流域的农业及自然地理环境，该研究选择化肥减施、植被缓冲带和河道植草3种管理措施进行生态环境效益模拟评估.化肥减施措施只布设在农业用地，植被缓冲带和河道植草措施布设在所有土地利用类型上，具体措施描述及参数设置见表2.管理措施的生态环境效益多以污染物的削减率来表示，削减率计算公式见式(6).
表2 最佳管理措施情景设置Table 2 BMPs scenarios setting
式中：R为污染物削减率，%；
LBAS为基准情境下模型输出的总氮或总磷负荷，kg；
LBMPs为布设管理措施后模型输出的总氮或总磷负荷，kg.
2.1 模型率定与验证
参数敏感性分析结果及参数取值如表3.径流量对径流曲线参数、地下水延迟系数、浅水层补给参数、植被吸水补偿系数、土壤饱和水力传导度、河道曼宁系数、植被截留量等参数较为敏感，表明这些参数对径流模拟精度影响较大.图3为2018年北运河上游沙河闸径流实测值与模拟值比较结果，模拟的径流与实测的径流较为一致，r2和NSE系数均在0.63以上，表明构建的SWAT模型能较好地模拟北运河上游流域的径流特征.表4为2018年79号子流域总氮、总磷负荷实测值与模拟值的对比结果，SWAT模型模拟的总氮、总磷负荷与实测值一致，误差均小于10%，模拟精度较高.总体上，SWAT模型能够准确地模拟北运河上游流域的氮磷负荷，可应用于北运河上游流域氮磷污染模拟研究.
表3 SWAT模型参数敏感性分析及取值Table 3 Parameter sensitivity analysis and calibration of SWAT model
图3 2018年北运河上游沙河闸径流实测值与模拟值比较Fig.3 Comparison between measured and simulated runoff at the Shahe Floodgate Station of the upper Beiyun River in 2018
表4 2018年79号子流域总氮、总磷负荷的实测值与模拟值Table 4 Measured and simulated TN and TP load of 79 subbasins in 2018
2.2 氮磷污染负荷空间分布特征分析
北运河上游2019年总氮、总磷负荷的空间分布特征如图4所示.该研究采用自然断点法对总氮、总磷负荷进行等级划分，结果(见图4)显示，2019年北运河上游产生的总氮、总磷负荷分别为126 444.22和12 394.76 kg，空间分布呈现东南高、西北低的特征.总氮、总磷负荷较高的子流域主要分布在南邵镇、沙河镇、百善镇、昌平镇、马池口镇等地区.这些地区主要土地利用类型为城镇用地、耕地和农村居民点，地势平坦、人口密集且活动频繁，导致氮、磷污染物更加容易产生和积累.而总氮、总磷负荷较低的子流域主要分布在大庄科乡、井家庄乡、八达岭镇、雁翅镇等地区.这些地区主要土地利用类型为林地和草地，植被覆盖度高、人类活动少，因此氮、磷污染物不易产生，流失量也很少.此外，模拟结果还显示研究区西南部的流村镇地区总氮、总磷负荷也较高.这主要是因为该地区的土地利用以城镇用地和果园为主，频繁的人类活动和果园施肥也使该地区产生了较多的氮、磷污染物.
图4 2019年北运河上游总氮、总磷负荷的空间分布Fig.4 Spatial distribution of TN and TP load in the upper Beiyun River of 2019
2.3 关键源区识别
图5为总氮、总磷污染强度空间分布情况.由图5可知，总氮、总磷污染强度高的子流域集中分布在研究区东南部，其中19号、20号、32号等13个子流域为氮污染关键源区，产生的总氮负荷占全流域总氮负荷的34.56%.而磷污染的关键源区为14号、48号、64号等6个子流域，产生的总磷负荷占全流域总磷负荷的13.77%.综合氮、磷污染关键源区识别结果，确定研究区非点源污染关键源区为14号、19号、20号等17个子流域.这17个子流域面积仅占流域总面积的13.16%，产生的总氮、总磷负荷却分别占全流域的39.16%和38.10%，符合关键源区的基本特征，即较小的面积贡献较多的污染负荷.
图5 北运河上游总氮、总磷污染强度的空间分布Fig.5 Spatial distribution of TN and TP pollution intensity in the upper Beiyun River
2.4 最佳管理措施污染物削减效果评估
图6为不同植被缓冲带和河道植草措施总氮、总磷削减率.不同面积比植被缓冲带对总氮、总磷的削减率均在36%以上，其中面积比为1/5的植被缓冲带总氮、总磷削减率最大，分别为39.20%和40.37%；
面积比为1/60的植被缓冲带总氮、总磷削减率最小，分别为36.81%和38.93%.结果显示，随植被缓冲带面积比增加，总氮、总磷的削减率无明显提升(削减率增加均小于1%).其原因可能是，当植被缓冲带面积比超过一定范围时，大部分的污染物在缓冲带前端已被过滤截留，此时植被缓冲带面积比的增加对污染物削减率提升较小.而植被缓冲带作为工程措施不仅占用工程用地，措施建设和维修保护成本也较高，因此在措施布设过程中需要考虑当地实际的自然地理及社会经济条件.
河道植草对总氮的削减率在7.94%～19.47%之间，对总磷的削减率在26.29%～50.90%之间(见图6).随着河道植草长度增加，总氮、总磷的削减率呈增加趋势，说明河道植草措施的污染物削减率与植草河道长度呈正相关.结果显示，植被缓冲带和河道植草措施对总磷、总磷都有较好的削减效果.这是因为植被缓冲带布设在河流与陆地间，能通过吸收、阻滞、转化和过滤等作用增加对氮磷污染物的截留量和吸收量；
河道植草措施通过降低河流流速、减少河道冲刷，也能有效地削减非点源氮磷污染[47-48].但植被缓冲带和河道植草措施对泥沙沉积物的截留，可能会导致河道堵塞等问题，因此需要对植被缓冲带和河道植草措施进行定期清理维护.
图6 不同植被缓冲带和河道植草措施对总氮、总磷削减率的影响Fig.6 Reduction rates of TN and TP in different filter strip and grassed waterway
表5为化肥减施50%对总氮、总磷削减率的影响.化肥减施措施对总氮、总磷削减率分别为0.3%和3.0%，该措施对污染物的削减率差.这主要是因为化肥减施措施主要布设在农业用地上，而该研究区非点源污染关键源区农业用地占地面积小(仅占9.62%)，导致化肥减施措施对污染物削减效果差.
表5 化肥减施50%对总氮、总磷削减率的影响Table 5 Reduction rates of TN and TP by 50% of chemical fertilizer
该研究基于SWAT模型对北运河上游流域总氮、总磷负荷进行模拟，识别研究区非点源污染关键源区，并评估不同BMPs对关键源区总氮、总磷的削减效果.结果表明，植被缓冲带和河道植草对总氮、总磷的削减效果要优于化肥减施，这与前人的研究结果[49]一致.但在已有的研究中，多将植被缓冲带和河道植草布设在耕地或果园等农业用地[50-51]，没有考虑这两种措施在治理城市非点源污染方面的效益.而该研究将这两种措施布设在包括城镇用地在内的所有土地利用类型上，结果显示，不同植被缓冲带对总氮、总磷的削减率都在36%以上，河道植草对总氮、总磷的削减率分别在7.94%～19.47% 和26.29%～50.90%之间，说明这2种措施对北运河上游流域农业及城市非点源污染都有着很好的污染物削减效果.
以往研究[52-53]结果表明，BMPs是非点源污染治理最有效的方法之一，而其效益评估研究一直是研究的热点.最佳管理措施效益评估研究，主要分为生态环境效益评估和经济效益评估两方面.生态环境效益评估主要通过模型模拟的方法，例如：Liu等[54]基于SWAT模型评估单个和组合BMPs针对活性磷和总磷的削减效果；
杜颖恩等[55]基于SWMM模型评估不同LID设施空间格局对径流控制和雨水径流污染负荷的影响；
HSPF模型可通过SPEC-ACTION和BMPRAC模块评估BMPs对污染物的削减效果[56].以往研究多将不同BMPs布局在全流域范围内进而评估污染物削减率.这可能导致资源利用效率不高，未能将非点源污染的治理效率提升到最大，特别是在资源有限条件下.而该研究识别了研究区非点源污染关键源区，并重点评估了布设在非点源关键源区不同BMPs对污染物的削减效果.这不仅能有效提升非点源污染治理效率，还能使资源调配更合理.但受模型内措施库和数据量限制，本研究选择的最佳管理措施较少.因此，在以后的研究中需要进一步完善最佳管理措施数据库，为非点源污染治理提供更优的措施选择.该研究受数据条件限制，未能对BMPs的经济效益进行评估.对BMPs进行经济效益评估也是降低治理成本的有效途径.经济效益评估的方法比较多样，例如：Panagopoulos等[57]采用成本效益比评估BMPs经济效益；
Haas等[58]在考虑成本效益比的基础上进一步比较了不同BMPs的成本收益比；
欧美等国家通过建立BMPs成本-效益数据库为用户提供参考[59].但是，目前缺少评估经济效益的指标体系，对经济效益评估没有较为统一的标准.除经济效益外，本研究对最佳管理措施的设置仅考虑降低污染物削减率，而未考虑到影响最佳管理措施布设的其他因素，如地形地貌、土地利用等[60-61].因此，未来研究应进一步完善措施筛选布设体系，充分考虑最佳管理措施布设的多方面影响因素，科学合理地进行最佳管理措施的筛选布局.
此外，模型模拟与关键源区识别也存在一定局限性.物理模型建模过程中主要利用实测数据率定模型参数，进而降低模型不确定性[62].同时，对于不同BMPs污染物削减率的模拟结果也需要相关实测数据验证.然而，受监测及试验条件的限制，本研究缺少全流域的实测水质数据和实测的BMPs污染物削减率数据用于模型率定验证.尽管本研究基于现有水质数据从年总氮、总磷负荷上对模型模拟结果进行了验证，且实际值与模拟值之间的拟合度较高，但是为了能进行更长时间序列、更精确地模拟，仍需要加强北运河上游流域水质监测工作，并通过试验验证不同BMPs的污染物削减率，以便开展更精确地分析计算，提供更精确的科学参考.非点源污染的产生受地形、土壤、土地利用、降雨和人类活动等多种因素影响，这些因素可能会导致关键源区识别的不确定性.但基于SWAT模型的单位面积负荷指数法，未考虑降雨空间分布不均、流域地形等特点，在识别关键源区上存在不确定性.因此，在未来研究中还需要进一步加强关键源区识别方法的创新，提高关键源区的识别精度.
a) 2019年北运河上游流域产生的总氮、总磷负荷分别为126 444.22和12 394.76 kg，空间分布特征为东南部高、西北部低，城镇用地、耕地和果园是总氮、总磷负荷的主要来源.
b) 北运河上游关键源区分布在东南部17条子流域，占流域总面积的13.16%，产生的总氮、总磷负荷分别占全流域的39.16%和38.10%.
c) 1/5面积比植被缓冲带的总氮、总磷削减率最高，分别为39.20%和40.37%；
2 km河道植草的总氮、总磷削减率最高，分别为19.47%和50.90%；
由于北运河上游流域关键源区范围内农地面积较小(仅占9.62%)，化肥减施措施对污染物削减效果较低.
