水庫水溫模擬預(yù)測中常見問題探討

李克鋒，梁瑞峰，脫友才，鄧云，李嘉

(四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室，四川成都　610065)

?

水庫水溫模擬預(yù)測中常見問題探討

李克鋒，梁瑞峰，脫友才，鄧云，李嘉

(四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點實驗室，四川成都610065)

針對影響水庫水溫模擬效果的常見問題，提出以水溫-氣溫相關(guān)法確定缺乏資料地區(qū)的現(xiàn)狀水溫，對水庫地形校正效果開展定量分析，討論了氣溫直減率對水氣界面長波輻射、蒸發(fā)、傳導(dǎo)的影響程度，對比了部分氣象經(jīng)驗公式的熱通量計算差異，以二維模型為例定性分析了零方程與雙方程紊流模式、靜水壓力假定與完整垂向動量方程的差異，提出了可通過熱量收支平衡、逐月徑流量與庫容關(guān)系等方法定量檢驗?zāi)M結(jié)果的合理性。

水庫水溫；地形；氣象；數(shù)學(xué)模型

水電資源作為一種清潔能源一直是國際社會能源行業(yè)的重點發(fā)展領(lǐng)域，但水電站的興建也極大地改變了河道水溫的自然變化過程，特別是在水能資源豐富的河流，巨型電站或者多個梯級水庫聯(lián)合運行時，對河道水溫的影響范圍更寬，對區(qū)域性的環(huán)境影響更大。水溫影響逐漸成為水電開發(fā)過程中環(huán)境影響評價的重要內(nèi)容[1-3]。

目前，國內(nèi)外水庫水溫影響研究一般從數(shù)值模擬、室內(nèi)試驗、原位觀測3個方面開展，其中數(shù)值模擬方法具有物理機(jī)制清晰、考慮因素全面、可有效結(jié)合理論分析與實驗研究成果等優(yōu)勢，越來越多地應(yīng)用于新建項目的環(huán)境影響評價中，其模擬精度直接關(guān)系到對工程環(huán)境影響的正確認(rèn)識。但水溫的模擬由于受模型選擇、參數(shù)取值、離散方法、邊界條件、初始條件等多種因素的影響，在實際應(yīng)用中如果處理不當(dāng)可能影響模擬效果。本文擬對水庫水溫模擬中常見的主要問題展開分析和討論，以期為相關(guān)研究提供參考。

1　現(xiàn)狀水溫的確定

現(xiàn)狀水溫作為數(shù)值模擬的輸入條件和水溫影響程度的對比基礎(chǔ)，是模擬的重要基礎(chǔ)資料。在平原河流或者具有布設(shè)水文儀器條件的河段，水溫作為常規(guī)水文監(jiān)測項目，其數(shù)據(jù)獲取尚不至于制約模擬預(yù)測。隨著水電開發(fā)逐漸向高海拔寒冷地區(qū)拓展，研究河段缺乏水溫觀測資料，難以確定模擬所需的現(xiàn)狀水溫，但通過氣象相關(guān)、鄰近流域監(jiān)測數(shù)據(jù)類比等方法，仍可得到合理的基礎(chǔ)水溫數(shù)據(jù)[4]。

氣象條件是影響天然水溫變化的重要因素。河段沿程有氣象監(jiān)測數(shù)據(jù)且周邊氣候相近的流域存在水溫監(jiān)測數(shù)據(jù)時，可采用水溫-氣溫相關(guān)方法確定該區(qū)域的現(xiàn)狀水溫。以西藏自治區(qū)境內(nèi)河流為例，怒江上游僅嘉玉橋、龍普、扎那水文站曾短期監(jiān)測水溫，現(xiàn)有數(shù)據(jù)難以滿足上游規(guī)劃梯級電站現(xiàn)狀水溫確定的需要，但可根據(jù)西藏境內(nèi)其他流域的氣溫、水溫關(guān)系進(jìn)行推算。

水溫-氣溫相關(guān)關(guān)系可表示為：

T水=a×T氣+b

(1)

式中，a為水溫隨氣溫的變化率；b為某處水溫受上游來流水溫的影響程度，℃。

西藏自治區(qū)雅魯藏布江、年楚河、拉薩河、尼洋河、怒江等干支流的水溫-氣溫相關(guān)系數(shù)見表1。從相關(guān)性來看，各水文站水溫均與氣溫呈高度相關(guān)關(guān)系(相關(guān)系數(shù)>0.8)；水溫隨氣溫的變化率平均為0.757 9；各干支流相關(guān)關(guān)系中的參數(shù)b均為下游大于上游，說明各河流下游水溫受上游影響程度存在差異性。因此，可采用式(1)計算怒江干支流的天然水溫。

此外，缺乏水溫觀測資料地區(qū)的現(xiàn)狀水溫的確定，也可直接移用鄰近流域流量相近站點的監(jiān)測資料，對一些監(jiān)測年限較短的資料，還需要與氣象數(shù)據(jù)建立相關(guān)關(guān)系后進(jìn)一步插補延長。

2　地形的概化與校正

由地形概化得到的網(wǎng)格是水庫庫容、面積等幾何尺寸在數(shù)學(xué)模型中的直接體現(xiàn)，因而地形精度一定程度上決定了預(yù)測精度。從大斷面、地形圖等概化后的水庫庫容曲線受庫區(qū)形狀不規(guī)則、測量精度、提取方法等影響，與實際庫容、面積曲線存在差異，因而在計算中應(yīng)對概化庫容、面積曲線做進(jìn)一步校正，消除地形誤差帶來的水庫調(diào)節(jié)能力變化。

浯溪口水利樞紐工程是江西省景德鎮(zhèn)市昌江干流中游一座以防洪為主，兼顧供水、發(fā)電等的綜合利用工程。該工程正常蓄水位為56 m，水庫總庫容為4.274×108m3，具有不完全年調(diào)節(jié)性能。圖1和圖2顯示了浯溪口水庫地形校正對水庫水溫的影響。浯溪口庫區(qū)地形系根據(jù)分辨率較低的地形圖得到，經(jīng)與設(shè)計院提供的高程-面積-容積曲線對比，從地形圖概化得到的正常蓄水位水面面積為17.66 km2，與實際水面面積誤差接近30%。對此地形不做校正時，采用立面二維水溫模型計算得到的平水年水庫下泄水溫最大升幅達(dá)6.0℃；經(jīng)對地形開展復(fù)核，根據(jù)設(shè)計方提供的復(fù)核后地形校正模型計算采用的庫容面積曲線(圖1)，正常蓄水位水面面積校正為24.70 km2，下泄水溫的最大升幅減小至3.3℃(圖2)。復(fù)核前后水面面積差異帶來水面熱通量計算存在較大差異，特別是水庫蓄水后期的高水位運行，面積差異更趨放大，顯著改變了水庫的冬季蓄熱能力，導(dǎo)致不合理的下泄水溫大幅升高現(xiàn)象。

圖1　地形校正前后浯溪口水庫水面面積Fig.1　Surface area of Wuxikou Reservoir before and after the topographical correction

圖2　地形修正前后浯溪口水庫下泄水溫Fig.2　Discharge water temperature of Wuxikou Reservoir before and after the topographical correction

3　氣象資料與計算公式的選擇

影響水庫水溫模擬的氣象要素一般包括太陽輻射、氣溫、云量、風(fēng)速風(fēng)向、相對濕度等。

太陽輻射和氣溫是影響水體熱量收支的最重要因素。河谷深切的山區(qū)水庫應(yīng)充分考慮兩側(cè)山體遮蔽對輻射到達(dá)水面時段的影響，特別是南北走向的河流。因地形遮蔽作用，水面接收的太陽直射輻射最大可減少25%左右[5]。

氣象站點的選擇也間接影響到水面熱通量的計算。金沙江下游的溪洛渡庫區(qū)有金陽、永善等氣象站，金陽站位于庫區(qū)中段，在地理位置上更有代表性，但金陽氣象站海拔1 450 m，與溪洛渡正常蓄水位600 m相差太大，而氣溫具有隨海拔上升而降低的特性，直接采用金陽站氣溫計算水氣界面熱通量將會帶來較大誤差。

圖3為金陽氣象站氣溫修正前后對溪洛渡庫區(qū)河段水氣傳導(dǎo)、蒸發(fā)、長波輻射的影響程度。在采用該地區(qū)山區(qū)氣溫直減率5℃/km[6]將1 450 m高程的氣溫修正至600 m高程后，該區(qū)域年均氣溫比下游永善站高出3.2℃，河段水氣界面熱通量年均增加達(dá)73.7 W/m2，也從側(cè)面說明氣象因素是金沙江溪洛渡段增溫率(0.43℃/100 km)高于向家壩段增溫率(0.07℃/100 km)的原因之一。

圖3　金沙江下游溪洛渡庫區(qū)河段氣溫修正對水氣界面熱量的影響Fig.3　Impact of air temperature correction on heat flux in Xiluodu Reservoir

此外，水面熱通量經(jīng)驗公式的選取也對計算存在較大影響，本文以大氣長波輻射為例進(jìn)行說明。大氣中的水珠、蒸汽、塵埃等可以吸收來自太陽的短波輻射等，使其自身增溫，然后又通過長波輻射形式輻射到水面，其強(qiáng)度取決于氣溫的高低。大氣長波輻射絕大部分可被水體吸收，約有3%被水面反射到大氣中。大氣長波輻射波長范圍為4～120 μm，輻射強(qiáng)度一般取決于氣溫和云量。

大氣長波輻射計算公式為：

φan=(1-γa)·σ·εa(273+Ta)4

(2)

式中，Ta為水面上2 m處的氣溫；γa為長波反射率，取0.03；σ為Stefan-Boltzman常數(shù)，5.67×10-8(W/m2·K4)；εa為大氣發(fā)射率。

大氣發(fā)射率的常用計算公式為：

(3)

(4)

(5)

式中，C為云層覆蓋率；Kε與云層高度有關(guān)，均值為0.17；ea為水面上空氣蒸發(fā)壓力；Aa為經(jīng)驗系數(shù)，取值范圍為0.5～0.7，常取值0.6；Ab為經(jīng)驗系數(shù)，取值范圍為0.031～0.076 mmHg-0.5，常取值0.031；eair為空氣中的水蒸氣分壓力，mmHg。

空氣中的水蒸氣分壓力可用下式表示：

(6)

式中，Td為露點溫度。

仍以浯溪口水庫庫區(qū)氣象條件來說明公式(3)、(4)、(5)的差異(圖4)。從圖4中大氣長波輻射在一年內(nèi)的變化可發(fā)現(xiàn)，采用不同的大氣發(fā)射率經(jīng)驗公式，冬季的差異小于夏季，但最小也相差36.7 W/m2；夏季差異較大，6月差異高達(dá)97.3 W/m2。

短波輻射、水體長波輻射、蒸發(fā)通量、傳導(dǎo)通量在選擇不同的計算公式時，也存在與大氣長波輻射計算類似的問題，在預(yù)測中宜根據(jù)研究區(qū)域的氣候特征對氣象公式作適用性驗證。在商業(yè)軟件中，由于熱通量計算模式相對固定，尤其應(yīng)注意公式的適用性。

圖4　浯溪口庫區(qū)采用不同大氣發(fā)射率公式的水面長波輻射通量Fig.4　Water surface longwave radiation flux under different air emission rate formulae in Wuxikou Reservoir

4　垂向動量交換的模擬

應(yīng)用于水庫水溫預(yù)測的數(shù)學(xué)模型主要從空間維度上分為一維、二維、三維?？v向一維模型可用于水深較淺的水庫，垂向一維模型用于深且短的水庫；二維模型主要有深度平均的平面模型和寬度平均的立面模型；三維數(shù)學(xué)模型則能充分反映水庫水溫的變化細(xì)節(jié)，特別是泄流孔口附近水流及水溫的三維特征[7]。由于三維水溫計算工作量巨大，目前主要用于近壩區(qū)電站進(jìn)口處的模擬，而一維模型適用性有限，因此主要就立面二維水溫模型作簡要對比。

水庫的流場分布與水庫的溫度分布具有強(qiáng)烈的相互作用，準(zhǔn)確模擬水庫溫度場的關(guān)鍵在于精確模擬浮力流流場，因此立面二維模型之間的差異主要是紊流模式和垂向動量方程的差異。

根據(jù)紊流模型采用的微分輸運方程的個數(shù)，可以分為零方程、單方程、雙方程3類。美國陸軍工程兵團(tuán)的CE-QUAL-W2模型推薦采用W2零方程模式計算渦粘系數(shù)：

(7)

式中，νt為垂向渦流系數(shù)；κ為范卡門常數(shù)；lm為混合長度；u為縱向流速；τwy為因風(fēng)而產(chǎn)生的橫向剪應(yīng)力；k為波數(shù)；τytributary為因支流入流而產(chǎn)生的橫向剪應(yīng)力；ρ為水體密度；C為常數(shù)；Δzmax為垂向網(wǎng)格間距的最大值；Ri為理查森數(shù)。

W2模式通過Ri考慮了垂向密度梯度帶來的紊動粘性變化，混合長則由網(wǎng)格垂向最大間距確定，這就使紊動粘性系數(shù)νt的確定依賴于網(wǎng)格劃分。W2可通過原位觀測率定相對合理的混合長，但率定難免將水溫實測中的不確定因素帶入其中，影響模型的適用性。

單方程模式增加紊動動能的精確輸運方程，其難點仍在于長度比尺的確定。雙方程模型解決了單方程模型長度比尺的確定難題，但在使用中如采用完全的不可壓縮假定，垂向溫差引起的密度流不能通過浮力項影響渦粘系數(shù)，也將導(dǎo)致浮力流模擬偏差。

立面二維模型垂向動量方程的差異主要在于是否采用完整的垂向動量方程。完整的垂向動量方程和靜水壓力假定動量方程分別為：

(8)

(9)

式中，w為垂向流速；p為壓強(qiáng)；T為水溫；ρa為參考密度；B為河寬；g為重力加速度；β為等壓膨脹系數(shù)；νe為分子粘性系數(shù)ν與紊動渦粘系數(shù)νt之和。

與完整的垂向動量方程相比，靜水壓力假定動量方程忽略了垂向流速的時變項、對流項、擴(kuò)散項，使得模型主要適用于流速的各向梯度變化不大的流動模擬，而當(dāng)模擬存在明顯上浮或者下潛現(xiàn)象的溫差異重流時，模擬結(jié)果存在偏差。

5　預(yù)測結(jié)果合理性分析

水庫水溫模擬預(yù)測的結(jié)果是否合理，需要通過對熱量平衡、延遲時間等進(jìn)行分析來確定。分析結(jié)果如表2、表3所示。

表2列出了金沙江下游烏東德水庫在平水年單獨運行時的熱量收支情況。從表2可知，烏東德水庫在6—9月入庫、出庫熱量占全年熱量收支的65%以上，水庫3—5月入庫、出庫熱量的全年占比雖小，但出庫熱量總體小于入庫熱量，反映了在此期間水庫下泄低溫水現(xiàn)象的存在。

表2　金沙江下游烏東德水庫平水年的熱量收支情況

表3　土溪口水庫典型平水年逐月徑流量與庫容的替換關(guān)系

對于調(diào)節(jié)性能較強(qiáng)的水庫，其下泄水流的低溫水降幅可通過來流到達(dá)壩前的延遲時間做輔助分析。表3為四川宣漢縣渠江支流前河土溪口水庫的逐月徑流量與庫容的替換關(guān)系，替換次數(shù)與延遲時間均采用當(dāng)月水位下的全庫容計算，即假定水庫整個斷面參與流動。表3中可見水庫來流量年內(nèi)分布極不均勻，典型平水年在3—5月來流量僅占全年徑流量的25%，庫中水體難以充分替換，11月至次年4月，當(dāng)月來流均不能完全替換當(dāng)月庫內(nèi)水體；水庫調(diào)節(jié)引起的低溫來流延遲出庫是較大低溫水降幅的主要成因。

6　結(jié)語

水溫模擬受模型選擇、參數(shù)取值、離散方法、邊界條件、初始條件等多種因素的影響，在實際應(yīng)用中如處理不當(dāng)，可能影響模擬效果。本文對水庫水溫模擬預(yù)測中常見的主要問題開展了討論，對直接影響模擬過程的現(xiàn)狀水溫的確定、地形的概化與校正、氣象資料的處理與計算公式的選擇、垂向動量交換的模擬進(jìn)行了分析，并提出了具體方法與處理思路，對采用熱量收支平衡、逐月徑流量與庫容關(guān)系等方法分析模擬結(jié)果的合理性進(jìn)行了闡述。

由于數(shù)值模擬包括前后處理、離散方法與過程、計算控制、驗證分析等多個環(huán)節(jié)，文中所涉僅為少數(shù)常見問題，實際處理中需要在更多方面加以關(guān)注，如垂向一維計算中綜合擴(kuò)散系數(shù)的確定、河流水溫計算中熱源項處理、二維計算中水庫泄流孔口喇叭口的概化、泄流孔口及周邊最小網(wǎng)格尺寸的確定、多維嵌套中的界面?zhèn)鬟f、分層取水過程中的動邊界處理、流場與溫度場的匹配驗證等。預(yù)測結(jié)果的合理性分析在水庫水溫預(yù)測中也至關(guān)重要。環(huán)評中可結(jié)合水庫入流和出流的熱量平衡分析、分時段的庫水交換次數(shù)分析等方法，總體上判斷預(yù)測結(jié)果的合理性。此外，還可將預(yù)測結(jié)果與具有可類比性的已建水庫水溫觀測資料進(jìn)行對比，來分析預(yù)測結(jié)果的合理性。

[1]宋策, 周孝德, 辛向文. 龍羊峽水庫水溫結(jié)構(gòu)演變及其對下游河道水溫影響[J]. 水科學(xué)進(jìn)展, 2011, 22(3): 421-428.

[2]李鐘順, 陳永燦, 劉昭偉, 等. 密云水庫水溫分布特征[J]. 清華大學(xué)學(xué)報(自然科學(xué)版), 2012, 52(6): 798-803.

[3]脫友才, 劉志國, 鄧云, 等. 豐滿水庫水溫的原型觀測及分析[J]. 水科學(xué)進(jìn)展, 2014, 25(5): 739-746.

[4]李克鋒, 郝紅升, 莊春義, 等. 利用氣象因子估算天然河道水溫的新公式[J]. 四川大學(xué)學(xué)報(工程科學(xué)版), 2006, 38(1): 1-4.

[5]張艷林, 程國棟, 李新, 等. 山區(qū)太陽輻射對水熱過程影響的敏感性分析[J]. 冰川凍土, 2012, 34(3): 650-659.

[6]翁篤鳴, 孫治安. 我國山地氣溫直減率的初步研究[J]. 地理研究, 1984, 3(2): 24-34.

[7]馬方凱, 江春波. 三峽水庫近壩區(qū)三維流場及溫度場的數(shù)值模擬[J]. 水利水電科技進(jìn)展, 2007, 27(3): 17-20.

Discussion on Common Problems in the Prediction of Reservoir Water Temperature Simulation

LI Ke-feng, LIANG Rui-feng, TUO You-cai, DENG Yun, LI Jia

(State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China)

Concerning the common problems affecting prediction results of reservoir water temperature, the method for the determination of natural water temperature in ungauged basins was put forward based on water-air temperature correlation and quantitative analysis for reservoir topographic correction was conducted. The impact of air temperature lapse rate on long-wave radiation, evaporation and conduction across water-air interface was discussed. The difference of heat flux calculated by some empirical meteorological formulas was compared. Qualitative analysis was conducted by comparing the zero equation turbulence model and double equation turbulence model, as well as investigating the difference between the hydrostatic symposium and complete z-momentum equation. A method for quantitative rationality verification of simulation results was put forward through the analysis on heat balance and correlation between reservoir capacity and monthly runoff.

reservoir water temperature; topography; meteorology; mathematical model

2015-11-22

國家自然科學(xué)基金(51379136；50979063；51279114)

李克鋒(1965—)，男，湖北石首人，研究員，博士，主要研究方向為生態(tài)與環(huán)境水力學(xué)，E-mail：kefengli@scu.edu.cn

10.14068/j.ceia.2016.03.015

X820.3；TV697.2

A

2095-6444(2016)03-0057-05