1 引言(Introduction)

随着人类活动加剧和极端天气的频繁出现, 水体中氮磷等营养物质含量持续增加, 造成全球水域富营养化问题严重(张恒等, 2012；Davidson et al., 2014).河流是流域氮磷营养盐输出的主要途径, 同时也是湖泊入湖污染物的主要来源之一, 研究河流的氮磷通量特征对水环境管理具有重要意义(毛战坡等, 2003；陈凌霄等, 2014).利用实测流量和水质数据计算河流污染物负荷是最接近“真值”的估算方法.目前河流污染物通量的研究方法主要包括流域模型法和统计模型法两大类(Alexander et al., 2002).其中流域模型法主要通过水文模拟和污染物迁移转化模拟来构建污染源排放与输出的响应关系, 常用的模型有SWAT、HSPF和AGNPS等(Hunter et al., 2008；刘博等, 2011).流域模型法对主要环境过程有较为详细的描述, 但对输入数据的质量和使用人员的专业经验有较高要求, 受资料和模型应用成本的约束, 很多情况下难以满足构建流域模型要求, 且模拟结果存在较强的不确定性(杨军军等, 2013；庞树江等, 2018).统计模型法主要包括平均法、流量负荷比例法和多元线性回归法(Dolan et al., 1981).其中平均法最为简单, 但计算误差较大.流量负荷比例法由Beale在1962年提出(Beale, 1962), 该方法使用污染物通量与流量的比例估算通量, 相对于平均法, 准确度较高, 但其仅能估算时间尺度较大时期的污染物通量, 无法在更细的时间尺度上估算通量变化.线性回归法则利用污染物浓度与流量和时间的多元统计回归关系模拟污染物逐日负荷, 该方法因对数据量的需求较小、模拟结果较为精确被广泛应用, 其中LOADSET(Load estimator)模型为其典型代表(Qian et al., 2007).水质监测频率对河流营养盐计算有显著影响(Webb et al., 1997；Johnes, 2007).Babitsch等(2020)指出小样本的使用可能导致污染物年负荷估算值偏低.Piniewski等(2019)研究发现, 高频水质监测数据能够改善模型模拟污染物负荷的性能, 尤其是在模拟总悬浮物和总磷负荷方面.Li等(2019)采用不同频率的水质监测数据评价了丰裕河流域的总氮年负荷特征, 发现两周一次的采样频率的评价效果相对较好.然而基于LOADSET氮磷通量估算的研究多采用低频水质数据(李娜等, 2017；陈炼钢等, 2019), 很少有研究考虑水质监测频率对其估算结果的影响(Toor et al., 2008), 特别是对云贵高原地区的入湖河流.

随着全球气候变暖, 极端气候事件发生的频率和强度发生显著变化, 对人类社会和生态系统产生重要的影响, 已经成为生态环境研究中不可忽略的部分(Moberg et al., 2005).已有研究表明气候变化会对地表水和地下水的水质有潜在影响(Lipczynska-Kochany, 2018).Tibby等(2007)通过分析澳大利亚湖泊水质与气候之间的关系, 发现水温、营养盐浓度与气温显著相关, 且降水也会对水质产生强烈影响.Salerno等(2018)通过模型模拟发现, 气候变化和城市化可能会导致河流水质严重恶化.在部分流域中气候变化对河流中氮磷负荷的影响超过了土地利用和人类活动(Morse et al., 2014).滇池作为云贵高原上最大的淡水湖泊, 具有供水、气候调节和文化娱乐等生态系统服务功能.由于氮磷营养盐长期超标, 滇池的富营养化问题十分严重, 蓝藻水华频繁爆发, 严重制约了区域生态环境的可持续发展(高伟等, 2019).宝象河是滇池的第二大入湖河流, 是滇池氮磷负荷的重要来源, 其TN和TP负荷分别占入滇池总量的10.66%和7.67% (杨文龙等, 2001), 研究宝象河污染物通量变化及其影响因素对滇池流域污染物总量控制具有重要意义.然而, 近年来对宝象河污染物通量的研究多基于低频的水质监测数据, 且对其驱动因子的研究也主要集中在人类活动(李娜等, 2012), 并未考虑低频采样是否适应于高原地区入湖河流的通量评估.

本研究基于逐日流量与周水质观测数据, 利用LOADSET模型建立了宝象河TN和TP通量的最优回归方程, 估算出TN、TP通量在不同时间尺度(日、季、年)的分布特征; 通过抽样部分水质数据, 评估了低采样频率对宝象河TN和TP通量估算的影响；最后, 基于宝象河逐日气温和降水数据, 利用RclimDex软件计算了流域13种极端气候指数, 分析了极端气候指数与宝象河氮磷通量变化的关系.

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

宝象河位于云南省昆明市滇池流域东北部(24°58~25°03′N, 102°41′~102°56′E；图 1), 属昆明古六河之一, 发源于官渡区东南部的老爷山, 干流在羊甫分洪闸分为南北两支分别汇入滇池, 河流全长41.7 km, 流域面积302 km², 占整个滇池流域的10.3%, 是入滇的第二大河流(郭怀成等, 2013).流域位于低纬度高原亚热带季风气候区, 年内温差相对较小, 多年平均气温14.7 ℃, 且干湿季分明, 年内降水主要集中在5—10月, 年平均降水量953 mm(苏斌等, 2019).宝象河干流设有水文站1个(在羊甫分洪闸以上), 具有连续的逐日水量观测资料, 为研究滇池入湖河流氮磷通量提供了良好的水文基础.

2.2 数据来源与样品采集

本研究使用了3种数据, 气象数据、河流水文数据和水质数据.其中, 气象数据采用了流域3个气象站(干海子、大板桥和宝象河)的逐日平均降水数据和昆明站的逐日最高气温和最低气温, 分别来源于云南省水文水资源局和国家气象科学数据中心(http://data.cma.cn/).河流水文数据采用了宝象河干海子水文站的逐日流量数据, 来源于云南省水文水资源局.气象水文数据的时间跨度为2018年1月1日—12月31日.水质数据(TN、TP)来源于课题组实验观测, 本课题组自2018年1月开展了在宝象河干海子水文站的水质监测, 采样频率为每周一次(周三或周四), 2018年共监测51次.研究中, 样品的前处理严格参照《地表水和污水检测技术规范》(HJ/T 91-2002), 利用采水器采样, 并将水样保存于聚乙烯塑料瓶中, 带回实验室分析水质.水质分析选用水杨酸钠法和磷钼蓝法, 利用CleverChem380全自动间断化学分析仪完成.流域水文、气象及水质监测点分布如图 1所示.

2.3 污染物通量估算模型 2.3.1 LOADSET模型

LOADSET(Load estimator)模型是美国地质调查局(USGS)基于FORTRAN编程语言开发的河流污染物负荷估算模型.该模型利用连续的流量数据和离散的水质浓度数据估算河流污染物通量(式(1)).根据AIC(Akaike information criterion)和SPPC(Schwarz posterior probability criterion)准则(Akaike, 1975；Schwarz, 1978), 从11种估算河流污染物通量的回归方程优选出拟合度最高的方程(Runkel et al., 2004), 用来模拟不同时间尺度(日、月、年)下河流的污染物负荷.

(1)

式中, 为瞬时污染物通量；为T时段的污染物通量；Δt为时间间隔；N为离散时间间隔数；Q和C分别为流量和污染物浓度随时间的变化数.

2.3.2 Beale比例法

Beale比例法作为一种污染物通量估算模型, 它是1962年由Beale (1962)提出的, Tin、Kendall和Stuart对其进行完善(Dolan et al., 1981), 其主要适用于流量数据较为完整而水质数据相对短缺的情况对污染物年通量的计算.Beale比例法如下所示.

(2)

式中, 为污染物负荷；μ_x为年平均流量；m_y为污染物测定日的平均污染物负荷；m_x为污染物测定日的平均流量；n为测定污染物浓度的天数.

其中S_xy和S_x²的计算公式如下：

(3)

(4)

式中, x_i为污染物测定日的流量；y_i为污染物测定日的污染物负荷.

2.3.3 平均法

平均法采用某一时间段内河流实测断面的污染物平均浓度和平均流量的乘积估算污染物通量(苏静君等, 2017), 具体计算公式如下：

(5)

式中, 为T时间段的污染物通量；Q_T为T时间段的平均流量；_T为T时间段的污染物平均浓度.

2.4 极端气候指数

本研究对极端气候指数的定义和计算采用世界气象组织“气候变化检测和指标”(Karl et al., 1999), 该标准被广泛应用于国内外极端气候的研究.世界气象组织气候委员会推荐使用的27个核心极端气候指数包括16个极端温度指数和11个极端降水指数(高妍等, 2014), 其计算可在RclimDex软件中完成, 该软件是加拿大气象中心基于R编辑器开发的用于计算多种极端气候指数的软件.由于部分计算降水和温度指数是基于年尺度计算出来的, 不适用年内数据的比较, 本研究计算了基于月尺度的DTR(平均温差)、RX1day(连续1日最大降水量)、RX5day(连续5日最大降水量)、TMAXmean(平均最高气温)、TMINmean(平均最低气温)、TN10p(冷夜日数)、TN90p(暖夜日数)、TNn(最低气温)、TNx(最低气温极大值)、TX10p(冷昼日数)、TX90p(暖昼日数)、TXn(最高气温极小值)、TXx(最高气温)13种极端气候指数.

3 结果与讨论(Results and discussion) 3.1 宝象河水文水质变化特征

宝象河2018年流量呈现单峰分布特征(图 2a), 具有显著的季节性差异(C_v=66.74%), 这与降水的变化趋势接近.在1—4月份降水少, 河流流量小, 且日流量波动较小.4月份后随着降水的增多, 流量开始逐渐增大, 8月份达到最高值, 最高可达9.14 m³·s^-1, 但日流量波动较大, 8月后流量逐渐减小, 于11月降到最低, 而12月河流在降水的作用下流量又呈现出上升趋势.

宝象河TN、TP浓度变化如图 2b所示, TN和TP的年均浓度分别为(4.60±0.83) mg·L^-1和(0.20±0.08) mg·L^-1.宝象河氮磷浓度具有显著的年内变化特征, 且TN和TP的浓度变化趋势接近, 两者存在显著相关性(r=0.56, p < 0.01).从箱线图可以看出, TN、TP浓度变化波动较大, 通过Mann-Kendall趋势检验(刘蓉等, 2016), TN和TP的统计量Z分别为-4.60和-5.90, 其绝对值均大于2.32, 通过了置信度99%的显著性检验, 表明宝象河的TN和TP的浓度在年内均呈现显著下降的趋势, 这与宝象河的流量时间分布密切相关, 降水在补充河流水量的同时, 稀释了水体中TN、TP浓度(涂安国等, 2010).

3.2 基于LOADSET模型氮磷逐日负荷模拟

本研究基于宝象河2018年的日流量数据和TN、TP浓度的周检测数据, 利用LOADSET模型估算了宝象河干海子断面氮磷的逐日负荷.根据AIC值最低的原则, 模型从11个污染物通量回归方程中优选出模拟宝象河TN、TP通量变化的最佳方程.LOADSET模型优选的宝象河TN、TP通量回归方程如下所示：

(6)

(7)

式中, L为污染物通量；Q为流量；lnQ=ln(Q)-ln(Q)中心值；dtime=十进制时间-十进制时间中心值；per代表用户自定义时间段(period).

结果表明, TN和TP通量回归方程的可决系数R²分别为94.02%和78.51%, 逐日通量的模拟结果如图 3所示.宝象河2018年TN和TP的日通量整体上均呈现出先上升后下降的趋势, 4月开始缓慢上升, 5月迅速上升, 8月出现峰值, 后逐渐下降, 与流量变化趋势相同.经过相关分析, 宝象河TN和TP的日通量与流量显著相关, 相关系数分别为0.98和0.95(p < 0.01).宝象河氮磷通量年内变化趋势与之前的研究结果一致(李娜等, 2012), 但本研究的年氮磷通量估算值(2018年)明显低于文献值(2008年), 一方面可能是由于流量监测点不同, 本研究的流量监测点位于中游的干海子, 而文献的流量检测点位于下游的宝丰村；另一方面可能是由于流域水环境治理和保护有所成效.2018年宝象河TN的年通量为270.49 t, 年内通量差异较大(C_v=51.11%), 其中夏季最多(110.31 t), 占年通量的40.78%, 秋季次之(73.23 t), 占27.07%, 春季和冬季TN的通量接近(43.44 t和43.54 t), 分别占16.06%和16.10%.TP的年通量为11.19 t(C_v=54.59%), 其中春、夏、秋、冬分别占年通量的18.04%、41.96%、23.05%和16.96%(2.02、4.70、2.58和1.90 t), 夏季TP通量最高.水体中的氮磷主要来自于城市生活和农业面源(徐晓梅等, 2016；苏斌等, 2019), 从氮磷通量年内分布可以看出, 夏季雨季来临, 城市面源和农业面源中的氮磷通过降水-冲刷作用进入河流, 造成宝象河氮磷通量增多, 强降水给河流带来的面源污染负荷掩盖了水量增加的稀释作用.

3.3 模型和水质采样频率对氮磷通量估算的影响

污染物通量估算结果的准确性及精度受估算方法及采样频率的影响.研究基于Beale比例法和年平均法分别估算了宝象河2018年TN和TP年通量, 并与LOADSET模型法估算结果进行比较(图 4).结果表明, Beale比例法对宝象河2018年TN、TP年通量的估算值分别为256.19 t和11.52 t, 较LOADEST的结果分别相差5.29%、2.98%；年平均法对TN、TP年通量的估算值分别为290.94 t和12.61 t, 与LOADEST的结果分别相差7.56%、12.68%.这两种方法与LOADSET模型的通量估算结果接近, 表明LOADSET模型对宝象河TN、TP通量估算结果相对可靠, LOADSET模型可用于宝象河TN、TP通量的估算.

采样间隔越大, 模型估算的通量误差越大(郝晨林等, 2012), 为研究不同采样频率对宝象河TN、TP通量估算的适应性, 研究以最接近真实值的高频采样污染物通量结果为标准值, 通过调整LOADSET模型的输入数据, 模拟了四种低频采样方式下(1次·月^-1, 采样时间分别在每月的第1周、第2周、第3周、第4周)宝象河TN、TP年通量(图 5), 并与高频采样方式下(4次·月^-1)年通量的估算值进行比较.研究结果表明, 低频采样1、2、3、4方式下, TN的年通量估算值分别为268.15、270.10、300.07、265.72 t, 相对于高频采样的误差分别为-0.87%、-0.15%、10.94%、-1.77%；TP的年通量估算值分别为10.89、12.44、9.87、10.79 t, 相对于高频采样的误差分别为-2.66%、11.12%、-11.82%、-3.63%.基于LOADSET模型低频水质采样的TN、TP年通量估算结果误差相对较小, 在研究宝象河TN、TP通量时可采用1次·月^-1的低频采样方式, 且在每月的第1周和第4周较为合适.

3.4 极端气候指数与氮磷通量的相关性分析

为评估主要极端气候指数与宝象河污染物通量的关系, 本研究利用RclimDex模型计算出2018年宝象河流域13个极端气候指数的月变化, 然后利用PASW Statistics 18分析了宝象河2018年12个月的TN、TP通量和极端气候指数的相关性(表 1).

相关分析结果表明, TN通量与低气温指数TMINmean、TNn、TNx(r=0.80、0.83、0.74, p < 0.01)高气温指数TMAXmean、TXn(r=0.60、0.69, p < 0.05)和降水指数RX5day(r=0.64, p < 0.05)呈正相关, 与平均温差DTR呈负相关(r=-0.87, p < 0.01).TP通量与低气温指数TMINmean、TNn、TNx(r=0.79、0.81、0.76, p < 0.01)高气温指数TMAXmean、TXn(r=0.61、0.66, p < 0.05)和降水指数RX5day(r=-0.69, p < 0.05)呈正相关, 与平均温差DTR呈负相关(r=-0.83, p < 0.01).研究发现, 宝象河TN和TP通量均与TMINmean、TNn、TNx、TMAXmean、TXn、RX5day、DTR显著相关, 但两者在相关性大小方面又有所差异, TN与TMINmean、TNn、TXn、DTR的相关性大于TP, 且TP与RX5day、TMAXmean、TNx的相关性大于TN.

温度升高和强降水会造成水体中营养物质和污染湖负荷增加(何霄嘉等, 2019).一方面当温度升高时会导致水体中的溶解氧减少, 水体在缺氧的还原性环境中, 沉积物中的氮磷更容易释放, 此外温度升高还会导致水体中微生物的活性增强, 也会促进河流沉积物中的氮磷释放到上层水体, 进而导致河流中氮磷负荷增加(Spears et al., 2006；Jiang et al., 2008).另一方面强降水通过冲刷地表侵蚀土壤, 不仅将城市和农业面源的污染物质带入河流, 降低水环境质量, 还会携带大量泥沙进入水体, 降低水体的纳污能力, 影响污染物的迁移转化(夏星辉等, 2012).综上所述, 宝象河TN、TP通量除受人为因素影响外(任玮等, 2015), 极端气候也是重要影响因素之一, 流域内TN、TP通量随连续降水增多、气温升高以及平均温差的减小而增多.

4 结论(Conclusions)

1) 2018年宝象河TN的年通量为270.49 t, 呈现显著的季节性差异, 其中夏季最多(110.31 t), 占年通量的40.78%, 春秋冬季分别占16.06%、27.07%和16.10%.而TP的年通量相对于TN较少, 为11.19 t, 其中春夏秋冬各个季节的TP通量分别占年通量的18.04%、41.96%、23.05%和16.96%, 夏季TP通量最多.

2) LOADSET模型、Beale比例法和年平均法对宝象河TN、TP年通量的估算结果接近, LOADSET模型估算结果可靠, 可用于宝象河TN、TP通量变化的研究；基于LOADEST模型的低频水质采样的氮磷估算结果与高频采样差异较小, 宝象河TN、TP通量估算受采样频率影响较小.

3) 宝象河TN、TP通量受连续5日最大降水量、平均最低气温、平均最高气温、最低气温、最低气温极大值、最高气温极小值和平均温差7种极端气候指数的显著影响, 其中平均温差与宝象河TN、TP通量呈现负相关关系, 在制定流域氮磷总量控制方案时应考虑极端气候的影响.