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

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

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

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

1 數(shù)學(xué)模型

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

(1)連續(xù)方程

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

(3)溫度

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

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

2 數(shù)值計算及分析

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

2.1 網(wǎng)格劃分

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

圖1 水庫概化及計算區(qū)域單元格

2.2 初始條件

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

2.3 邊界條件

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

圖2 入流流量過程線

圖3 出流流量過程線

2.4 模型驗證

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

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

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

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

2.5 模擬結(jié)果及分析

圖7 壩前斷面的溫度分層

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

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

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

圖8 輸水洞斷面水溫分布

3 結(jié)論

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

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

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

［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］李云生，劉偉江，吳悅穎，等.美國水質(zhì)模型研究進展綜述［J］.水利水電技術(shù)，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］.四川大學(xué)學(xué)報，2009，41(1):30－33.

［7］張士杰，彭文啟.二維水庫水溫結(jié)構(gòu)及影響因素研究［J］.水利學(xué)報，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］申滿斌，陳永燦，劉昭偉，等.岸邊排放污染物濃度場三維渾水水質(zhì)模型研究［J］.水力發(fā)電學(xué)報，2005，4(3):93－98.

［10］郭磊，高學(xué)平，張晨，等.北大港水庫水質(zhì)模擬及分析［J］.長江流域資源與環(huán)境，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］.水動力學(xué)研究與進展，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］李錦繡，杜斌，孫以三.水動力條件對富營養(yǎng)化影響規(guī)律探討［J］.水利水電技術(shù)，2005(5):23－28