暴雨径流对新安江入库总磷负荷量的影响

李慧赟，王裕成，单 亮，罗潋葱，朱广伟,5，许 海，史鹏程,

1. 中国科学院南京地理与湖泊研究所，江苏 南京 210008

2. 杭州市生态环境局淳安分局，浙江 杭州 311700

3. 杭州市生态环境科学研究院，浙江 杭州 310014

4. 云南大学，云南 昆明 650500

5. 中国科学院大学，北京 100049

我国是世界上水库最多的国家[1-2]. 根据2019年《全国水利发展统计公报》数据，我国已建成各类水库98 112座，总库容达8 983×108m3[3]. 随着经济社会的快速发展，人类对水资源的需求不断增加，水安全问题已成为国家重大安全问题[4]. 水库的主要功能也从发电、防洪、灌溉逐渐转变为城市供水. 控制水体营养盐水平特别是降低水体磷水平，以控制水库富营养化、消除蓝藻水华，成为我国主要的湖库富营养化控制策略[5]，这也是北欧、北美等发达国家在20世纪六七十年代经大量辩论得出的结论[6-7].

磷是构成生命有机体的重要元素，是湖库富营养化状况和蓝藻水华情势变化中最重要的营养盐[8]. 然而，在湖库富营养化治理和蓝藻水华防控工作中，对磷的控制并非易事[9]. 以长三角地区重要战略饮用水源地千岛湖为例，尽管浙江、安徽两省为响应习近平总书记关于重视千岛湖水资源保护的号召，加强点源污染负荷治理，截至目前黄山市累计关停污染企业220多家、整体搬迁企业90多家，但千岛湖仍出现水体总磷(TP)浓度年际波动较大、河流区和过渡区TP高值波动的现象[10]. 千岛湖外源TP负荷呈现较大的波动性. 这些问题均与全球变暖导致的我国东南季风区暴雨事件频率和强度增加有关[11]. 然而，目前针对暴雨径流对大型水库外源磷输入影响的定量核算研究相对较少，暴雨前后水库内部水动力条件发生的变化对水库水体磷赋存的影响机制尚不清楚.

综上所述，该文旨在分析近60年千岛湖流域暴雨量和频率的时间变化特征，构建千岛湖流域水文模型，结合水质高频自动在线监测数据，定量计算暴雨径流对千岛湖主要入库河道(新安江)TP输入的影响.研究结果能够揭示千岛湖水环境对水文过程变异的响应，为千岛湖水环境安全保护与管理提供理论与技术支撑.

1 材料与方法

1.1 研究区域

千岛湖(原名新安江水库)位于浙江省西部与安徽省南部交界的淳安县境内，是我国第一座自行设计建造的大型水电站拦蓄新安江干流而形成的深水水库[12]. 水库坝高105 m，库区水面纵长150 km，平均水深31.13 m，最大水深100 m[13]. 当设计水位(黄海标高)为108 m时，相应水库面积为573.33 km2，蓄水量为178.6×108m3. 千岛湖是典型的山谷特大型水库，新安江电站以上流域面积为10 442 km2. 最大入库河流新安江发源于皖赣交界处怀玉山脉，流经安徽省黄山市的休宁县和屯溪区等地，并于歙县街口镇汇入千岛湖，其多年平均入库流量为225 m3/s，约占水库多年平均总入库流量的63%. 千岛湖及其上游流域位于我国亚热带季风气候区，多年平均气温为17.3 ℃，多年平均降雨量为1 733 mm. 千岛湖降雨量季节性分布存在明显差异，春、夏季降雨量充沛，降雨量约占全年总降雨量的74%.

千岛湖最初以发电为主，兼顾防洪、灌溉等其他效益. 随着经济社会的快速发展，下游地区对供水需求和水库生态环境保护的要求日益提高. 2019年9月，千岛湖正式向杭州市供水，成为杭州市800多万人口的饮用水源地，其水质的安全保障事关重大.自1960年建库以来，千岛湖素以水质良好著称，并于2012年被纳入国家良好湖泊生态环境保护试点.然而，作为一个入库流量大、面源负荷高的山区水库，千岛湖水质存在较大的时空波动性[10]. 近年来，千岛湖水环境已发生显著变化，具体表现为透明度逐渐下降、TP浓度波动超标、蓝藻化趋势明显等. 千岛湖虽全年大部分季节水质稳定，全库85%以上监测点TP达到GB 3838-2002《地表水环境质量标准》Ⅰ～Ⅱ类水平，但街口(新安江入库断面，见图1)等个别监测点TP为Ⅲ～Ⅳ类水平，且受暴雨影响剧烈. 暴雨径流给千岛湖水环境质量保护带来较大压力.

1.2 数据来源

1.2.1气象水文数据

在中国气象数据网(http://www.nmic.cn)下载得到千岛湖流域各国家气象站(屯溪站和淳安站，见图1)1961-2021年逐日降雨量和蒸发量数据；通过《中华人民共和国水文年鉴》[14]摘录、整理得到自2001年以来千岛湖上游各水文站(屯溪站和渔梁站)逐日流量数据. 以上数据用于构建新安江模型.

1.2.2常规水质监测数据

图 1 千岛湖流域地理位置Fig.1 Location of Qiandaohu Reservoir Basin

为了核算千岛湖最大入库河流新安江(街口断面)外源TP周年负荷量，选取2020年5月-2021年4月为典型年，按照常规水质监测方法，于典型年的每月下旬，对街口断面的TP浓度进行监测. 水样分表层(水下1 m)、中层(水下15 m)和下层(水下28 m)3层采集，现场采用保温箱冷藏后及时带回实验室进行分析. TP浓度采用GB 11893-1989《水质总磷的测定 钼酸铵分光光度法》[15-16]测定. 取表层、中层和下层水样TP浓度的平均值作为该断面的TP浓度值，参与外源负荷量的计算. 以上数据用于计算街口断面外源TP的周年逐月负荷量.

1.2.3高频水质在线监测数据

通过街口水质自动监测站监测得到2020年1月1日-2021年5月16日高频TP浓度数据，监测频率为1条/(4 h). 获取TP浓度高频数据共计2 694条(缺测108条)，采用三次多项式方法进行插值补缺后将时间单位折算到日. 筛选得到街口断面与常规水质周年监测同时期(2020年5月1日-2021年4月30日)逐日TP浓度系列，与逐月人工监测数据进行对比，结果如图2所示. 由图2可见，插补校正后的逐日TP浓度序列具有较高的精度. 以上数据用于计算街口断面外源TP的逐日负荷量.

图 2 2020年5月1日-2021年4月30日街口断面TP浓度和降雨量的逐日变化Fig.2 Monitoring data of TP concentration and rainfall of Jiekou Section from May 1, 2020 to April 30, 2021

1.3 研究方法

1.3.1暴雨的定义

暴雨是一种短时间内产生较强降雨的天气现象.根据国家质量监督检验检疫总局、国家标准化管理委员会批准发布的GB/T 28592-2012《降水量等级》[17]，当12 h降雨量为30.0～69.9 mm或24 h降雨量为50.0～99.9 mm称为暴雨，暴雨之上又划分大暴雨(12 h降雨量为70.0～139.9 mm或24 h降雨量为100.0～249.9 mm)和特大暴雨(12 h降雨量≥140.0 mm或24 h降雨量≥ 250.0 mm)两个量级. 该文将暴雨、大暴雨和特大暴雨均列为研究范围，统称为暴雨.

1.3.2新安江模型的构建

该文采用新安江模型模拟千岛湖主要入库河道的水文过程. 该模型由Zhao所在课题组[18-19]于1973年编制千岛湖洪水预报方案时设计开发，其特点是采用流域蓄水容量分布曲线考虑产流面积内各点蓄水的不均匀性. 新安江模型的蒸散发计算采用三层模型，产流计算采用蓄满产流模型，总径流利用自由水蓄水库结构分为地表径流、壤中流和地下径流三部分，三水源均按线型水库结构计算河网总入流，河网汇流采用延迟滞时法[19-20]. 新安江模型包含4个计算子模块，共14个参数. 目前，该模型已经在国内外湿润、半湿润地区得到广泛应用[21-22].采用通用模式搜索(generalized pattern search，GPS)优化算法自动优选模型参数[23]. GPS算法是在特殊方向集上抽取目标函数，通过比较函数值的大小，找出下降方向进而解决最优化问题[24]. 该文以Nash-Sutcliffe效率系数(N SE)最大化作为参数优选的目标函数，并在率定过程中加入水量平衡误差( WBE) 线性约束[25]. NSE和 WBE的计算公式分别见式(1)(2)：

式中：Si为模拟流量，m3/s；Oi为 实测流量，m3/s；O为实测流量平均值，m3/s；n为实测值或模拟值的个数. 其中， N SE反映了模型的整体效率，其值越接近1，表明模型的适用性越高； WBE 满 足| WBE|≤5%.

1.3.3暴雨径流量的模拟方案

基于该文构建的新安江模型，提出了计算通过街口断面的暴雨径流流量的方法，具体步骤为：①计算暴雨前降雨产生的径流流量(Q0)，即将暴雨当天及其后日期的降雨量设置为0，运行模型计算街口断面流量；②计算暴雨发生后的径流量(Q)，即将暴雨当天降雨量设置为实际暴雨量，暴雨后的日降雨量设置为0，运行模型计算街口断面流量；③计算本次暴雨产生的径流量(Qs)， 即Qs=Q-Q0.

1.3.4暴雨径流携带外源TP负荷量的计算方法

千岛湖主要入库河流TP负荷量的计算方法[26-27]如下：

式中：W为通过街口断面的TP负荷量，t；C为TP浓度，mg/L；Q为河道流量，m3/s；t为时间，d；η为时间换算系数，取值为 8 6400×10-6.

2 结果与分析

2.1 模型率定与验证

由于街口断面缺乏历史实测流量数据，故分别采用千岛湖上游流域屯溪站和渔梁站逐日实测流量数据对新安江模型进行率定和验证. 其中，模型率定期为2001-2015年，验证期为2016-2020年. 屯溪站模拟流量和实测流量拟合程度较好，率定期 N SE为0.94， WBE 为 0.8%，验证期N SE 为0.93， WBE为-3.4%；渔梁站模拟流量和实测流量拟合程度稍差，率定期NSE 为0.91， WBE 为 1.6%，验证期 NSE 为0.90，WBE为-4.9%. 其中，屯溪站模拟流量与实测流量对比如图3所示. 以上结果说明，新安江模型在对千岛湖流域2个水文站所在断面的流量模拟中，均得到较高的模拟精度，模型在千岛湖流域有良好的适用性. 鉴于屯溪站模拟精度较好，因此采用该站实测数据率定所得参数作为模型最优参数，用于模拟街口断面逐日流量，即新安江入库流量.

图 3 千岛湖流域上游屯溪站模拟流量与实测流量的对比Fig.3 Comparison between simulated streamflow and observed streamflow at Tunxi Station of Qiandaohu Reservoir Basin

2.2 千岛湖流域暴雨时间变化特征分析

分别统计千岛湖流域1961-2020年不同量级降雨量及降雨频次，小雨、中雨、大雨和暴雨雨量分别为294.6、459.8、477.2和482.2 mm，分别占年降雨量的17.2%、26.8%、27.8%和28.1%；小雨、中雨、大雨和暴雨平均频次分别为100.8、28.9、13.9和6.2 次/a.采用Mann-Kendall秩次相关检验法[28-29]分别对千岛湖流域1961-2020年不同量级降雨量和降雨频次系列的趋势性进行检验，结果如表1所示. 近60年来，千岛湖流域暴雨雨量和暴雨频次均在0.001显著性水平上检验到上升趋势，暴雨量平均每年上升6.42 mm，暴雨频次平均每年增加0.08次. 小雨雨量在0.1显著性水平上检验出上升趋势，小雨雨量平均每年上升0.40 mm，小雨频次无显著变化趋势. 另外，中雨、大雨雨量和降雨频次均无明显变化趋势.

表 1 千岛湖不同量级降雨量和降雨频次多年变化趋势Table 1 Multi-year variation trend of rainfall of different magnitude and rainfall frequency in Qiandaohu Reservoir

在水质周年监测时期(2020年5月1日-2021年4月30日)，分析千岛湖上游流域不同量级降雨量和降雨频次可知，典型年小雨雨量为331.7 mm，降雨103次；中雨雨量为527.0 mm，降雨32次；大雨雨量为559.0 mm，降雨16次；暴雨雨量为779.3 mm，降雨9次. 典型年暴雨雨量显著高于多年平均值，为多年平均暴雨量的1.6倍，暴雨次数与多年平均暴雨次数相比，增加42.9%. 从暴雨的季节性分配来看，夏季暴雨量占年总降雨量的90.7%，春季暴雨量占年总降雨量的9.3%. 以上结果可以看出，选择2020年5月1日-2021年4月30日作为典型年，研究暴雨径流对千岛湖街口断面外源TP负荷量的影响是合适的.

2.3 暴雨径流对新安江外源TP负荷量影响的定量核算

2.3.1基于常规水质监测的新安江外源TP负荷量核算

采用新安江模型模拟得到新安江街口断面逐日流量并折算到逐月，结合逐月常规水质监测所得街口断面TP浓度，计算逐月通过街口断面的TP负荷量，结果如图4所示. 由图4可知，2020年5月-2021年4月街口断面月均流量为304 m3/s，为多年平均流量的1.4倍. 月均流量最高值为1 264.1 m3/s (2020年6月)，最低值为16.7 m3/s (2021年1月). 街口断面TP浓度与流量具有较好的相关性(R=0.66). 典型年全年通过街口断面的TP负荷量为517.1 t，TP逐月负荷量变化范围为1.1～200.3 t，其中暴雨量集中的6月和7月TP入库负荷最高，6-7月TP入库负荷占全年TP入库负荷的69.6%，其次是中雨和大雨量集中的5月和9月. 可以看出，TP负荷量年内分布与河道流量年内分布具有较强的一致性，说明通过街口断面的TP负荷量受降雨量影响显著. 这是由于千岛湖地处亚热带季风气候区，降雨量年内分配不均匀，春季及梅雨季降雨量大，而秋冬季节干旱少雨. 雨量的波动性导致入库流量年内变化大，水体TP浓度明显呈现雨季高、枯季低的特征. 相应地，TP入湖负荷呈现典型的季节特征.

图 4 2020年5月—2021年4月街口断面逐月流量、TP浓度和负荷量模拟值Fig.4 Simulated monthly streamflow, TP concentration and TP loading at Jiekou Section from May, 2020 to April, 2021

2.3.2基于高频水质监测和新安江模型的暴雨径流定量核算

采用该文提出的暴雨径流量模拟方案，逐日模拟由暴雨产生的通过街口断面的径流量，结果如图5所示. 典型年通过街口断面进入千岛湖的总径流量为96×108m3，其中由暴雨产生的径流量为46×108m3，占总入库径流量的47.9%. 由夏季暴雨产生的径流量占年总径流量的43.8%，由春季暴雨产生的径流量占年总径流量的3.6%.

图 5 新安江模型模拟所得街口断面在实际降雨条件下与去除暴雨条件下的流量Fig.5 Simulated streamflow of Jiekou Section under actual rainfall conditions and removal of rainstorm conditions

2.3.3暴雨径流携带TP负荷对新安江总入库负荷量的贡献

图 6 千岛湖流域逐日降雨量、街口断面逐日模拟流量、逐日TP浓度监测值和逐日TP负荷量计算值Fig.6 Daily rainfall of Qiandaohu Reservoir Basin, simulated daily streamflow at Jiekou Section, daily TP concentration monitoring value and calculated daily TP load

采用新安江模型对典型年街口断面逐日流量进行模拟，结合街口水质自动监测站的TP浓度高频观测数据，得到街口断面TP逐日入库负荷量，结果如图6所示. 通过逐日流量模拟，可以看出街口断面在典型年最高流量为4 437 m3/s，出现在2020年7月7日，当日千岛湖流域降雨量达165.7 mm. TP日均浓度为0.048 mg/L，最高值为1.278 mg/L，出现在2020年7月8日. TP日均负荷量为4.1 t，最高值为438.6 t，出现在2020年7月8日，与TP最高浓度出现日期一致. 在整个典型年，TP全年总负荷量为1 506 t. 与基于常规水质月监测值的计算结果相比，采用高频监测数据计算所得TP全年总负荷量是采用逐月监测结果计算所得TP全年负荷量的2.9倍. 其中，在强降雨集中的6-7月，TP入库负荷占全年TP入库负荷的91.3% .

基于暴雨径流量核算结果和高频观测数据，计算由通过街口断面的暴雨径流所携带的TP负荷量.典型年通过街口断面的暴雨径流所携带的TP负荷量为1 046 t，占全年TP总入库负荷量的69.4%. 其中，由暴雨径流携带的TP负荷量为1 029 t，占全年TP总入库负荷量的68.3%. 由此可见，短时间、高强度的暴雨径流对TP入库负荷量的影响显著. 在充分考虑暴雨过程影响的情况下，对TP入库负荷估算所得数值远大于忽略暴雨过程影响情况下计算所得数值(人工水质监测通常避开暴雨天气). 这主要是因为，相对于中小降雨而言，暴雨冲刷易造成大规模土壤侵蚀，同时由于暴雨径流的携沙能力较强，可以将大量含磷污染物输送进入水库.

3 讨论

该文借助新安江模型和水质高频自动监测数据，对千岛湖街口断面外源TP负荷量进行了逐日计算，定量区分了通过该断面的暴雨径流量及其携带的TP负荷量，并分别估算二者对街口断面周年总径流量和TP总入库负荷量的贡献. 在以往研究中，由于技术手段和采样条件的限制，关于磷营养盐时空变化特征分析往往基于逐月观测或短期高频观测[30-31]. 然而，暴雨事件具有很大的随机性，不同特征的暴雨形成的径流过程不同. 监测频率低、监测时间不连续且易受天气条件制约等缺点，使得传统监测很难精准捕捉不同强度、频率和发生时间条件下水体磷的快速迁移和状态变化过程[32]. 该文采用的以人工观测数据作为校正基础，以高频监测网络与水文数值模型高效融合为主要计算平台的新方法，有效提升了千岛湖TP外源负荷量的估算精度，加深了千岛湖磷污染过程的了解，为水库安全保障与富营养化治理提供支撑.

千岛湖是河流建坝形成的人工储水体. 大收大支、水位落差大、沿河流故道至大坝出水口存在明显单向流是水库不同于大多数湖泊的典型水文特征. 水库水体中的磷营养盐主要来自外源输入，而降雨径流携带的面源污染入湖是千岛湖这种流域人类活动相对较少、植被覆盖度高的流域最主要的磷补给方式. 暴雨冲刷会造成大量磷随径流进入水库. 对于某一特定流域，外源磷营养盐的汇入量主要取决于降雨强度及降雨发生频率[33-34]. 暴雨强度决定着淋洗和冲刷地表磷污染物能量的大小，暴雨频率和数量决定着稀释磷污染物的程度，直接影响水库外源磷通量和赋存量.如Zhang等[35]联合运用大气环流模式与流域水文模型预测未来气候变化对我国石头口门水库入库径流量与面源磷负荷的影响，结果表明，夏季降雨量增大引起的径流量增加对水库TP浓度的贡献较大；张倚明等[27]基于野外观测数据，采用数值模拟方法得到不同强度降雨条件下千岛湖街口断面外源磷负荷量，其中暴雨TP负荷量占断面全年总负荷量的33%.

暴雨通过改变径流影响水库的物理过程，进而影响水体磷浓度的时空变化特征. 对于大型深水水库，除对外源磷输入的影响外，暴雨强度、频率及暴雨发生时间变化对水库底部缺氧状态持续时间和缺氧区面积影响较大，进而对水库内源磷释放通量产生较大影响. 水库分层期间，暴雨后浑水异重流的潜入会提高水库底部水温，削弱水体热分层结构稳定性，导致水库水体提前混合，蓄积在水库底部的磷营养盐被输送到上部水体，从而为翌年藻类的大量生长提供了充足的营养盐. 因此，未来研究应面向国家水安全保障重大战略需求，结合三维水动力水质模型，动态追踪暴雨期间水库水体磷赋存的时空变化，深入理解暴雨径流对贫中营养型水库富营养化的作用过程和作用机理，为我国水库有害藻类水华的形成及防控提供有力的科技支撑.

4 结论

a) 新安江模型在率定期和验证期 NSE均在0.90以上，| WBE|均 在5%以下( WBE为水量平衡误差)，说明模拟结果可靠，能够较好地反映千岛湖主要入库河流新安江的径流变化特征.

b) 通过典型年对新安江街口断面逐日流量和TP负荷量的模拟计算，发现TP全年总负荷量为1 506 t，是采用逐月监测结果计算所得TP年负荷量的2.9倍. 采用高频监测网络与水文数值模型高效融合的方法能够有效提升水库TP外源负荷量的计算精度.

c) 千岛湖流域暴雨雨量和频次均呈显著上升趋势，暴雨过程对千岛湖水环境的影响增强. 如典型年暴雨总量占年总降雨量的35.5%，暴雨雨量和频次均明显高于近60年平均值. 新安江街口断面暴雨产生的径流量及其携带TP负荷量分别占该断面年总入库径流量和TP总入库负荷量的47.9%和69.4%.

d) 暴雨径流对千岛湖水库外源TP负荷产生较大影响，给水库富营养化管理带来启示，如加强对暴雨后形成TP浓度峰值河段的农业面源和生活污水治理的投入，控制临湖面城镇与农业土地开发强度，建设面源污染拦截工程等.