劉有志，相建方，陳文夫，廖建新，程 恒，黃海龍

（1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；2.中國三峽建設(shè)管理有限公司，北京 100038）

1 研究背景

超高拱壩是一種受溫度和水壓荷載影響較大的空間殼體結(jié)構(gòu)，其中上游水庫水溫分布直接影響到大壩運行期穩(wěn)定（準穩(wěn)定）溫度場的分布[1-2]。特別是對于壩體基礎(chǔ)約束區(qū)，上游庫底水溫將直接影響到大壩的基礎(chǔ)溫差和溫度控制設(shè)計標準的制定，而上游水庫水溫的分布及變化過程也會對拱壩全生命周期的運行工作性態(tài)產(chǎn)生影響，因而研究水庫水溫的時空分布和演化規(guī)律對拱壩安全意義重大。

蘇聯(lián)和美國在1930年代即開始重視水庫水溫的研究，前蘇聯(lián)在現(xiàn)場試驗研究中做了大量細致的工作，美國在水溫數(shù)學(xué)模型的建立和工程實踐應(yīng)用方面的研究一直在不斷推進，其后的發(fā)展主要沿這兩個方向進行。在數(shù)值方法的研究中，Raphael J M[3]、Orlob J T 等[4]、Huber W C 等[5]和Brooks N H等[6]的研究和應(yīng)用推動了數(shù)值方法的發(fā)展，對于高壩水庫分層的現(xiàn)象有了初步的認識；Edinger 等[7-8]開發(fā)了橫向平均的LARM模型用于庫水溫和水質(zhì)的計算，經(jīng)過模型不斷的改進和功能的完善，發(fā)展成了現(xiàn)在的能夠應(yīng)用于河流、湖泊、水庫和海灣的二維橫向平均水動力學(xué)和水質(zhì)模型CE-QUAL-W2[9-10]，為解決計算過程中的數(shù)值耗散問題，該模型還多次進行了算法的改進[11]；與此同時，數(shù)值方法中還出現(xiàn)了三維計算方法，其中以集成了水動力學(xué)模塊、泥沙輸運模塊、污染物運移模塊和水質(zhì)預(yù)測模塊等多個模塊的EFDC為代表[12]。中國1980年代開始大規(guī)模興建混凝土壩時，主要通過類比同類工程的水庫水溫實測資料，來擬合得到了水庫水溫的預(yù)測經(jīng)驗公式，如中國混凝土拱壩設(shè)計規(guī)范[13]中推薦的水庫水溫計算公式，就是朱伯芳提出的經(jīng)驗公式[14]。但實踐表明，該方法主要適用于100 m 級的水庫大壩，且無法考慮大壩成庫過程中水溫變化過程的模擬對于垂向分層明顯的水庫，胡平等[15]結(jié)合工程實際開發(fā)了一維垂向計算方法，在一些水庫的計算和預(yù)測中達到了較好的效果，但大壩成庫過程中水溫變化過程的模擬過程仍無法實現(xiàn)。近年來，隨著200 m 級以上高庫大壩建設(shè)的逐漸增多，經(jīng)驗公式法已難以適用于高壩水庫水溫的模擬，數(shù)值模擬方法得到越來越多專業(yè)人員的關(guān)注和研究，如蔣紅[16]、鄧云[17]、張俊華等[18]、章若茵等[19]對溪洛渡等大型水庫的水溫或流沙異重流模擬預(yù)測方法進行了研究，這些算法以水庫建設(shè)前的預(yù)測或者室內(nèi)實驗為主，未及進行正式蓄水后的跟蹤研究，對于這些算法和理論缺乏系統(tǒng)的驗證和反饋。

現(xiàn)有的文獻資料表明，壩高200 m 以上的水庫還未進行過從開始蓄水到運行初期庫水溫溫度演化規(guī)律的專門研究。以溪洛渡水庫為例，當水庫水深較深、庫容較大時，原來的研究成果不能直接用于這類水庫，需要將模型和實際的蓄水過程、氣象條件及水文條件等多種影響因素相結(jié)合來分析。因此，狹長型高庫大壩的水庫水溫模擬問題在實踐中依然非常復(fù)雜，如何有效模擬此類水庫的真實水庫水溫分布及變化規(guī)律仍是業(yè)內(nèi)關(guān)注的重點難題。

2 水庫水溫計算的3種方法簡述

當前大型水庫水溫計算方面的研究較多，但是并沒有根據(jù)河道的形狀對水庫類型進行細分。筆者通過多座水庫的對比研究發(fā)現(xiàn)，二維算法更加適用于河道狹長，且水庫寬高比小于1的水庫，該類水庫往往水深較深，一維算法適用于水庫上游表面開闊的水庫，壩體多為重力壩結(jié)構(gòu)，而經(jīng)驗公式方法具有廣泛適用性，但是缺點是精度不夠，無法滿足精確模擬的要求，且無法考慮大壩成庫過程中水溫的變化過程。對于本文研究的水深較深的狹長型河道上的水庫而言，CE-QUAL-W2 在模擬水庫水動力學(xué)過程和密度梯度作用方面有著先天的優(yōu)勢。

2.1 CE-QUAL-W2CE-QUAL-W2是一個考慮縱向和垂向的二維水動力學(xué)與水質(zhì)模型。由于該模型假設(shè)橫向平均，因而最適宜模擬像本文這樣的狹長河道型水庫的縱向和垂向的水質(zhì)梯度。該模型之前已經(jīng)被用于河流、湖泊、水庫等水體的計算中。該模型能夠計算水體自由表面、垂向和縱向的流速，以及溫度。本文選用ULTIMATE 算法求解水動力輸運方程，該算法能夠減少不真實的物理過程，更加符合本文的計算。同時，在求解過程中，CE-QUAL-W2 還提供了一個可控制垂直對流解法時間權(quán)重的參數(shù)THETA，當取0時是完全的顯式格式，當取1時是完全的隱式格式，當取0.5時對應(yīng)的是Crank-Nicholson 解法，本文按照建議選取THETA為0.55。該模型還提供了多種湍流模型，經(jīng)過對比驗證，本文選用k-ε模型。

由于按照靜水假設(shè)進行計算，垂向的動量方程沒有被考慮在內(nèi)，該模型可能會在有明顯垂向加速度的情況下得到不準確的結(jié)果，比如在初期的快速蓄水階段等過程中，因而仍需要在實際模擬中校正。關(guān)于CE-QUAL-W2的理論研究部分，在說明書上已經(jīng)有詳細的介紹[9-10]，在此不進行贅述，在此僅列出主要控制方程。

X方向動量方程（X-momentum）：

Z方向動量方程（Z-momentum）：

連續(xù)性方程（continuity）：

狀態(tài)方程：

自由表面方程：

式中：U為水平速度，m/s；τx為x方向上的平均剪切應(yīng)力；W為垂直速度，m/s；τz為z方向上的平均剪切應(yīng)力；B為通道寬度；ρ為密度；q為壓力；η為水面高程；為考慮溫度、總?cè)芙夤腆w或鹽度、無機懸浮物的密度函數(shù)。

2.2 一維算法關(guān)于將一維算法作為對照算法，在此做簡要描述?？扇〕鲆粋€微元，研究其熱量運動，計算模型簡圖見圖1。

圖1 一維水庫水溫計算模型簡圖

（1）垂直向：單位時間內(nèi)由下面進入的流量為Qy，帶進的熱量為cρQyT ；單位時間內(nèi)由上面流出的流量為Qy+dQy，帶走的熱量為：

故單位時間內(nèi)凈帶進熱量為：

（2）水平向：單位時間內(nèi)入庫帶進熱量為cρiqiTidy；出庫帶走熱量為cρq0Tdy ；故單位時間內(nèi)凈剩熱量為：

（3）由短波輻射熱：

自下邊離去的輻射熱：

自上面進入的輻射熱：

留下的凈輻射熱為：

（4）由擴散作用：

下邊進入：

上邊流出：

凈流入為：

（5）水體升溫吸熱：

式中：c為水的比熱；ρ為水的密度；T為水的溫度；qi為入庫水流單位高度的流量；Ti為入庫水流的溫度；q0為出庫水流單位高度流量；ρi為入庫水流的密度；R（y）為高度y 處的短波輻射熱；k為輻射熱的衰減系數(shù)；A（y）為y 處的水庫面積；Dm為水分子擴散系數(shù)；E為水的紊動擴散系數(shù)。由熱量平衡可知：

則有：

2.3 經(jīng)驗公式算法朱伯芳院士方法是經(jīng)驗公式法的代表，在中國混凝土拱壩設(shè)計中被廣泛采用，本文對其進行簡要介紹[13]。

庫水溫度T（y，τ）是水深y和時間t的函數(shù)，可按下列方法計算：

任意深度的水溫變化

任意深度的年平均水溫

水溫相位差

式中：y為水深，m；τ為時間，月；ω為溫度變化的圓頻率，ω=2π/P；P為溫度變化的周期，12個月；ε為水深y 處的水溫滯后，即對于氣溫的相位差，月；T（y，τ）為水深y 處在時間為τ時的溫度，℃；Tm（y）為水深y 處的年平均水溫，℃；Ts為表面年平均水溫，℃；A（y）為水深y 處的溫度年變幅，℃；τ0為氣溫最高的時間，月；d、f、γ為經(jīng)驗擬合系數(shù)，分別取為2.15、1.30、0.085。

3 大型水庫水溫現(xiàn)場反饋與分析

3.1 模型的建立溪洛渡水電站工程屬于中國西南地區(qū)金沙江梯級開發(fā)中的第三梯級，最大壩高285.5 m，壩頂高程610 m，水庫正常蓄水位600 m，死水位540 m，汛期限制水位560 m，距上游白鶴灘水電站約200 km。大壩2013年初開始下閘蓄水，2013年6月水位蓄水至540 m，2014年10月首次蓄水至600 m，2015年6月水位回落至540 m，2015年10月水位再次到達600 m。

庫區(qū)河谷較窄，河道寬度在600 m 左右，且河道較深，為典型的U 形河道，庫區(qū)地形圖如圖2所示。根據(jù)實際河道地形信息進行網(wǎng)格的劃分，縱向分為56 段，每段長度4 km，垂向分為72層，每層深度4 m，由于計算范圍較大，網(wǎng)格數(shù)目較多，為減少計算量，本文采用的計算網(wǎng)格是滿足網(wǎng)格對計算結(jié)果無影響條件下的最大尺寸·網(wǎng)格，縱向、垂向和俯視網(wǎng)格如圖3所示，其中圖3（b）為近壩段垂向截面網(wǎng)格。

3.2 初始邊界與參數(shù)選擇本次模擬計算的輸入文件中，河道地形信息采用溪洛渡地區(qū)的實際測深數(shù)據(jù)，氣象資料采用近壩區(qū)氣象站的實測氣溫、露點溫度、云量、風(fēng)速和風(fēng)向，太陽輻射對庫水溫度的影響按照當?shù)貙嶋H經(jīng)緯度位置、實際云量和山體高程等遮蔽信息計算。來水流量采用主河道水文站測得流量，出水流量按照電站實際的調(diào)度過程給出，該水庫以發(fā)電為主要目的，出水口的位置僅考慮電站引水口，由于以上輸入信息均按照實際給出，因而計算過程中的水位與實際水庫的水位過程一致。

溪洛渡庫區(qū)底部堆渣區(qū)實測資料顯示，庫區(qū)底部水溫在15℃左右，庫底以下10 m 處的巖石溫度在21℃左右，溫度高于一般地質(zhì)條件下的溫度，水和巖石接觸面的溫差為6℃，兩者之間存在一定熱交換，庫底的熱交換必須考慮在內(nèi)，熱交換系數(shù)Kw應(yīng)取為比通常情況下的參數(shù)稍大[8-9]，可取為0.4 W·m-2·℃-1。

圖2 庫區(qū)地形圖

圖3 縱向和垂向網(wǎng)格

3.3 現(xiàn)場監(jiān)測數(shù)據(jù)對于溪洛渡這樣的不完全年調(diào)節(jié)水庫而言，來流流量、出流流量、來流溫度和氣溫、輻射等氣象條件是影響水溫分布的直接因素，其中來流流量、來流溫度是影響庫水溫的主要因素。為滿足計算的需要，也為保障數(shù)值模擬計算精度，本次反饋分析同時考慮了干流和支流上的水溫和流量資料，在原水文站中增設(shè)溫度計測量來水溫度，同時增設(shè)了兩個臨時測站來測量支流的流量和溫度情況，這些測站均于2012年布置完成并能夠取得測值。此外，氣象條件由于影響作用范圍有限，直接采用壩址區(qū)測站的數(shù)值。

獲取水庫水溫的實測值是本次研究的重要內(nèi)容，為獲得準確的壩前水溫值，在大壩的31個壩段中選取6#壩段、10#壩段、16#壩段、22#壩段和27#壩段共5個典型壩段來實地測量其壩面上游的水溫，可正常使用的溫度計共108 支，測點最低高程373 m，最高高程604 m 高程，具體布置位置如圖4所示。溫度計采用高精度電阻測溫計，測得水溫精度能夠達到0.01℃。從2012年3月1日開始取得測量初值，每12 h 取一次讀數(shù)。

溫度換算公式：

式中：R0為0℃時的電阻，電阻值為46.60 Ω；Rt為任意時刻實測溫度下阻值；K為溫度系數(shù)，單位為℃/Ω。

圖4 壩前溫度計布置位置示意

圖5 多年平均水溫與壩前實測水溫對比圖

圖6 工程設(shè)計出、入庫流量與實際出、入庫流量對比

與工程可行性研究階段的多年統(tǒng)計數(shù)據(jù)相比，進、出庫流量，來水水溫等資料均出現(xiàn)明顯的變化。圖5—6是初步設(shè)計階段和現(xiàn)階段來水溫度、水庫進、出流量等關(guān)鍵因素的對比，可以看出，大壩實際蓄水時上游來水水溫在某些月份明顯高于多年平均統(tǒng)計河水水溫，實際流量與設(shè)計流量過程也有較大差別。

3.4 計算結(jié)果對比分析圖7所示的壩前水溫溫度計的實測值和計算值變化過程線對比可知，采用實際監(jiān)測資料對庫區(qū)水溫模擬的計算結(jié)果基本上能夠反應(yīng)出實際水溫變化的趨勢，多數(shù)情況下與實測值吻合較好（大壩壩前堆渣至370 m 高程，2012年3月15日基坑開始進水，2012年7月到400 m 高程，2013年1月水位450 m 高程，2013年6月水位540 m 高程，2014年10月蓄水至600 m 高程）。相對而言，在蓄水的前期水位較低、來水湍流現(xiàn)象突出、不同來水水溫混摻過程較為劇烈時，模擬結(jié)果與實測值的誤差值相對較大。2013年水位升高之后的計算基本上能夠反映真實情況，平均誤差不超過1℃。每年的平均誤差可按照下式來計算：

得到下表：

表1 蓄水后3年平均誤差 （單位：℃）

從表1可知，2013—2015年平均誤差逐年減小，可以推斷隨著水位蓄至正常水位后，水流的波動劇烈程度減小，同時各種因素的可控制性增加，計算方法中對于二維計算的基本假設(shè)和相應(yīng)的湍流模型更加接近真實狀態(tài)。

由上述分析可知，所建的河道模型能夠反映溪洛渡庫區(qū)水溫的變化規(guī)律，模型的適用性和可信度較高，隨著水位的相對穩(wěn)定和計算時間的增加，計算的平均誤差會進一步減小。

4 結(jié)果與討論

基于上文的計算分析成果，可對該類型水庫從蓄水到初期運行階段的水溫演化規(guī)律做進一步的分析研究。

圖7 典型高程溫度過程線

（1）3種數(shù)值方法的差異。如下圖8中所示，由2015年10月19日實測值與3種計算方法計算結(jié)果的對比可知，經(jīng)驗公式獲得的水溫分布情況與實際值差異較大，表明對于高拱壩而言采用該公式是不合適的；設(shè)計階段采用基于多年平均統(tǒng)計來水水溫、流量等參數(shù)的一維算法獲取的水溫分布規(guī)律與實際水溫分布規(guī)律類似，但數(shù)值誤差較大；本文采用邊界條件精確處理后的CE-QUAL-W2模型，并且考慮了蓄水過程，所獲得的水溫分布與變化規(guī)律與實際情況吻合更好。圖中實測值是2015年10月19日，水位600 m時的實測值，方法1為朱伯芳經(jīng)驗公式法，方法2為設(shè)計階段一維算法，方法3為本文采用的CE-QUAL-W2模型算法。

（2）水溫垂向分布情況對比分析。本文關(guān)于水庫水溫的跟蹤模擬計算從2012年3月圍堰拆除開始到2015年10月截止，歷時3年8個月。各個典型時刻沿高程方向水溫分布情況如圖9所示，圖9（a）為2013年6月20日首次蓄水至540 m 高程，由于在該階段的蓄水過程中，50 d 內(nèi)水位上升了100 m，庫容增加了53 億m3，該過程水體充分摻混，因而未出現(xiàn)明顯的分層現(xiàn)象；圖9（b）為蓄水進行了18個月后，第一次蓄水至最高水位600 m 高程，水溫的垂向分布已經(jīng)出現(xiàn)了明顯的分層，與以往認為水庫水溫分層需要很多年的認知有所不同，同時由圖可知，剛蓄水至最高水位時，計算與實際的溫度值基本吻合，但分層位置仍稍有差異，隨著成庫后水體的進一步穩(wěn)定，這種差異逐漸減小，如圖9（c）（d）所示，可見數(shù)值模擬的過程與實際的物理過程能夠基本吻合。其中，圖9（a）（d）中470～530 m 高程之間實測值波動較大的原因是受到汛期孔口泄流的影響。

圖8 實測值與3種方法對比

圖9 典型時刻不同壩段上游溫度測量值與計算值對比

（3）蓄水初期前3年水庫水溫變化規(guī)律。由計算結(jié)果可知，溪洛渡水庫水溫基本在蓄至正常水位后第二年3月份開始出現(xiàn)分層，隨后水溫的空間分布隨流量和溫度等因素的變化而變化，蓄水初期水庫水溫變化情況見圖10。由圖可知，水體的溫度變化基本上可以分為8個典型區(qū)域：在接近水面區(qū)域，由于受來水、氣溫和太陽輻射的影響，溫度隨時間變化波動較大，溫度年變幅約13℃，典型溫度過程線如圖11中536 m 高程過程線所示；在接近水面區(qū)域以下，水庫中存在大面積的溫度過渡區(qū)，在上下層溫度傳導(dǎo)過程中起到過渡的作用，對比圖11中4 條過程線可以明顯看出溫度層層傳遞并逐漸減弱的過程，從水庫水位上升至560 m 高程后的次年春季開始，536 m 高程溫度的極大值出現(xiàn)在7、8月份，475和431 m 高程溫度的極大值出現(xiàn)時間逐漸后移，并且量值減小，374 m 高程的溫度過程線幾乎為水平，過渡區(qū)域溫度年變幅基本在3.2～10.2℃左右，高程越低變幅越小；從2015年以后，庫底水溫基本處于相對穩(wěn)定狀態(tài)，年變幅小于1℃，該區(qū)域溫度基本穩(wěn)定在某個溫度值附近，上部邊界的影響很難傳播到該區(qū)域，典型情況如圖11中374 m 高程過程線所示。

（4）水庫水溫長期運行規(guī)律。為了解溪洛渡多年運行后水庫水溫可能的變化趨勢，模擬從最初蓄水一直計算至2040年12月的計算工況，其中，2016—2020年的輸入數(shù)據(jù)分布采用2013年和2014年的數(shù)據(jù)（即來水水溫一年低，一年高，其他參數(shù)保持不變），2021—2040年后，均采用2014年10月—2015年10月的實測數(shù)據(jù)，得到如下圖12所示的結(jié)果。

圖10 蓄水后水庫水溫變化云圖

圖11 16#壩段不同高程上游表面溫度計溫度變化過程線

溪洛渡的水體特點是水深較深，同時庫容較大，這在一定程度上有利于水體的穩(wěn)定，尤其對于下層而言，除與巖基發(fā)生熱交換之外基本不受外部影響。底部水體在趨穩(wěn)的同時，由于受上層和基礎(chǔ)傳導(dǎo)的熱量的影響，當上游來水溫度稍高時，溫度會略有回升，而當1月份上游來水溫度低于底部溫度時，1月份的低溫水會直接沉入底部，使底部溫度重新回到較低的狀態(tài)，即底部庫水溫度將隨著上游來水溫度的變化輕微波動，但總體處于相對穩(wěn)定狀態(tài)，而上部水體溫度則將如圖12所示呈現(xiàn)出不同程度的波動。

從多年的預(yù)測結(jié)果可以看出，溪洛渡水體會在短期內(nèi)出現(xiàn)分層并有一個相對穩(wěn)定的狀態(tài)，但是極易受年際上游來水水溫變化的影響，沒有絕對的穩(wěn)定狀態(tài)，尤其是上部淺層區(qū)域。

由上述分析可以推測，接下來的運行階段，若不遇到極端年份，各高程的水溫將會隨時間每年呈近似周期性變化，并且隨著深度的增加年變幅減小。上部520～600 m 高程范圍內(nèi)，隨季節(jié)和水位的改變，溫度變化較為明顯，變動范圍在15～27℃之間，年變幅為12℃；450～520 m 高程的70 m 范圍內(nèi)，水體傳熱作用層層減弱；450 m 高程至庫底的范圍內(nèi)基本上常年穩(wěn)定在15℃左右，接近冬季最低河水溫度。

圖12 374m、450m和520m 高程溫度長期預(yù)測過程線

5 結(jié)論

本文主要取得以下幾個方面的結(jié)論：（1）本文在高壩大庫中引入二維橫向平均的水動力模型CE-QUAL-W2，按照真實的邊界條件對該水庫從開始蓄水到初期運行歷時近4年的溫度變化過程進行了模擬分析，得到的計算結(jié)果與實測數(shù)據(jù)吻合較好，驗證了該模型對于溪洛渡水庫的適用性。（2）從蓄水到初期運行的整個過程中，當水位較淺時，水體的水溫變化規(guī)律依然表現(xiàn)為天然河流的特性，當水位逐漸增加至水庫正常蓄水位的過程中，水體的溫度開始逐漸表現(xiàn)為水庫的特性。研究發(fā)現(xiàn)，成庫后約半年時間后，溪洛渡水庫能形成4個典型的區(qū)域，每個區(qū)域的溫度年變幅不同，其中過渡區(qū)區(qū)域最大，年內(nèi)不同特征位置的水溫極大值出現(xiàn)時間逐漸后移并減小，原來出現(xiàn)在高溫季節(jié)的水溫極大值后移至年末，至接近水庫底部時年變幅幾乎為零；此外，水庫長期運行的預(yù)測研究表明，庫底溫度除受上部水體的微弱影響之外，還主要受冬季低溫水體的影響，如果遇到極端年份，則有可能會打破庫底水溫原有的平衡。（3）中國還有一批類似溪洛渡拱壩的其他300 m 級的特高拱壩在建、擬建或規(guī)劃建設(shè)，如白鶴灘、烏東德及松塔等電站。本文多年水庫水溫跟蹤仿真分析的成果及水溫變化分布規(guī)律研究可以為這些電站的大壩溫控防裂及結(jié)構(gòu)設(shè)計提供有益的借鑒。