程 海 燕，脫 友 才，許 維 忠，楊 顏 菁，尹 國 武，李 傳 運

(1.四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065； 2.重慶市鯉魚塘水庫開發(fā)有限公司，重慶 404100)

0 引 言

為恢復或保護水庫的生態(tài)健康和功能，水庫的水流狀態(tài)和熱狀態(tài)已成為水環(huán)境管理的關鍵驅動因素，而水溫是其中一個基本的生態(tài)變量[1-2]。由于熱分層的形成會直接惡化底部水質，尤其會加速底部水體中溶解氧的消耗[3]，降低水生生物棲息地質量，對穩(wěn)定分層型水庫而言，水庫水溫的分布及變化尤為值得關注。水庫水溫變化受到水體間能量傳輸、水庫表面水-氣交換和庫底部水體-河床熱交換的影響，其中水庫表面水-氣熱交換是垂向分層的主要驅動力之一[4]。因此，水庫流動特性、區(qū)域水文、運行方式及氣候的差異使得水庫間熱狀態(tài)演變規(guī)律存在差異性[5-8]。IPCC《第五次氣候變化評估報告》指出：近年來，全球大部分地區(qū)的氣候變化過程以升溫為主要特征[9]。氣候變暖勢必將對穩(wěn)定分層型水庫的熱狀態(tài)、水質以及水生生物產(chǎn)生進一步影響。

在全球氣候變暖的趨勢下，湖庫熱力特性和冰期發(fā)生改變，這些改變因地域的不同而存在差異[10-11]。例如，在氣候變化條件下，溫帶地區(qū)湖泊表層水體升溫趨勢會較寒帶地區(qū)更加顯著，使得該區(qū)域水庫熱狀態(tài)對氣候變化的敏感性更強[12-15]。同時，全球湖庫也存在一些共性特征，由于氣候變暖，湖庫表面熱源的增加導致湖庫表層水溫普遍呈增溫趨勢，進而對水體混合狀態(tài)產(chǎn)生影響[11，16]。水庫和湖泊都是大型蓄水體，水庫的結構和人為運行管理使得水庫與湖泊的物理環(huán)境存在本質上的差異，包括驅動起因、幾何形態(tài)、入流條件和水動力學特征等[17]。目前國內外針對氣候變化對湖庫熱狀態(tài)的影響研究，研究對象以湖泊居多[13-19]，對有入流和人類調控影響的大型水庫熱混合狀態(tài)的討論則相對較少。

隨著全球水資源供需矛盾的加劇，水庫供水成為緩解供水壓力的最主要途徑，在中國國民經(jīng)濟發(fā)展中發(fā)揮著重要作用[20]。水庫主要供水途徑包括農業(yè)灌溉供水、工業(yè)供水、人居生活供水以及水產(chǎn)養(yǎng)殖供水等，此類途徑對水質要求較高，供水型水庫水質安全因此至關重要。水庫內水質狀況的改變可能影響供水型水庫供水安全。當水庫氮、磷含量過高時，庫表水體氧氣過飽和，庫表浮游生物和水生生物大量繁殖，庫表水體出現(xiàn)富營養(yǎng)化現(xiàn)象；而庫底溶解氧含量較低，庫底容易形成厭氧條件，促進庫底沉積物中營養(yǎng)鹽的釋放，增加水體富營養(yǎng)化的風險，惡化水庫水質[21-23]。國內外目前關于未來氣候變化對穩(wěn)定分層型供水水庫水溫的影響研究還較少，而該類型水庫在中國分布較為廣泛，氣候變化對水庫水溫的影響可能改變水庫水溫分層特性，進而對水庫水質安全造成威脅，對供水水庫而言，這種影響尤為值得關注，同時也是本文討論的重點。

本文選取重慶市鯉魚塘水庫為研究對象，構建了考慮水動力、水溫耦合作用的立面二維水溫數(shù)學模型(該模型在國內外深水型狹長水庫中已得到廣泛運用[24-25])。在此基礎上，采用為期一年的現(xiàn)場原型水溫數(shù)據(jù)對模型進行了參數(shù)率定，進一步地，基于GFDL-ESM2M2b大氣環(huán)流模型的輸出結果開展了4種未來氣候變化情景下2021～2099年鯉魚塘水庫的水溫模擬預測，由此評估分析了氣候變暖對垂向水溫結構、下泄水溫過程、熱穩(wěn)定性、水庫混合模式，以及水庫水質的影響，以期為西南山區(qū)供水水庫未來運行、管理和水環(huán)境風險防控提供科學依據(jù)。

1 研究區(qū)域概況

鯉魚塘水庫(N31.38°，E108.38°)位于重慶市開縣境內，地處長江三峽區(qū)段小江流域的二級支流桃溪河上游，是一座以灌溉、供水為主，結合發(fā)電等綜合利用于一體的水利樞紐工程。水庫主庫源頭為麻柳河，另有兩個支庫源頭分別為紫水河和水田河，如圖1所示。水庫正常蓄水位450.00 m相應庫容0.98億m3，總庫容1.02億m3，水庫面積3.14 km2，回水長度為11.35 km。大壩最大壩高為103.8 m，取水口底板高程為396.5 m。壩址多年平均徑流量1.77億m3，年庫水替換次數(shù)為1.73次，為穩(wěn)定分層型水庫。

2 研究方法

2.1 模型構建

2.1.1模型方程

CE-QUAL-W2模型是由美國陸軍工程兵團水道實驗站研發(fā)的寬度平均立面二維水動力學和水質模型，適用于河流、河口、湖泊和水庫，尤其是對相對狹長的深水型水體的水流、水溫和水質模擬有較好適用性[26]。鯉魚塘水庫為深水型水庫，壩址處最大水深可達90 m左右，且形狀狹長，故采用 CE-QUAL-W2 模型對鯉魚塘水庫水溫進行模擬。該模型遵循質量守恒、動量守恒和能量守恒，滿足Boussinesq假定，控制方程包括連續(xù)性方程、動量方程、狀態(tài)方程、自由水面方程及熱輸運方程[26]。

連續(xù)性方程

(1)

動量方程

x方向：

(2)

z方向：

(3)

狀態(tài)方程

ρ=f(Tw，ΦTDS，ΦSS)

(4)

自由水面方程

(5)

熱運輸方程

(6)

式中：B為水體寬度，m；u為縱向流速，m/s；w為垂向流速，m/s；q為側向單位體積凈入庫流量，m3/s；η為水位，m；α為河道傾角，rad；ρ為水體密度，kg/m3；Ax、Az分別為縱向和垂向紊動渦流黏滯系數(shù)，m2/s；ux為支流流速在x方向的分量，m/s；f(Tw，ΦTDS，ΦSS)為密度函數(shù)，自變量為水溫、鹽度、懸浮物濃度，當只計算水溫而不計算污染物時，密度函數(shù)僅為水溫的函數(shù)；Bη為水面寬度，m；Φ為側向平均條件下熱量濃度，J/m3；Φ=ρCpT，Cp為水的比熱，J/(kg·℃)，T為水溫，℃；Dx、Dz分別為縱向和垂向的離散系數(shù)，m2/s；qΦ為單元控制體側向熱量出入流的速率，J/(m3·s)；SΦ為熱源項，J/(m3·s)。

2.1.2模型網(wǎng)格劃分與驗證邊界

計算區(qū)域為鯉魚塘全庫區(qū)，麻柳河主庫被劃分為73×97(縱向×垂向)個矩形單元網(wǎng)格，相應的紫水河支庫被劃分為43×97(縱向×垂向)個，水田河支庫被劃分為29×97(縱向×垂向)個，單元網(wǎng)格縱向尺寸為50～100 m，垂向尺寸為1 m。圖2為網(wǎng)格分布圖。

驗證邊界條件采用源自水庫日調度運行記錄的入、出庫流量和水庫水位，實測日均入庫水溫(見圖3)以及臨近氣象站開州氣象站氣象資料(包括日均氣溫、太陽輻射、露點、云量、風速)。為探求入流水溫與氣溫的相關關系，為模型的參數(shù)率定與驗證提供基礎數(shù)據(jù)，于2019年11月23日至2020年12月2日對庫區(qū)及下泄水溫進行了原型觀測，具體觀測信息見表1，各水溫觀測斷面位置見圖1。其中A1斷面為干流入庫(麻柳水文站)水溫連續(xù)觀測斷面，A2斷面為水庫下泄水溫連續(xù)觀測斷面，B為壩前50 m垂向水溫觀測斷面。初始時刻溫度場利用2019年11月23日測得的庫區(qū)斷面垂向水溫插值得到。研究時段內，麻柳河、紫水河、水田河干支庫入庫流量按流域集水面積分配，其分配比為8∶1∶1。3條河流來流流量相對較小，水溫與氣溫相關性較為一致，故視主庫與支庫入流水溫相同。

表1 各斷面觀測信息Tab.1 Observation information of each section

2.2 未來氣候變化模式

為了解氣候變化條件下水庫熱狀態(tài)變化的可能范圍，本文對4種典型濃度路徑(representative concentmtion pathways，RCPs)的未來情景下水庫水溫進行預測模擬。每一種RCP代表一大類當前關注的多種氣體的綜合排放情景，各典型濃度路徑特征如表2所列[9]。4種未來情景氣候數(shù)據(jù)通過對美國海洋與大氣管理局(NOAA)地球物理流體動力學實驗室(GFDL)的GFDL-ESM2M2b大氣環(huán)流模型的輸出結果進行降尺度處理得到，降尺度模型采用ISIMIP2b模型。由于本文主要考慮氣候變暖帶來的影響，將研究區(qū)域不同未來情景下(2021～2099年)的氣溫數(shù)據(jù)作為氣象邊界之一，各情景下年均氣溫的變化如圖4所示，太陽輻射、露點、云量和風速均采用多年平均值。

由于鯉魚塘入庫流量較小，入庫水溫與氣溫存在較高的相關性(R2=0.99)，各未來情景下的入流水溫分別通過其與氣溫的相關關系推求得到(見圖5)。

2.3 統(tǒng)計方法

2.3.1誤差分析

本文選用平均絕對相對誤差(MARE)以及納什效率系數(shù)(NSE)用于量化模型模擬值相對實測值的偏差和吻合度，其計算公式如下：

(7)

(8)

2.3.2趨勢性分析

利用 Mann-Kendall 趨勢檢驗法對鯉魚塘水庫各熱狀態(tài)指標進行趨勢檢驗。通過判斷Z的絕對值來確定趨勢檢驗結果，|Z|≥1.64、|Z|≥1.96、|Z|≥2.57時分別表示通過了信度90%，95%和99%水平的顯著性檢驗。Z為正值(負值)時，表示增加(減少)趨勢，其計算公式為

(9)

式中：Z為檢驗值；S為檢驗統(tǒng)計量；n為樣本數(shù)；S的計算公式詳見文獻[26]。

2.3.3水庫熱穩(wěn)定性參數(shù)

采用簡化熱穩(wěn)定性指數(shù)(SI)來量化水庫熱穩(wěn)定性的變化[18]。SI用于衡量水庫穩(wěn)定分層時混合所在水柱所需要的能量，其數(shù)值越大(小)，表明打破水庫分層所需的能量越高(低)，水庫熱穩(wěn)定性越強(弱)。其計算公式如下：

(10)

3 結果分析

3.1 模型驗證結果

模型各主要參數(shù)取值如表3所列。其中縱向渦黏系數(shù)和縱向擴散系數(shù)對庫區(qū)水動力模擬較為重要，而表面太陽輻射吸收系數(shù)、純水中太陽輻射衰減系數(shù)、河床熱交換系數(shù)、風遮蔽系數(shù)和動態(tài)光遮蔽系數(shù)等對水庫水溫特性模擬較為重要，在參考了同氣候帶水庫水溫研究的基礎上，以上參數(shù)通過實測水溫數(shù)據(jù)進行了率定取值[27]。圖6比較了2020年4月11日、5月23日、8月11日和12月2日水庫壩前斷面垂向水溫的實測值與模擬值。結果顯示：4月為升溫期，水庫表層20 m范圍內水溫開始分層；5～8月受氣溫升高和取水口取水影響，水庫存在明顯分層結構；12月氣溫下降，表層溫躍層消失，水庫由單溫躍層結構逐漸過渡為垂向同溫結構，4～12月實測表層水溫變化范圍為17.5～31.9 ℃，庫底水溫維持在12.1～13.1 ℃之間；模擬與實測結果的表、底層水溫及其變溫層變化趨勢一致，2020年4月11日、5月23日、8月11日和12月2日，實測與模擬水溫的平均絕對相對誤差(MARE)分別為1.4%，2.4%，2.9%和3.1%，納什效率系數(shù)(NSE)分別為0.95，0.97，0.99和0.82，整體吻合度較好，模型可信度較高。

表3 模型主要參數(shù)取值Tab.3 Values of main parameters of the model

圖7對比了觀測時段內模擬和實測的供水取水口下泄水溫過程。對比結果顯示：模擬與實測的下泄水溫變化趨勢一致，且兩者水溫變化速率吻合較好。平均絕對相對誤差(MARE)為2.0%，納什效率系數(shù)(NSE)為0.98。

以上驗證結果表明，在參數(shù)取值合理的情況下，本文構建的模型能夠合理地模擬出鯉魚塘水庫在當?shù)貧夂驐l件下水溫的變化規(guī)律。

3.2 氣候變化對水庫水溫及下泄水溫的影響

圖8為各未來情景下預測指標趨勢性檢驗結果。隨著溫室氣體排放濃度的增加，氣溫升溫趨勢越顯著，RCP 2.6～RCP 8.5情境下的氣溫增溫速率從0.04 ℃/10 a(|Z|=1.05)增長到0.40 ℃/10 a(|Z|=8.53)，其中RCP 2.6情景下氣溫變化趨勢不顯著。通過對各情景下年均氣溫和年均表層水溫相關性進行分析，發(fā)現(xiàn)庫表水溫與氣溫相關性較高(見圖9)，各情景下相關性判定系數(shù)分別為0.87，0.86，0.92以及0.94，可見水庫表層水溫對氣候變暖的響應敏感。鯉魚塘水庫2021～2099年除了RCP 2.6情景外，RCP 4.5～RCP 8.5情景下年均表層水溫均有顯著上升趨勢，均通過了99%置信度的顯著性檢驗，且各情景下庫表水溫升溫趨勢與氣溫升溫趨勢較為一致(見圖10)，|Z|分別為3.82，7.53和8.61，升溫速率分別為0.11，0.26 ℃/10 a和0.37 ℃/10 a。各情景下表層水溫年均波動幅度分別為19.9，19.7，20.2 ℃和20.5 ℃，而氣溫年均波動幅度分別為28.9，28.7，29.5 ℃和30.1 ℃。這是因為隨著氣溫的升溫，庫表水-氣熱交換使得水庫表層熱通量增大，庫表水溫也隨之呈增溫趨勢。同時，由于水的比熱容較空氣大，各年表層水溫年內波動幅度小于氣溫波動幅度，在氣溫增溫顯著的情景中水溫的增溫率略小于氣溫，且增溫顯著性越高，二者增溫率差異越大。

如圖8所示，各情景下庫底水溫變化趨勢與氣溫和表層水溫變化趨勢相似，RCP 4.5～RCP 8.5情景下庫底水溫呈顯著上升趨勢，變化率由0.12 ℃/10 a升為0.29 ℃/10 a。由于庫底水溫主要受來流和中上層水體熱傳導影響，而水溫垂向熱傳導在水溫分層期還會受變溫層的抑制作用，因此各情景下庫底水溫升溫顯著性較庫表水溫升溫趨勢顯著性低，庫底水溫對氣溫變化的敏感性較庫表水溫對氣溫變化的敏感性略弱。表層與底層水溫的升溫趨勢表明，在氣溫顯著升溫情景下，庫區(qū)整體水溫均存在顯著增溫趨勢。

下泄水溫對水庫下游生態(tài)至關重要，通過對下泄水體年均水溫的趨勢檢驗分析(見圖8)，發(fā)現(xiàn)溫室氣體排放濃度由低到高的各情景年均下泄水溫增溫率分別為0.05，0.10，0.25 ℃/10 a和0.34 ℃/10 a，變化趨勢仍與氣溫變化趨勢較為一致。圖11為各情景下年均入流水溫與下泄水溫相關關系。RCP 2.6～RCP 8.5情景下相關性判定系數(shù)分別為0.87，0.86，0.81和0.90，可見下泄水溫與入流水溫相關性較高。隨著入流和氣溫的增溫，水庫整體熱源的增加，下泄水溫隨之增加。在氣候變暖的趨勢下，隨著表層和來流水體的增溫，庫區(qū)整體水溫存在增溫趨勢，下泄水溫對氣候變暖有所響應，由于水庫下泄水體為水庫中下層水體，其對氣溫變化的敏感性略高于庫底水溫對氣溫的敏感性。

3.3 氣候變化對水庫熱穩(wěn)定性及混合的影響

不同情景下的年均和主要分層期(5～11月)平均SI值變化趨勢如圖12所示。RCP 2.6與RCP 4.5情景下的主要分層期SI值變化趨勢不顯著，而在RCP 6.0和RCP 8.5情景下，主要分層期平均SI值均呈現(xiàn)顯著上升的趨勢(|Z|分別為4.96和5.61)，增長速率分別為31.2 kg/(m2·10 a)和43.9 kg/(m2·10 a)。在未來氣候情景下，由于水庫表層水體的顯著增溫，引起水庫表底溫差的擴大，分層期斜溫層對熱量向底部水層傳導的抑制作用增大，水庫熱穩(wěn)定性明顯增強。水庫熱穩(wěn)定性的增大意味著水庫混合阻力的增大，不利于水庫水體的置換。

水庫的水溫分層周期的長度與水庫水體的混合密切相關，關注氣候變暖條件下分層周期長度的變化有助于了解氣候變暖對水庫混合的影響。參考Tahoe湖的水溫分層時期判斷方法，本文認為水庫分層時期從連續(xù)3 dSI超過600 kg/m2的第一天開始，到SI連續(xù)3 d低于600 kg/m2的第一天結束，而一年中分層周期的長度是該年分層結束的儒略日期與該年分層開始的儒略日期之間的差值[18]。隨著氣候的變暖，中高濃度情景(RCP 4.0～RCP 8.5)下水庫水溫開始分層的時間均存在顯著提前的趨勢(見圖8)，分別提前0.52，0.84 d/10 a和1.01 d/10 a。同時，RCP 6.0和RCP 8.5情景下的水庫分層結束時間呈顯著延遲趨勢，均通過99%置信度檢驗，且延遲速率均大于1.0 d/10 a，RCP 4.5情景下分層結束時間受氣溫的波動影響而不存在顯著變化。從而，RCP 6.0和RCP 8.5情景下水庫水溫分層周期相應地延長，尤其在RCP 8.5情景，其分層周期延長速率達到2.02 d/10 a(見圖8)。圖13對各情景下各年的分層周期長度進行了統(tǒng)計分析，發(fā)現(xiàn)除RCP 4.5以外，各情景下分層周期的中位數(shù)基本大于均值，說明分層周期偏大的年份多于分層周期偏小的年份，這是由于分層時間的顯著延長造成的；而RCP 4.5情景下較多的異常值是由氣溫的波動對分層結束時間的影響所造成的。高濃度情景下，氣溫的變暖趨勢使得水庫表層出現(xiàn)較高溫度水體的時間提前，此時表底溫差較大，SI超過判斷閾值，分層開始時間因此呈顯著提前趨勢。由于水庫熱穩(wěn)定性的增加，水庫混合阻力增大，冬季水庫上層水體冷卻發(fā)生垂向對流混合達到垂向同溫狀態(tài)所需的時間更長，水庫分層結束時間因此呈延遲趨勢，水庫分層周期長度因此呈延長趨勢。

3.4 氣候變化對水庫水質的影響

對多數(shù)水庫而言，水庫水質存在藻類暴發(fā)的風險，而藻類的生長受營養(yǎng)鹽、水動力、光照、水溫、溶解氧等環(huán)境條件的影響。

其中，庫區(qū)水溫的整體升高會改變水中某些關鍵化學過程的反應速度[28]，且有益于雜食性魚類快速繁殖，取代浮游動物成為優(yōu)勢種群，打破水生生態(tài)系統(tǒng)平衡，間接增加了藻類過度繁殖暴發(fā)的風險[29]。較高溫度下，水華藻類——藍藻容易形成優(yōu)勢群種大量暴發(fā)[30]，水庫庫表較其他庫區(qū)水體水溫高，且升溫趨勢最為顯著，可能增加水庫水華暴發(fā)的風險，因此庫表水溫的升高對水庫供水是不利的。為探求氣候變化對水華暴發(fā)風險的影響，參考與鯉魚塘水庫同屬亞熱帶季風氣候的太湖(其水體溫度為24～30 ℃，是藍藻暴發(fā)的必要條件[31])，選取24 ℃為水華可能暴發(fā)的臨界水溫。各情景下庫表水溫每年初次到達24 ℃的時間變化趨勢性檢驗結果如圖8所示。庫表水溫每年首次達到24 ℃的時間在RCP 2.6、RCP 6.0和RCP 8.5情景下呈現(xiàn)顯著提前的趨勢(提前時間變化率為0.31～2.76 d/10 a)。氣溫顯著升高引起的庫表水溫顯著升高導致了這種時間顯著提前的趨勢，對于水庫管理而言，這種趨勢意味著每年可能需更早地關注水庫水華暴發(fā)的防治，水庫水華暴發(fā)風險有一定升高。

其次，水庫水溫分層導致溶解氧分層，水庫熱穩(wěn)定的增加將進一步加劇溶解氧的分層強度。同時，水溫分層周期的延長使得庫底低溶解氧水體置換周期延長，不利于水庫水體的垂向交換。以上一系列變化，可能對水庫其他水質因子產(chǎn)生影響，進而影響到水生態(tài)系統(tǒng)的垂向結構和演變趨勢，將導致庫底水質進一步變差，對鯉魚塘水庫這類供水型水庫而言值得關注。

4 結 論

(1) 氣候變暖使得水庫水溫和下泄水溫顯著升高。各氣候變化情景下，年均庫表水溫、庫底水溫和下泄水溫的增溫趨勢與氣溫變化存在不同程度的響應關系。隨著溫室氣體排放濃度的增加，氣溫、表層水溫、庫底水溫和下泄水溫升溫顯著性越大。下泄水溫與入流水溫相關性較高。庫表水溫和下泄水溫對氣溫變化的敏感性略高于庫底水溫對氣溫變化的敏感性。

(2) 氣候變暖使得水庫熱穩(wěn)定性顯著提高，水庫水體垂向混合能力變弱。在氣候變暖趨勢下，高濃度排放情景(RCP 6.0和RCP 8.5)下的水溫主要分層期的SI值呈顯著上升趨勢，該趨勢意味著水庫混合阻力的增大，不利于水庫水體的置換。同時，水庫水溫分層開始時間呈現(xiàn)顯著提前趨勢，而分層結束時間呈現(xiàn)顯著延遲趨勢，年內水溫分層周期長度相應地呈顯著延長趨勢。分層周期的延長使得變溫層對水庫水體垂向混合的阻隔影響增大。

(3) 氣候變暖后水庫熱狀態(tài)的變化對水庫水質可能產(chǎn)生進一步不利影響。氣候變化引起水庫水溫的升高，進而可能打破水生生態(tài)平衡，特別是庫表水溫的顯著升高增加了藻類暴發(fā)的風險。水庫水溫分層穩(wěn)定性的增加，加劇了水庫溶解氧的分層，阻礙了水庫垂向水體的混合和置換，也將成為庫底水質進一步變差的原因之一。

本次研究中僅考慮了氣溫因子的變化，實際上，水庫表面水-氣熱交換通量的變化，還與太陽輻射、風速、相對濕度等氣象因子相關，且氣候變化中這些因子也在同步發(fā)生變化。同時，氣候變化也將引起流域徑流過程和水庫運行方式的改變。在下一步工作中，還需更全面考慮這些因素的綜合影響作用，以便更科學、準確量化氣候變化對水庫熱狀態(tài)的影響閾值。