暴雨徑流對新安江入庫總磷負荷量的影響

李慧赟，王裕成，單 亮，羅瀲蔥，朱廣偉,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濃度峰值河段的農業面源和生活污水治理的投入，控制臨湖面城鎮與農業土地開發強度，建設面源污染攔截工程等.