磨盤山水庫溫度場的數值模擬

伍悅濱，徐 瑩，田 禹，張海龍

(1.哈爾濱工業大學城市水資源與水環境國家重點實驗室，哈爾濱150090，ybwu@hit.edu.cn; 2.哈爾濱工業大學土木工程學院，哈爾濱150090;3.哈爾濱供水工程有限責任公司，哈爾濱150076)

水庫蓄水帶來了防洪、發電、灌溉，供水等綜合效益，同時也引起了河流水文、泥沙、生態等各方面的環境影響，隨著社會經濟的發展，大量工業廢水和城市污水未經處理排入河流、水庫，超過了水體的自凈能力，河流，水庫污染日益嚴重，影響著人類的生存和發展.水溫是水質的一個重要影響因素，因此進行水庫水溫分析和預測，是水庫環境評價、規劃、治理和管理的基礎.

磨盤山水庫位于拉林河干流上游五常市沙河子鄉沈家營村的上游1.8 km附近，拉林河是松花江右岸的一級支流，流域總面積19200 km2.除供水外，磨盤山水庫兼有防洪、灌溉及其他綜合利用效益.水庫的建成將改變庫區及大壩下游河段的水文和水環境狀況，進而對水庫下游河段的灌溉和生活用水產生重大影響.

1 數學模型

EFDC(Environmental Fluid Dynamics Code)是由美國國家環保署資助開發，用于模擬河流、湖泊、水庫、海灣、濕地和河口等水環境系統中流場、溫度場、濃度場及生態過程的三維數值計算模型［1－2］.對水庫水溫進行模擬的數學模型主要采用Mellor等［3－5］提出的控制方程.

(1)連續方程

(2)曲線正交坐標中x向和y向動量方程

(3)溫度

式中:u和v分別是沿著正交曲線水平方向x向和y向的流速分量;mx和my是平面坐標轉換因子，m=mxmy.

其中，z為無量綱垂向坐標，w*為物理垂向速度;水深H=h+ζ，h代表實際垂向的底高程;ζ代表水位;動量方程中f表示科里奧利(Coriolis)參數; Av為垂向紊動粘性系數;Qu和Qv表示動量方程中的源匯項;θ為水溫;在溫度方程中的Qθ項是外源輸入輸出項.連續方程中的QSS和QSW表示河床與水體之間的泥沙及水的交換通量，其正值表示通量由河床入水體中，反之為負值.三角函數δ(0)表示這類通量進入水體的最底層.

2 數值計算及分析

磨盤山水庫為典型的河道型水庫，干流回水長約13 km.地形總趨勢為東南高、西北低，呈不對稱的“U”字型，地面坡降約為3.8%.磨盤山水庫壩高45 m，正常蓄水位318 m，相應的水庫總庫容355.92×106m3，調節庫容為292.95×106m3.汛期限制水位為 317.00 m，重復利用庫容為27.27×106m3.設計洪水位為318.77 m，相應的設計總庫容為375.98×106m3.校核洪水位為320.89 m，相應的校核總庫容為443.02×106m3.

2.1 網格劃分

庫區地形資料是進行水庫模擬計算的基礎，首先根據水庫地形高程圖，按照模型要求劃分計算網格.計算中時間步長為10 s，運行周期為1 d，即86 400 s，模型運行時間設定為30個周期.控制方程的數值解法為有限體積法和有限差分法的結合，對于模擬對象，沿流動方向采用正交曲線網格離散，而在垂直方向上則采用分層方法進行離散，進而引入了層厚的概念.根據磨盤山水庫的地理特點，將其水平方向上劃分1 052個矩形網格，每個單元格大小為85 m×65 m.考慮溫度分層效應，垂直方向上分為8層，占水深比例分別為12.5%，12.5%，20%，20%，20%，5%，5%，5%.水庫概化及計算區域單元格見圖1.

圖1 水庫概化及計算區域單元格

2.2 初始條件

拉林河流域內現有沈家營、五常等13處水文站，根據磨盤山水庫工程地點和資料情況，選用的主要分析測站是拉林河的沈家營測站.測站觀測項目有降水、蒸發、氣溫、氣壓、日照、溫度和風速等以及磨盤山庫區的日氣溫變化.初始水位水深根據地形圖得到，初始水溫根據水庫管理區測量資料以及《水利水電工程水文計算規范》(SDJ214－83)(試行)條例中給出.地面坡降為3.8%，初始計算流速設為零.

2.3 邊界條件

根據庫區管理站測量的2007年入庫流量和出庫流量作為上、下游邊界條件，入庫流量過程線見圖2，出庫流量過程線見圖3.

圖2 入流流量過程線

圖3 出流流量過程線

2.4 模型驗證

通過對磨盤山水庫枯水期、平水期和豐水期的模擬，得到了水庫出流溫度變化曲線，并與實際出流溫度變化進行了比較，結果如圖4～6所示.模擬結果與實測值基本吻合，說明計算模型及方法的可信性.

圖4 2007年4月壩前水溫實測與模擬值對比

圖5 2007年6月壩前水溫實測與模擬值對比

圖6 2007年9月壩前水溫實測與模擬值對比

2.5 模擬結果及分析

圖7 壩前斷面的溫度分層

本文分別模擬了2007年枯水期(四月)、平水期(六月)、豐水期(九月)磨盤山水庫的水溫變化規律.壩前水溫在枯水期、平水期、豐水期的溫度場計算結果分別見圖7.水溫在4月份為升溫期，入流水溫和表面熱通量逐漸增加，上層水體水溫增長迅速;在4月末，表層水溫已達到10℃，但是靠近底部的水溫幾乎沒有變化，仍保持低于6℃的低溫，在水面附近出現了較大的溫度梯度，形成了較明顯的水溫分層.6月份為高溫期，入流水溫和氣溫都達到最高.在表面很薄的水體中，由于水氣熱交換形成了表層溫躍層［6－7］，表層水溫保持在15℃作用，底部的水溫緩慢上升了2～3℃，但仍然保持低于8℃的溫度，在整個垂向斷面上的溫差達到了接近10℃的水平，在這一時期，水體上部分保持穩定的溫度分層現象.9月份為降溫期，氣溫和入流水溫逐漸下降，水體向大氣散失熱量，冷水下沉使表層溫躍層消失，水庫在整體垂直方向形成水溫分層，并開始出現層與層之間的垂向摻混，即秋季翻滾現象［8－10］.

水庫的這種分層結構對水質具有重要的影響.對于磨盤山水庫而言，由于其溫度分層作用，不同深度的水質狀況將有所不同.在枯水期，水庫水溫混合基本均勻，導致水中的營養物質沿垂向均勻混合.又加之冬季水中生化反應緩慢，消耗較少的營養物質，至春季時沉積了大量營養物質，為藻類快速生長提供了物質基礎.平水期處于夏季，水庫上部分水體水溫分層較為明顯.此時溫躍層在垂直方向上具有“屏障”作用［11－13］，水中的營養物質在垂直方向上的交換量被大大的削減.在豐水期，汛期來臨，洪水挾帶著沖起的河底淤沙，涌進水庫，使大量的亞磷酸鹽進入水流.水庫水溫的穩定分層被破壞后，“屏障”消失，水庫底層的沉淀物和有害物質會翻滾到表層，對水體造成污染.

圖8是枯水期、平水期和豐水期在水庫輸水洞斷面的模擬水溫分布剖面圖.磨盤山水庫的輸水洞位于右岸，進口位于右岸壩軸線上游約250 m處.水體自輸水洞進入長直管線，向凈水廠輸送原水.由模擬結果可知，水庫在枯水期和平水期沿深度方向上水溫分布穩定，水溫較低，可以直接選取水面以下10 m的水庫水作為飲用水水源;在豐水期，水溫結構發生變化.隨著氣溫降低，表層水溫也下降，相應的水密度增大，導致表層水開始下沉，底層水上升，在整個深度范圍內表層和底層發生了混合，這種秋季翻滾現象可能會引發水質問題［14－15］.

圖8 輸水洞斷面水溫分布

3 結論

1)通過實測值與模擬值的對比驗證，應用EFDC建立的水庫水動力學模型，能夠很好的模擬水庫中水溫的時空分布.因此基于該模型進行水體中流場、溫度場、濃度場及生態過程的數值模擬是進行水庫環境保護、管理及運行的有力工具.

2)磨盤山水庫屬于分層型水庫，水溫出現季節性分層現象，庫表分層現象最為明顯，至庫底分層減弱，庫底水溫變化平緩.枯水期垂向分層現象不明顯，平水期形成表層溫躍層，對于垂向傳質具有一定程度的屏障作用，豐水期出現秋季翻滾，對水質影響較大.

3)水庫取水口深度應該隨季節調整.在枯水期和平水期沿深度方向上水溫分布穩定，大部分水體水溫較低，水深10 m以下均可布置取水口;而在豐水期出現的秋季翻滾現象，形成水體垂向的自循環流動，可能引發水質問題.因此，豐水期至冬季時段，取水口位置應有所調整，并建議凈水廠強化水處理工藝.

［1］HAMRICKJ M.A Three-dimensional environmental fluid dynamics computer code:Theoretical and computational aspects［M］.Gloucester，Massachusetts:Virgina Institute of Marine Science，the College of William and Mary，1992:1－63.

［2］李云生，劉偉江，吳悅穎，等.美國水質模型研究進展綜述［J］.水利水電技術，2006(2):68－73.

［3］MELLOR G L，BLUMBERG A F.Modeling vertical and horizontal diffusivities with the sigma coordinate system［J］.Mon Wea Rev，1985，20(5):851－875.

［4］MELLOR G L，YAMADA T.Development of turbulence closure model for geophysical fluid problems［J］.Rev Geophysical Space Phys，1982，20(6):851－875.

［5］ARAKAWA A，LAMBV R.Computational design of the basic dynamical processes of the UCLA general circulation model［M］//CHANG J.Methods in Computational Physics.New York:Academic Press，1977:173－265.

［6］梁瑞峰，李嘉，李克峰，等.立面二維水庫水溫模型的并行研究［J］.四川大學學報，2009，41(1):30－33.

［7］張士杰，彭文啟.二維水庫水溫結構及影響因素研究［J］.水利學報，2009，40(10):1254－1259.

［8］OUGOLNITSKY G A，USOV A B.Water quality control in watercourses［J］.Water Resources，2003，30(2): 226－232.

［9］申滿斌，陳永燦，劉昭偉，等.岸邊排放污染物濃度場三維渾水水質模型研究［J］.水力發電學報，2005，4(3):93－98.

［10］郭磊，高學平，張晨，等.北大港水庫水質模擬及分析［J］.長江流域資源與環境，2007，16(1):11－16.

［11］WOOL T A，DAVIE S R，HUGO N R.Development of three-dimensional hydrodynamic and water quality models to support TMDL decision process for the Neuse River Estuary，North Carolina［J］.Journal of Water Resources Planning and Management，2003，129(4):295－306.

［12］JIN K R，HAMRICK J H，TISDALE T.Application of a three-dimensional hydrodynamic model for lake okeechobee［J］. Journal of Hydranlic Engineering ASCE，2001，126:13－14.

［13］陳景秋，趙萬星，季振剛.重慶兩江匯流水動力模型［J］.水動力學研究與進展，2005，12:829－835.

［14］JOHNSON B H，KIM K W，HEATH R E，et al.Validation of three-dimensional hydrodynamic model of chesapeake bay［J］.J Hyd Engrg，2003，119:2－20.

［15］李錦繡，杜斌，孫以三.水動力條件對富營養化影響規律探討［J］.水利水電技術，2005(5):23－28