松濤水庫水溫數值模擬研究

方 晴，馬劍波，沈文龍，熊 絲，陸寶宏

(1.河海大學水文水資源學院，江蘇 南京 210098；2.江蘇省秦淮河水利工程管理處，江蘇 南京 210022；3.江西省水利規劃設計研究院，江西 南昌 330029)

流域內水庫的建設會形成流動相對緩慢的大體積水域，使天然河道的水溫分布規律發生變化，加上與太陽輻射、入流、出流、庫底河床等的熱量傳輸，形成了特有的水庫溫度場。調節性能強、庫水交換周期長的大型水庫，沿水深方向多出現分層現象。

國外對水庫水溫的研究工作開展較早，美國在20世紀30年代便開始了水溫監測，并在60、70年代提出了數值模擬方法，使水庫水溫研究步入了新時代。我國在20世紀50年代起步，對部分水庫進行了監測和研究，70年代提出一些經驗公式，80年代后引入并改進國外數學模型，開展了大量關于水庫水溫的模擬研究[1]。但研究區域多集中在亞熱帶和溫帶地區，海南等熱帶地區的水庫水溫分布規律研究甚少。海南島屬典型的海洋性熱帶季風氣候，太陽輻射強、日照時間長，氣溫較高且年變幅小，異于我國大多數地區的氣候特征。因氣候條件是影響水溫分布的重要因素，故該地區的水溫變化規律與其他地區存在較大差異。

本文以海南省松濤水庫為例，采用垂向一維模型模擬了水庫的水溫分布，分析了該水庫的水溫結構和壩前垂向水溫分布規律，對當地有關部門制定保護措施、進行流域規劃的環境影響評價具有一定的科學指導意義。

1 研究區域概況

松濤水庫位于海南省南渡江上游段，儋州市和白沙縣境內，屬南渡江流域(見圖1)。控制流域面積1 496 km2，屬多年調節大(1)型水庫，總庫容33.45億m3，多年平均入庫水量14.9億m3，年庫水替換次數為0.45，庫面寬闊，平均水面寬2 460 m，最寬處達7 500 m，縱向尺度相對較小，庫區流速緩慢，整體形成湖泊型水庫。

圖1 松濤水庫位置示意

2 研究方法

垂向一維模型忽略溫度、密度等特性在橫向和縱向上的變化，有輸入條件簡單、計算穩定的特點，對于調節性能較強、水溫在橫縱向上變化甚微的湖泊型水庫，模擬垂向水溫效果較好[2]。松濤水庫為湖泊型水庫，水面寬闊，垂向變化明顯，符合一維的假定，故采用垂向一維水溫模型——DYRESM模擬水溫分布，探究熱帶地區湖泊型水庫水溫演變規律。

2.1 模型原理

DYRESM運用拉格朗日分層理論，將水庫簡化為一個具有多水平層且各層屬性一致的系統。模型可模擬、預測水庫不同時間尺度的熱量特征值，對環境因子和水庫屬性特征的長期變化進行敏感性檢測，已成功應用于多個湖泊及水庫的溫度、鹽度、密度預測[3- 6]。

模型以一維假定為前提，引入湖泊指數LN[7]和地球自轉的比率參數R[8]來校核一維假定的有效性。模型運行主要遵守水量平衡、熱量平衡兩大定律，考慮因水體垂向移流、水氣界面熱交換、分子和紊動擴散、入流、出流引起的熱傳遞及水體吸收的太陽輻射。

2.2 模型構建

模型需輸入地形資料、動力變量、經驗參數和狀態變量。其中，地形資料包括水庫水位-水面面積關系、入流河口處詳細情況、取水口、大壩及溢洪道高程；動力變量包括氣象條件、入流流量及水溫、出流流量等；經驗參數多為學者野外觀測、實驗室試驗或率定后確定[9]；狀態變量如溫度、鹽度、水位，用于模型初始化和驗證。

2.2.1氣象資料

本次模擬時間為2016年3月2日～2017年3月1日，時間步長取一日，故氣象數據取日均值或日累積量。氣象數據選自距水庫最近的儋州氣象臺站同期觀測資料，輻射數據選自同緯度的海口臺站同期觀測資料。

2.2.2地形、初始剖面資料

地形文件主要輸入水庫的地理位置和形態特征。

初始剖面資料需輸入一組自庫底至水面的高度及相對應的水溫關系數組。采用松濤水庫2016年3月2日實測水溫(見圖2)作為初始輸入資料，由圖2可知，水庫于3月初已出現水溫分層，溫躍層厚8.81 m，水溫遞減率為0.13 ℃/m。

圖2 模型輸入初始水溫分布

2.2.3入、出庫流量資料

入庫流量選取2016年3月2日至2017年3月1日的松濤水庫實測入庫資料(見圖3a)，其中8月17日～19日和10月13日～19日期間入庫流量突增。入流水溫為南溪河入庫斷面處福才水文站多年平均水溫。出流資料為松濤水庫實測出流過程，見圖3b。

2.2.4參數、配置文件資料

根據已有的研究成果及模型經驗設置各參數，見表1和表2。

圖3 模擬時段內日平均入、出庫流量過程

表1 模型參數推薦值

表3 模型率定采用的參數

表4 水溫垂向分布特征值統計結果

表2 配置文件參數初始值

2.3 模型率定

DYRESM是一個基于過程的參數化模型，涉及參數較多，且參數具有一定通用性，采用模型默認參數即可得到較滿意的模擬結果[9]。研究表明，該模型模擬中小型水庫水溫結構時精度較高，對于較大、深水庫，一些參數需重新率定，如垂向混合系數C、允許最大層厚、最小層厚等[10-12]。考慮到國內應用DYRESM模型模擬水溫結構的研究較少，尤其是熱帶地區，故模型默認參數不一定完全適用，本次基于實測水溫數據對參數重新率定。

率定時以模擬值與實測值的絕對誤差、均方根誤差最小為原則，所有參數在其適用范圍內反復調試，最后采用參數值見表3。

采用該組參數對2016年4月26日和6月13日水溫進行模擬，結果見圖4和表4。

圖4 松濤水庫水溫結構分布

由圖4可知，模型基本能反映松濤水庫的水溫變化趨勢和結構特點，兩組模擬值與實測水溫數據總體吻合，擬合度較好，確定性系數分別為1.00和0.97，均方根誤差分別為0.46 ℃和1.25 ℃，最大絕對誤差出現于溫躍層，分別為1.15 ℃和3.25 ℃。

由表4可知，實測水溫與模擬值相比，滯溫層的擬合度最高，溫躍層擬合度最低。滯溫層的溫差不超過0.05 ℃，且位置基本吻合；表溫層水溫模擬值偏高，厚度偏大，但根據已有研究成果[3]，該模擬差異在可接受范圍內；模擬的溫躍層中心位置上移，溫度梯度偏大，4月26日溫躍層中心上移2.25 m，溫度梯度偏大0.33 ℃/m，6月13日溫躍層中心上移0.81 m，溫度梯度偏大0.36 ℃/m。

綜上，率定后的DYRESM模型模擬效果較好，總體能揭示水庫水溫分布特征及變化趨勢，但模擬結果與實際值仍存在一定差異，且主要存在于溫躍層。經分析，產生差異的主要原因為①模擬出的溫躍層偏陡是DYRESM模型研究水庫水溫普遍存在的一個現象[11-13]，有些學者認為可能是模型結構存在一定問題。②DYRESM模型假定消光系數為恒定值，但實際卻隨時間和水深變化，其取常數時對模擬水溫有一定影響[10]。③溫躍層深度受內部波運動及晝夜溫差的影響，內部波運動會造成溫躍層的暫時性偏移。④研究表明[11]，DYRESM模擬中小型水庫水溫結構效果更佳，因為其水平溫度差異性小，易滿足一維假定，風應力、熱源輸入的空間變化對大型水庫均會造成水溫在水平面上的差異。⑤本文利用的氣象資料源于儋州氣象站，氣象監測數據尤其是風速能否代表湖面實際氣象條件存在一定爭議[12]。

2.4 模型驗證

松濤水庫僅有三日實測水溫數據，已用于模型輸入和率定，且該地區同類型水庫水溫資料稀缺，故無法移用其他水庫驗證模型，本文利用松濤水庫模擬水位過程與實際水位過程的擬合程度驗證模型，見圖5。由圖5可知，模擬結果的擬合度較高，均方根誤差僅為0.18 m，說明模型能較準確地模擬松濤水庫水量、水位變化過程，用于模擬該庫水溫變化是合理的。

圖5 松濤水庫模擬期內實測、模擬水位過程

3 水庫水溫結構分析

3.1 經驗法判別水庫水溫結構

采用以下3種經驗法初步判別松濤水庫的水溫結構。

(2)密度佛汝德數法。密度佛汝德數定義為慣性力與浮力的比值，即

(1)

式中，L為水庫縱向長度；H為平均水深；V為總庫容。松濤水庫回水長度約53 km，總庫容為33.45億m3，平均水深51.3 m，多年平均入庫流量51.4 m3/s，計算得Fr=0.005 1，可知松濤水庫水溫結構為穩定分層型。

(3)寬深比法。該方法的表達式為

R=B/H

(2)

式中，B為水庫水面平均寬度；H為水庫平均水深。松濤水庫平均寬度2 460 m，平均深度51.3 m，計算得R=47.95，可知松濤水庫的水溫結構為過渡型。

綜上，不同經驗方法得出的結論不同。相比較而言，α-β指數法和密度佛汝德數法綜合考慮了水文條件和水庫特征，而寬深比法僅考慮了水庫的水深和水面寬，故前者更優，且SL278—2002《水利水電工程水文計算規范》[14]推薦使用α-β指數法。但經驗法判定水溫結構時僅考慮了水庫特征、入庫徑流量，導致判定結果可能有所偏差，而數值模型法理論嚴密，能較準確地模擬出水庫水溫的時空變化。

3.2 水溫分布年內變化規律

模擬期內松濤水庫的水溫垂向分布變化過程見圖6、7，雙溫躍層結構如圖8所示。

圖6 松濤水庫模擬期內水溫垂向變化過程

圖7 水庫各月水溫垂向分布

圖8 水庫雙溫躍層結構分布

由圖6～8可知，松濤水庫水溫呈顯著的穩定分層結構，垂向水溫變化有明顯的季節性：6、7月平均水溫最高，水面最高溫度達33.83 ℃；2月最低，庫表與庫底水體摻混均勻，平均水溫19.20 ℃；3月～翌年1月為水庫的分層期，2月為混合期；全年內40 m以下至庫底水體水溫不隨時間改變，常年保持在19 ℃左右。

水庫水溫在年內的具體變化規律為：

(1)升溫期。2月份進入升溫初期，水溫為全年最低，上下層水體由于摻混水溫呈均勻分布，垂向溫差僅0.40 ℃；3月份隨著氣溫與太陽輻射的逐漸上升，庫表開始從外界吸收熱量，表層水體水溫上升，出現厚度較小(約3 m)、深度較淺、溫度梯度僅為0.55 ℃/m的溫躍層；4、5月進入升溫末期，水體攝入的熱能大幅上升，且水庫已吸收儲存了大量熱能，分層強度持續加劇，表層溫度升高，同溫層厚度增大，溫躍層下移，溫度梯度達1.7 ℃/m；6、7月氣溫與太陽輻射達到頂峰，此時溫度結構為全年最穩定狀態，表層平均水溫達33 ℃，溫度梯度為1.6～1.88 ℃/m。

(2)降溫期。8、9月進入降溫初期，水庫從外界吸收的熱量開始減少，但氣溫、太陽輻射仍處于較高水平，水溫穩定分層，僅表層水溫較6、7月略低；8月18日出現雙溫躍層，19日中部均溫層厚度變大，該雙溫躍層穩定存在至9月20日，此后上層溫躍層逐漸減弱，直至10月初完全消失，水庫恢復單溫躍層結構。10月14日水庫再次出現雙溫躍層，15日中層均溫層加厚，至19日達到最厚，此時水庫不斷向外界釋放熱量，表層水溫逐漸下降，加之中部同溫層因入流減少、湍流擴散減弱而逐漸變薄，上部溫躍層結構受到擾動，11月5日上部溫躍層消失，水庫恢復單溫躍層結構。隨后氣溫與太陽輻射大幅下降，表層水體冷卻加快、密度增加，上層冷水與下層暖水不斷摻混，表溫層厚度增大，逐漸破壞30 m深的溫躍層結構，溫度梯度減小，溫躍層下移，至1月底該溫躍層完全消失，水庫恢復上下層水體同溫狀態。由于該溫躍層深度大(水面以下30 m)，而水庫水溫與外界氣溫溫差逐漸減小，水體向外界散熱速率減慢，故該溫躍層存在時間長。

圖9為模擬期內水庫表層水溫、入流水溫和氣溫的變化過程，可知全年內三者關系為氣溫<入流水溫<表層水溫，不符合水溫滯后于氣溫的規律。可見熱帶地區的水庫水體長時間處于強太陽輻射和高氣溫的環境下，吸熱多，水溫常年維持在一個較高水平，同時也反映了除氣溫外，其他條件如太陽輻射、風速、水汽壓、入流水溫等也可能會對水溫產生影響。

圖9 水庫表層水溫、入流水溫與氣溫變化過程

3.3 雙溫躍層結構原因分析

松濤水庫分別在8月18日和10月14日出現雙溫躍層，原因分析如下：兩次雙溫躍層的出現均與前一天入庫流量(見圖3a)猛增有關。模擬期水庫的年平均入庫流量為68.62 m3/s ，8月17日入庫流量1 950 m3/s，入庫水溫小于庫表水溫，入流水體沿水庫向下潛入與其密度相同的水層，該部分水體紊動擴散能力大，摻混作用強，形成一定厚度的中層同溫層，與頂層同溫層、下層同溫層形成雙溫躍層結構。8月18日入流流量再次加大，達2 860 m3/s，19日中層同溫層厚度進一步增大；同樣10月13日、14日、18日入庫流量分別為1 280、1 000、4 390 m3/s，造成了10月～11月雙溫躍層的出現。

兩次雙溫躍層持續的時間不同，前一次持續時間長(約40 d)，后一次持續時間短(約20 d)，其原因可能是8月～10月氣溫、太陽輻射還處于下降初期，水庫散熱速率慢，且該時期屬流域的主汛期，入庫流量較大，中層同溫層穩定存在，只能靠表層水體緩慢冷卻來破壞雙溫躍層結構；而10月中下旬后入庫流量大幅度減少，中部同溫層無法維持強烈的湍流作用，加上水庫散熱速率增快，故水庫雙溫躍層結構持續時間縮短。

4 結 論

本文以海南省松濤水庫為例，采用DYRESM構建垂向一維水溫模型，得出以下結論：

(1)松濤水庫為典型的分層型，3月～翌年1月為水庫分層期，6、7月分層最明顯，垂向最大溫差達13.90 ℃，8月份進入降溫期，水庫水體出現翻滾，分層強度開始減弱，1月底分層現象完全消失，2月為混合期。

(2)滯溫層水體受環境影響較小，全年保持19.0 ℃左右；庫表水溫與入流水溫、氣溫的關系為庫表水溫>入流水溫>氣溫，與其他地區水溫滯后于氣溫的規律不同；模擬期內由于入庫流量突然增大，導致水庫出現兩次雙溫躍層結構。