1 引言(Introduction)

氮含量升高被认为是导致富营养化、低氧和有毒藻类等水生态系统问题的重要原因(Scavia et al., 2006; Lundholm et al., 2010; Moline et al., 2010; Li et al., 2014).目前, 我国的主要河流与湖泊面临着严峻的氮污染形势(Hong, 2012; 荣楠等, 2016; Tong et al., 2019; Zhang et al., 2019), 且随着经济社会发展有进一步恶化的可能(Gu et al., 2015), 因此, 如何有效控制氮淋溶流失成为我国乃至世界氮素营养盐管理的重要挑战.氮的环境淋溶流失过程受到人为氮输入负荷强度、气象条件、土壤质地和土地利用等多种因素影响(Howarth et al., 1996; Chen et al., 2014; Huang et al., 2014), 对参数和计算方法要求较高, 这也是氮淋溶流失评估的难点.目前, 在区域尺度上针对氮淋溶流失评估的方法主要包括流域模型、统计模型和输出系数模型3种, 其中, 流域模型是基于水量、能量和物料平衡原理模拟氮在环境中的迁移转化过程, 需要大量的实测水文、泥沙、水质数据率定参数, 因此, 对基础数据要求高, 且计算过程十分复杂, 在大区域尺度上应用存在困难, 目前常用的流域模型如HSPF模型(Filoso et al., 2010)、SWAT模型(闫雪嫚等, 2017)、MIKE模型(闫雪嫚等, 2017)等均有氮淋溶流失模拟功能; 相对于流域模型, 统计模型对氮通量的计算多是基于经验公式, 优点是计算过程相对简单, 数据需求较少, 代表性的模型有Global NEWS模型(Mayorga et al., 2010)和NANI相关法(Swaney et al., 2012), 但应用该类模型主要模拟的是河川径流中的氮输出量, 而不是氮的淋溶流失量; 第3种方法是输出系数法(刘瑞民等, 2008), 该方法在计算上最为简单, 但存在的主要问题是未充分考虑氮淋溶流失的复杂过程, 计算结果的不确定性较高.在较大的空间尺度上, 目前尚缺乏直接核算氮淋溶流失量特别是非点源淋溶流失量的方法, 不利于区域氮污染的有效管理.长江经济带覆盖11个省市, 面积占我国国土面积的21%, 人口和生产总值均超过全国的40%, 是我国经济社会发展的重点地区.长江流域的主体位于长江经济带内, 相关研究表明, 长江水体的氮浓度和通量处于上升趋势(Li et al., 2014; 李茂田等, 2014; Wang et al., 2014), 长江经济带的氮淋溶流失控制对长江及区域内其他水系的水环境质量改善具有重要意义.目前, 针对长江经济带的氮淋溶流失及其水环境影响的研究主要关注长江的氮输出负荷变化及其影响因素(Li et al., 2014; Wang et al., 2014; Chen et al., 2016; Liu et al., 2018), 研究表明, 流域的氮输入负荷、水文和气象条件是长江氮输出通量的重要驱动因子.由于长江经济带是基于行政单元提出的, 与长江流域在边界上不重合, 针对长江流域的研究方法不能直接应用于行政区域, 相关研究成果难以服务区域的氮素管理.因此, 开展长江经济带的氮淋溶流失与水环境影响评估具有十分重要的意义.

本研究从人为和天然氮输入角度构建区域氮淋溶流失模型, 首先核算长江经济带各地区的氮淋溶流失系数, 结合氮输入负荷评估结果, 计算长江经济带各区域氮淋溶流失量, 并分析研究区氮淋溶流失量的时空变化特征; 同时, 选择典型流域, 开展氮水环境淋溶流失与通量负荷关系研究, 验证氮淋溶流失量的计算成果; 最后, 构建氮淋溶流失的水环境影响模型, 评估长江经济带的氮淋溶流失效应.本研究提出的氮淋溶流失研究方法可为大尺度区域氮淋溶流失核算提供方法借鉴, 研究成果能够为长江经济带氮素环境管理决策提供指导和依据.

2 材料与方法(Materials and methods) 2.1 研究区概况

长江经济带的地理范围包括上海、江苏、浙江、安徽、江西、湖北、湖南、重庆、四川、贵州、云南等11个省市(图 1), 总国土面积为205×10⁴ km².2016年常住人口为5.91×10⁸人, 其中, 城镇人口为3.37×10⁸人, 占比57%.地区生产总值为33.72×10¹²元(当年价), 三次产业结构为8 : 43 : 49, 以第三产业和第二产业为主, 第一产业占比较小.长江经济带涉及110个地级及以上城市, 是我国经济社会发展的重点区域.

2.2 基于氮输入的淋溶流失量评估模型

氮从陆域输入到水环境, 中间需要考虑氮的流失和淋溶过程, 通常采用淋溶流失系数描述这一过程.因此, 淋溶流失系数是影响计算结果可靠性的重要参数.在灰水足迹的计算模型中, 提供了一种估算非点源氮淋溶和流失系数的方法.该方法考虑了影响氮淋溶与流失过程的8个因素(表 1), 根据每个因素的权重和评分, 计算出研究区的氮淋溶流失系数.由于该方法需要的计算参数相对容易获得, 且考虑的影响因素比较全面, 能够适应多种气象、土壤和人为干扰特征区域.本研究采用该方法估算非点源的淋溶流失系数.氮的陆域输入(Total Nitrogen Inputs, TNI)分为人为源和天然源, 本研究中人为源的输入负荷采用人类活动净氮输入模型(Net Anthropogenic Nitrogen Inputs, NANI)(Howarth et al., 1996)核算, 天然源输入(Natural Nitrogen Inputs, NNI)基于NPP与植被固氮的关系(Thornton et al., 2007; Koven et al., 2013)估算.

陆域的总氮输入量TNI计算公式为:

(1)

式中, TNI为总氮输入量(kg·km^-2); NANI为人类活动净氮输入量(kg·km^-2); NNI为天然氮输入量(kg·km^-2).

基于以上计算思路, 本研究提出了一种基于氮输入的区域氮淋溶流失量计算模型:

(2)

(3)

式中, NR为陆域氮水环境淋溶流失量(kg·km^-2); α为非点源污染物的淋溶流失比例, 取值范围为0~1;N_FERT为氮肥施用量(kg·km^-2); NR_HC为农村生活氮消费量(kg·km^-2); N_LC为畜禽养殖氮消费量(kg·km^-2); N_LP为畜禽产品氮含量(kg·km^-2); N_DEPO为大气NO_y沉降量(kg·km^-2); N_ANF为固氮农作物氮固定量(kg·km^-2); NNI为天然氮输入量, 主要指非农作物植被固氮量(kg·km^-2); RWTN为农村生活污水处理厂总氮削减量(kg·km^-2); β为点源污染物入河系数, 本研究取1.0;NU_HC为城镇生活氮消费量(kg·km^-2); IN为工业总氮排放量(kg·km^-2), 根据环境统计数据可直接计算; SWTN为城镇污水处理厂的总氮削减量(kg·km^-2), 根据污水处理厂水量和水质统计数据可直接估算; α_min为淋溶流失比例的最小值, 氮取0.01(Franke et al., 2013); α_max为淋溶流失比例的最大值, 取值为0.25(Franke et al., 2013); s为研究区域环境因子(土壤质地、大气氮沉降、气候)和农业实践(固氮作用、施肥量、通过收割和管理实践导致的氮去除等)的评分; W为评分的权重; i为影响因子.以上人类活动相关的各组分氮负荷详细计算公式与及其参数参见文献(高伟等, 2016).根据文献(Franke et al., 2013), α计算公式中的参数可参考表 1.

2.3 基于实测径流通量的氮淋溶流失反演模型

为验证氮淋溶流失模型的结果, 选取水文和水质数据充分的流域作为验证区域.由于氮污染物在河流迁移过程中会因反硝化、取水和水库滞留等作用发生损耗, 因而河流的氮径流输出通量小于淋溶流失量.流域氮淋溶流失量与径流通量的关系可用图 2表示.

氮在河流输送过程中发生反硝化作用、消耗性取水和水库滞留作用, 氮的流失负荷和径流输出通量负荷可采用如下公式计算(Dumont et al., 2005):

(4)

式中, NR为陆域氮水环境淋溶流失量(kg·km^-2); RE为河流径流输出通量(kg·km^-2); 下标N表示氮; L_den为氮在河流迁移过程中发生反硝化作用损失的比例; R_rem为消耗性取水导致的氮损失比例; D_din为氮在流域湖库中的滞留比例.

① 河流径流输出通量:流域控制断面的污染物通量计算需要河流流量与污染物浓度数据, 其中, 河流流量可根据水文站点实测值计算, 水质数据一般来源于环境监测站.由于水质与水量数据在监测频率上不一致, 导致通量计算存在较大的不确定性.目前基于实测的流量和污染物浓度进行通量负荷计算的常用方法有插值法、Beale比率法、外推法和多元回归法等(Dolan et al., 1981; Lee et al., 2016).根据本研究的数据获取情况, 选择插值法进行通量计算.

② 河流输送过程参数:湖库滞留作用导致的氮损失比例D_din可采用如下公式估算(Seitzinger et al., 2002):

(5)

式中, DEPT为湖库平均深度(m); RT是水力停留时间(a); i为流域中的湖库的编号; Q为流域总流量(m³·a^-1); Q_i为湖库i滞留的流量(m³·a^-1).

当在流域的上下游存在梯级湖库时, 综合的D_dinc可以用如下公式计算:

(6)

式中, D_dinc为梯级湖库的综合滞留氮损失比例, n为梯级湖库的数量; 其他符号同前.

河流反硝化作用损失的氮比例L_den可采用公式(7)估算(Seitzinger et al., 2002; Dumont et al., 2005), R_rem可采用公式(8)估算(Dumont et al., 2005).

(7)

(8)

式中, A为流域面积(km²); Q_irr为流域灌溉耗水量, 为灌溉水量与退水量之差(m³·a^-1); Q_div为河道的跨流域调出水量(m³·a^-1); Q_nat为河流的天然径流量(m³·a^-1).

2.4 氮淋溶流失的水环境影响评估

为评估氮淋溶流失的水环境影响, 本研究引入水污染指数(Water Pollution Level, WPL)评价模型(Mekonnen et al., 2015), 定义为区域的灰水足迹与径流量的比值, 可用于表征水环境的承载能力.当WPL大于1时, 表明区域水资源对污染物的吸收能力小于污染物进入水体的量, 将导致水体污染.水污染指数的定义如下(Mekonnen et al., 2015):

(9)

(10)

式中, GWF为氮淋溶流失导致的灰水足迹(m³·a^-1); Q_act为区域的实测径流量(m³·a^-1); C_max是水体最大可接受氮浓度(或水环境质量标准)(mg·L^-1); C_nat为水体总氮浓度的天然背景值(mg·L^-1), 一般采用历史追溯法、区域对照法和地下水水质借用法取值, 如果难以使用这3种方法, 灰水足迹指导手册还提供了全球不同污染物质的推荐值, 根据数据可得情况, 本研究采用区域对照法, 使用研究区上游河段的水环境监测值.

2.5 数据来源

本研究所需的数据具体来源如下:①用于计算人为氮输入的人口数量、化肥施用量、畜禽养殖量和粮食产量等数据来源于长江经济带各省市的统计年鉴, 统计时间尺度为2011—2015年; ②用于计算天然氮输入的初级生产力数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn), 本研究收集了长江经济带2000—2010年逐年1 km空间分辨率的NPP数据, 该数据是采用GLM_PEM模型计算获得, 基于ArcGIS的空间统计功能得到各地级市的数据; ③大气氮沉降数据来源于文献(Jia et al., 2014), 该成果计算了全国范围的2000—2009年大气氮沉降的年均值; ④中国土壤侵蚀空间分布数据和中国土壤质地空间分布数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn); ⑤水文水质数据、长江流域各水文站点的流量数据来源于《中华人民共和国水文年鉴第6卷:长江流域水文资料》, 时间为2006—2015年, 长江经济带和流域水资源量和利用相关数据来源于同期地区水资源公报, 水质断面TN浓度来源于国控监测断面常规监测数据, 指标为TN年均浓度, 时间为2006—2015年; ⑥长江流域的水库参数来自于全球库坝数据库(Lehner et al., 2011), 湖泊参数根据公开发表文献(王铭玮等, 2011)调查获得.

3 结果与讨论(Results and discussion) 3.1 长江经济带氮输入负荷

根据总氮输入负荷评估结果, 2011—2015年长江经济带的年均TNI为17480 Gg·a^-1, 折合单位面积强度为8465 kg·km^-2, 是中国2010年的1.8倍(高伟等, 2014).输入强度最高的地区位于东部和中部的农业区域, 总体呈现高低交错的分布格局(图 3).长江经济带的氮输入以人为输入源为主, 90%的城市人为氮输入占总氮输入的75%以上, 天然源占比超过50%的只有四川省阿坝藏族羌族自治州和甘孜藏族自治州及云南省的迪庆藏族自治州和怒江傈僳族自治州.从人为氮输入负荷看, 2011—2015年长江经济带的年均NANI值为15621 Gg·a^-1, 单位面积强度为7610 kg·km^-2, 占TNI的89.9%, 东部、中部和西部之间呈现显著的区域差异, 同时也表现出省内的高度空间异质性, 导致空间异质性的主要原因为氮输入驱动因子的差异, 如人口密度(Han et al., 2014).从单位面积NANI强度看, 安徽省铜陵市、湖北省鄂州市、江西省新余市等地区强度最高, 主要来源于化肥施用的贡献.由于大气沉降和天然植被固氮输入的计算时段与研究期不完全一致, 可能对计算结果带来一定的不确定性, 但由于大气氮输入和天然氮输入不是长江经济带的主要输入源(大气氮沉降和天然源输入分别占总氮输入的8.1%和10.6%), 且本研究的数据是基于多年平均结果.因此, 时期不同带来的计算误差相对较小.

3.2 氮淋溶流失因子空间分布

根据长江经济带的降水、土壤类型、大气沉降等特征, 得到影响氮淋溶流失的因子评分和氮淋溶流失系数.基于氮淋溶流失系数计算公式, 可以得到长江经济带的氮淋溶流失系数变化范围为0.119~0.161, 平均值为0.143.从空间分布看, 长江经济带的氮淋溶流失系数呈现显著的空间异质性特征, 长江以南地区氮淋溶流失系数高于北部(图 4).

3.3 氮淋溶流失量评估

基于氮的输入负荷和氮淋溶流失系数, 采用氮淋溶流失量计算公式得到长江经济带各地区的氮淋溶流失量(图 5).从总量来看, 2011—2015年长江经济带的氮淋溶流失量为3399 Gg·a^-1, 单位面积流失强度为1656 kg·km^-2, 占总氮输入负荷的19.4%, 占人为输入的21.8%, 这一比例与我国主要流域径流负荷与人为输入比值范围(18%~20%)(张汪寿等, 2014)较为一致.长江经济带的氮淋溶流失高风险区域位于江苏省、安徽省、湖南省和四川省中部等地区, 最高为7533 kg·km^-2.在氮淋溶流失量计算中, 点源的削减量是重要参数, 本研究收集了2014—2015年长江经济带城市污水处理厂削减量数据.2015年长江经济带的总体污水处理率达到92.6%, 各地区污水处理率达到87.7%~96.7%, 处于较高水平.2006—2015年期间, 上升最快的是贵州省, 从2006年的29.5%上升到95.2%.东部地区是削减量最高的区域, 最高的是江苏省、浙江省、上海市等, 且2014—2015年间所有地区的总氮削减增长率也十分可观, 最高的贵州省达到213%.从氮削减总量看, 2015年总氮的去除量达到50.69×10⁴ t, 较2014年增长了46.3%, 增长率高于污水处理量的增长速度, 说明近年来污水处理工艺改进力度较大.扣除污水处理厂的削减量后, 2014—2015年点源淋溶流失量出现下降趋势.相对于未扣除污水处理厂的情况, 点源的占比有所下降, 从33%下降到23%.此外, 点源占淋溶流失总量的比例也出现下降趋势, 从2014年的25%下降到2015年的23%.因此, 污水处理厂对点源削减的贡献较大, 在点源排放量上升的背景下, 扭转了点源淋溶流失量增长的趋势, 氮的总淋溶流失量也加速下降.然而, 由于收集到的污水处理厂数据序列长度较短, 仅覆盖2014—2015年2年, 难以反映2011—2013年的情况, 计算的年均点源削减量可能存在高估.

3.4 基于实测水文水质数据的氮淋溶流失反演与验证

根据长江流域43个水质站点和46个水文站点的总氮和流量观测数据, 获得各站点氮多年平均负荷(水质为2006—2016年平均, 流量为2006—2015年平均).从上游到下游, 随着流量增长, 各站点的氮负荷增加显著.氮的平均负荷为265 Gg·a^-1, 最小为0.43 Gg·a^-1, 最大为1490 Gg·a^-1.根据径流负荷与流域淋溶流失量的函数关系, 主要的计算参数包括消耗性取水系数、河流氮反硝化移除系数和水库滞留系数等.基于2006—2017年长江流域及西南诸河水资源公报, 2006—2017年流域的耗水率平均值为42.5%.综合考虑地表水供水量、耗水率和天然径流量(取多年平均水资源量), 得到长江流域总体和各地区的消耗性取水系数.2006—2017年长江流域二级水资源区地表水消耗性取水系数的变化范围为0.40%~40.29%.

根据河网氮反硝化移除系数计算公式, 可以估算出长江流域各地区的河网反硝化作用对氮负荷的去除作用, 各子流域河流的反硝化作用去除比例变化范围为10.1%~64.4%.根据数据库和文献提供的参数, 可以计算出长江流域的水库和湖泊对氮滞留系数.经计算, 长江流域各水库的氮滞留系数在3.3%~96.5%之间.

基于实测数据的氮淋溶流失和基于淋溶流失系数计算的氮淋溶流失量存在显著的相关性, 氮的线性相关系数为0.82(p < 0.01).因此, 本研究提出的从氮输入角度结合淋溶流失系数核算氮淋溶流失量的方法能够重现氮流失空间分布规律, 也表明本研究对氮的估算结果具有一定的可靠性.然而, 从线性关系的斜率可以看出基于实测反演的氮淋溶流失量明显大于基于淋溶流失系数估算的结果.可能原因有两个:①本研究估算的值偏小, 由于污水处理厂数据收集到2014和2015两个年份, 其余3年数据缺乏, 而近年污水处理量和处理效率变化较大, 使用近年数据可能高估点源削减负荷; ②基于反演的结果偏大, 由于流量和水质数据有限, 采用年通量核算方法的不确定性较高, 加之目前仅考虑大型湖库, 中小型湖库的影响被忽略, 这些都可能导致反演结果误差.由于影响氮淋溶流失的因素较多, 受当前数据和资料限制, 难以充分解析其原因, 为进一步认识长江流域氮淋溶流失的变化规律, 急需加强氮通量的精细化核算和驱动力研究.

3.5 氮淋溶流失量对水环境影响

根据氮水污染指数计算结果, 2011—2015年长江经济带的平均WPL为3.0, 氮的灰水足迹是本地水资源量的3倍, 已经超过水环境承载能力.从行政单元看, 统计的130个计算单元平均WPL为5.4, 变化范围为0.2~29.2, 长江经济带的氮污染程度呈现显著的空间差异, 总体趋势为东部较西部严重, 北部较南部严重.仅云南南部、四川西部和江西省及湖北省部分地区共9个地级行政单元的水污染指数小于1, 仅占总量的6.9%.长江经济带93.1%的地级单元氮水足迹均超过本地水资源量, 特别是安徽省和江苏省的北部地区, 氮的水环境超载程度十分严重.

4 结论(Conclusions)

1) 2011—2015年长江经济带的总氮年均输入强度为8465 kg·km^-2, 在地级单位上呈现高度的空间异质性.人为源是长江经济带氮输入的主要来源, 90%的城市人为氮输入占总氮输入的75%以上, 其中化肥施用的贡献最高, 控制化肥施用强度对长江经济带氮污染控制具有十分重要的意义.

2) 长江经济带总氮输入的19.4%被淋溶流失到水环境中, 平均淋溶流失强度达到1656 kg·km^-2, 非点源是氮淋溶流失的最大来源, 但点源的削减对区域氮淋溶流失的影响较大, 已经扭转了长江经济带点源淋溶流失量增长的趋势.

3) 典型流域的实测氮通量反演结果与淋溶流失核算量存在显著的线性相关关系, 表明基于氮输入核算的氮流失量具有一定的合理性, 但不确定水平仍然较高, 需要加强基于高分辨率的流量和水质监测数据对长江经济带的氮通量研究.

4) 长江经济带的氮淋溶流失量导致的灰水足迹已经超过地区水资源总量, 在地级行政单元上呈现显著的空间差异, 其中污染最严重的城市水污染指数达到29.2.长江经济带的氮素水污染污染格局总体趋势为东高西低, 北高南低, 长江经济带的氮污染控制必须充分考虑区域水环境承载力的高度异质性.