朱子晨，胡澤建，張永強(qiáng)，遲萬清，劉建強(qiáng)，邊淑華，熊叢博

（國家海洋局 第一海洋研究所，山東 青島 266061）

近年來我國逐步加強(qiáng)海岸帶開發(fā)利用活動(dòng)，在沿海地區(qū)興建了許多熱電廠，其具有就近利用海水作為循環(huán)冷卻水源等優(yōu)勢，但同時(shí)向海灣排放溫排水，造成海水溫度時(shí)間、空間分布變化，引發(fā)水體中溶解氧濃度降低，富營養(yǎng)化程度加深，藻類生物資源種群及數(shù)量變化等問題。針對溫排水海區(qū)海水溫度變化的研究，對搞清溫排水熱污染狀況等課題具有重要意義。目前，溫排水?dāng)?shù)值研究大多數(shù)基于二維模型，韓康等對近海海灣溫排水造成的環(huán)境影響情況使用二維數(shù)值模型進(jìn)行了模擬［1］，研究了多種工況下溫排水海區(qū)溫升場的特征；張細(xì)兵等對感潮河段溫排水環(huán)境影響進(jìn)行了二維數(shù)值模擬［2］，并使用實(shí)體模型進(jìn)行了比對，研究了當(dāng)?shù)販嘏潘^(qū)溫升場特征。然而二維數(shù)值模型無法反映溫度垂向分布不均勻的特征。曹穎等也有對半封閉海灣進(jìn)行過三維數(shù)值研究［3］，得到了溫排水影響的三維結(jié)構(gòu)，但現(xiàn)階段就此類問題而言，結(jié)合觀測資料的研究較為缺乏［4］。因此，本文通過建立三維海洋數(shù)值模型，結(jié)合海上觀測資料，對受溫排水影響的海區(qū)海水溫度場特征進(jìn)行研究。

1 三維非結(jié)構(gòu)網(wǎng)格有限體積海洋數(shù)值模型

1.1 水動(dòng)力模型

FVCOM模型是美國Massachusetts Dartmouth州立大學(xué)陳長勝領(lǐng)導(dǎo)的研究小組開發(fā)的河口、海岸模型，采用有限體積法，垂向采用σ坐標(biāo)。σ坐標(biāo)變換如下：

σ坐標(biāo)下海水運(yùn)動(dòng)方程、連續(xù)方程及狀態(tài)方程為：

在σ坐標(biāo)系下，水平擴(kuò)散項(xiàng)可以使用如下定義進(jìn)行近似：

式中，σ、z分別為σ坐標(biāo)系與直角坐標(biāo)系下垂向坐標(biāo)；ζ為潮位，H為平均海平面下水深值，D為總水深；t為時(shí)間，u、v、ω分別為x、y、σ方向的流速分量；ρ與ρ0為海水的密度與參考密度，g為重力加速度，f為科氏力參量，Am和Km為水平及垂向渦動(dòng)粘滯系數(shù)。

1.2 溫度對流擴(kuò)散及熱通量方程

σ坐標(biāo)系下，三維溫度對流擴(kuò)散模型為：

式中，T為溫度；Kh、Ah分別為垂向、水平熱擴(kuò)散系數(shù)；＾H為短波輻射熱通量。

1.3 數(shù)值計(jì)算方法

FVCOM在計(jì)算過程中，使用內(nèi)外模分開求解的方式，外模采用有限體積計(jì)算方法，使動(dòng)量方程、連續(xù)方程在三角形網(wǎng)格區(qū)域內(nèi)積分后，由改進(jìn)的4階龍格庫塔法求解。內(nèi)模采用顯、隱式結(jié)合的差分格式，局部流速變換使用迎風(fēng)格式，對流項(xiàng)采用2階龍格庫塔時(shí)間推進(jìn)格式，垂向擴(kuò)散采用隱式求解。溫度計(jì)算方法與潮流類似。模型在海陸邊界采用動(dòng)邊界處理方式。

2 模型的建立與驗(yàn)證

2.1 模型構(gòu)建

模型水深數(shù)據(jù)來自解放軍司令部航海保證部2005版10月第一版湄洲灣1：40000海圖（圖號(hào)14171），計(jì)算開邊界選取于湄洲灣外開闊海域。開邊界采用水位驅(qū)動(dòng)，輸入水位由崇武、平海鎮(zhèn)連續(xù)一個(gè)月觀測資料調(diào)和分析所得到的調(diào)和常數(shù)推算并插值得到［5］。

水深地形、觀測站位位置及計(jì)算區(qū)域非結(jié)構(gòu)網(wǎng)格剖分見圖1。

圖1 計(jì)算區(qū)域水深地形、觀測站位位置及網(wǎng)格剖分Fig.1 Bathymetry in computational domain，locations of observation stations and mesh dissection of finite element

模型主要參數(shù)設(shè)置如下：

時(shí)間步長：5s； 總運(yùn)行步長數(shù)目：276 480（即計(jì)算時(shí)間15d）；

上午，我到老師辦公室領(lǐng)書去，呀，徐老師真是一點(diǎn)也沒變，依然是那么和藹可親，那么美。她看見了我，笑瞇瞇地說：“以和呀，想我了沒？”我甜甜地回答：“想?！崩蠋熜α?，露出了兩個(gè)小酒窩。

底部應(yīng)力阻力系數(shù)：0.04； 底摩擦層高度：0.003m；

網(wǎng)格節(jié)點(diǎn)總數(shù)：11 664； 網(wǎng)格單元總數(shù)：22 071；

垂向σ坐標(biāo)分層：5； 電廠排水口循環(huán)冷卻水排量：44.23m3／s；

電廠排水口溫度：35.0℃。

2.2 計(jì)算結(jié)果驗(yàn)證

采用圖1中所示的5個(gè)站位連續(xù)24h海流同步觀測資料垂向平均流速、流向，1個(gè)驗(yàn)潮站24h連續(xù)潮位觀測資料，對模型計(jì)算結(jié)果進(jìn)行驗(yàn)證，觀測時(shí)間為2007－08－14T19：00—08－15T19：00，驗(yàn)證結(jié)果（圖2）表明，模型較好地模擬出了觀測期間各站海流情況。

圖2 各站位物理參數(shù)驗(yàn)證曲線Fig.2 Validation curves of physical parameters at individual stations

3 垂向溫度不均勻特征分析

3.1 湄洲灣海區(qū)流場特征與溫排水影響范圍

觀測資料與數(shù)值模型計(jì)算結(jié)果均顯示，湄洲灣海區(qū)潮流場呈明顯往復(fù)流態(tài)，流速較大，漲落急時(shí)刻流速均達(dá)到60cm／s以上。中央深槽附近落潮流速達(dá)到110cm／s，且落潮流速大于漲潮流速。由于表層溫排水影響范圍較中層、下層廣，因此以表層海水溫度場平面分布代表溫排水影響范圍［6］，圖3所示的為與潮流矢量圖相對應(yīng)的漲、落急時(shí)刻表層海水溫度分布。漲落急時(shí)刻溫度分布情況一定程度上反映了溫排水影響范圍，但溫排水最大影響范圍發(fā)生在漲、落憩時(shí)刻。

圖3 漲、落急時(shí)刻潮流矢量圖及海水溫度分布Fig.3 Vector diagrams of tidal current and distribution of water temperature at maximum flood and maximum ebb

3.2 觀測資料與計(jì)算結(jié)果中海水溫度變化趨勢分析

從圖4中5站同步溫度觀測資料來看，距離溫排水排放口較遠(yuǎn)的L2、L4站除轉(zhuǎn)流時(shí)刻外，24h內(nèi)溫度垂向分布較均勻，表層與0.6H層、底層基本一致，沒有明顯差異，僅隨潮流場呈周期性變化，這與近岸海域水深相對較淺，垂向混合均勻有關(guān)。但上述站位普遍在轉(zhuǎn)流時(shí)刻（漲憩、落憩時(shí)刻）1h后出現(xiàn)不同程度的垂向溫度不均勻特征，這種特征持續(xù)時(shí)間約為3h。在距離排放口較近的L1、L3站同樣在轉(zhuǎn)流時(shí)刻出現(xiàn)垂向溫度不均勻分布，持續(xù)時(shí)間較長，約為6～7h。與海上觀測資料相對應(yīng)，數(shù)值模型同樣模擬出了這一特征。此外，L5站觀測資料說明，在漲憩時(shí)刻后也出現(xiàn)約持續(xù)3h的溫度垂向分布不均特征。

圖4 各站位計(jì)算、觀測海水溫度Fig.4 Water temperatures calculated and observed at individual stations

3.3 特征潮時(shí)表底層海水溫差值平面分布

為了表現(xiàn)海區(qū)垂向水溫不均勻程度，使用數(shù)值模型計(jì)算區(qū)域表層水溫與底層水溫之差，得到表底層水溫差在海區(qū)的平面分布，如圖5所示。

圖5 各潮時(shí)表底層海水溫差平面分布Fig.5 Horizontal distribution of temperature difference between the surface and the bottom waters at each time of tide

3.4 斷面溫度變化特征

由于漲憩時(shí)刻S1－S2斷面溫差較大，因此以漲憩時(shí)刻該斷面垂向溫度分布代表垂直于海區(qū)主流向斷面的變化趨勢。另外，以漲急時(shí)刻S1－S3斷面及落急時(shí)刻S1－S4斷面溫度垂向分布代表與海流流向相近斷面的變化情況。

從S1－S2斷面垂向溫度分布情況來看，在排水口附近、中央深槽內(nèi)，溫度垂向上差異較大。S1－S3斷面則表明，溫排水沿漲潮流方向運(yùn)動(dòng)，且高溫水團(tuán)集中在海水表層，底層溫度相對較低。S1－S4斷面溫度垂向分布則說明，表層高溫海水沿落潮流方向運(yùn)動(dòng)較底層快，隨著海水向水深增加的區(qū)域運(yùn)動(dòng)，垂向溫差逐漸顯現(xiàn)。

圖6 斷面位置圖Fig.6 Locations of the sections

圖7 S1－S2斷面漲憩時(shí)刻溫度垂向分布Fig.7 Vertical distribution of water temperature along Section S1－S2at high water level

圖8 S1－S3斷面漲急時(shí)刻溫度垂向分布Fig.8 Vertical distribution of water temperature along Section S1－S3at maximum flood

圖9 S1－S4斷面落急時(shí)刻溫度垂向分布Fig.9 Vertical distribution of water temperature along Section S1－S4at maximum ebb

3.5 垂向溫度分布不均勻的成因分析

從觀測資料及數(shù)值模擬結(jié)果來看，電廠溫排水是湄洲灣海區(qū)溫度垂向分布不均勻的主要誘因。除此以外，這種不均勻分布特征在空間上與水深地形具有較大關(guān)系，在時(shí)間上與潮時(shí)具有相關(guān)性。

3.5.1 與排水口相對位置對垂向溫度分布的影響

數(shù)值模型計(jì)算結(jié)果表明，出現(xiàn)垂向溫度不均勻的區(qū)域在漲潮、落潮過程中，分別集中出現(xiàn)在排水口北側(cè)、南側(cè)，由于排水口附近平面溫度梯度最大，容易在潮流作用下形成高溫水舌［7］，且高溫海水密度較低，趨于在海水表層運(yùn)動(dòng)，因此排水口附近始終存在垂向溫差，且溫差高于其他區(qū)域。海上觀測資料基本與模型計(jì)算結(jié)果一致，距離排水口較近的L1、L3站位，相比于L2、L4、L5站，溫差更加顯著，且這種溫度在垂向上分布不均勻的特征持續(xù)時(shí)間較長。

3.5.2 水深地形與潮流場共同作用下垂向海水溫差在平面上的分布

圖7～圖9表明，表底層海水溫度差在平面上的分布與水深地形有密切聯(lián)系，在不考慮鹽度影響的條件下，海水溫度越高，其密度越低，低密度海水趨向于在海水上層運(yùn)動(dòng)，而下層海水溫度變化，主要依靠熱擴(kuò)散與熱傳導(dǎo)，變化相對緩慢。同時(shí)，由于潮流場作用，表層高溫水團(tuán)的分布范圍沿流向增大，在更大范圍區(qū)域體現(xiàn)了各層海水溫度變化的差異。因?yàn)樯鲜鲞\(yùn)動(dòng)趨勢，當(dāng)水深較淺處的高溫水團(tuán)經(jīng)過水深梯度較大的區(qū)域，向深水區(qū)域運(yùn)動(dòng)時(shí)，高溫海水趨向于在上層運(yùn)動(dòng)，而深水區(qū)水溫變化相對緩慢，將顯著出現(xiàn)水溫在垂向上的不均勻分布。斷面S1－S2漲憩時(shí)刻垂向溫度分布說明，當(dāng)斷面與流速基本成垂直趨勢時(shí)，由于高溫水體向水溫較低的深水區(qū)擴(kuò)散，將產(chǎn)生垂向溫差，差異相對較弱，而斷面S1－S3、S1－S4則說明，高溫水團(tuán)沿流向向水深較深處運(yùn)動(dòng)時(shí)，將產(chǎn)生顯著的垂向溫度差異。

湄洲灣海區(qū)島嶼與水道交錯(cuò)，有多處水深明顯變化區(qū)域，且水動(dòng)力條件較強(qiáng)。羅嶼南側(cè)水域位于中央深槽以東，自中央深槽至該水域，水深迅速變化，因此該區(qū)域表底層存在顯著水溫差異，從觀測資料來看，L3站位于水道，與數(shù)值模型相吻合。中央深槽向東南延伸，經(jīng)過盤嶼、大竹島時(shí)分為三條支水道，分別為大竹島西南側(cè)水道、盤嶼北側(cè)水道及大竹島嶼盤嶼之間的支水道，三條支水道與各自北側(cè)海區(qū)均具有明顯的水深梯度。排水口附近水深較淺，溫排水在隨落潮流向南運(yùn)動(dòng)時(shí)，首先經(jīng)過盤嶼北側(cè)支水道，由于水深發(fā)生較大變化，海水溫度垂向差異明顯。當(dāng)溫排水隨落潮流繼續(xù)向南運(yùn)動(dòng)并到達(dá)大竹島南側(cè)水道及大竹島嶼盤嶼之間支水道后，由于水道內(nèi)水深顯著變化，高溫水團(tuán)在表層形成鋒面，與底層溫度變化緩慢的海水形成溫度差異。但因上述兩條支水道距離排水口較遠(yuǎn)，溫排水僅在落憩時(shí)刻前后能夠傳遞至水道口門，該區(qū)域出現(xiàn)垂向溫度不均勻特征的時(shí)段較短。對于中央深槽北部而言，由數(shù)值模型得到的L2站各層溫度變化曲線也說明，約在落憩時(shí)刻1h后，由于漲潮流作用，該站位出現(xiàn)了表層水溫高于底層水溫的模擬結(jié)果。

4 結(jié) 論

通過數(shù)值模型與觀測資料相結(jié)合的研究，本研究得出了受溫排水影響下，湄洲灣局部海區(qū)垂向水溫不均勻特征的分布情況與變化規(guī)律：

1）由于湄洲灣溫排水熱源的影響，近岸淺海水域原本混合均勻的水體，垂向上可能產(chǎn)生溫度差異，當(dāng)潮流場作用程度與水深梯度較小時(shí)，這種垂向溫度不均勻特征，在距離排水口越近的區(qū)域越顯著，持續(xù)時(shí)間也越長。

2）溫排水隨潮流運(yùn)動(dòng)，當(dāng)潮流流向與等深線交叉，且水深梯度明顯時(shí)，高溫水團(tuán)將于很長一段時(shí)間內(nèi)持續(xù)在上層海水運(yùn)動(dòng)，水溫在垂向上呈現(xiàn)出顯著不均勻特征。

3）垂向溫度不均勻特征的分布范圍隨流場持續(xù)作用積累，并表現(xiàn)為沿流向運(yùn)動(dòng)，在轉(zhuǎn)流時(shí)刻前分布范圍積累到最大面積。

（References）：

［1］ HAN K，ZHANG C Z.Numerical simulation of warm water disposal in Sanya Power Plant［J］.Marine Environmental Science，1998，17（2）：54－57.韓康，張存智.三亞電廠溫排水?dāng)?shù)值模擬［J］.海洋環(huán)境科學(xué)，1998，17（2）：54－57.

［2］ ZHANG X B，JIN K，LIN M S.2DNumerical simulation for warm water drainage in tide river reach［J］.Journal of Yangtze River Scientific Research Institute，2006，23（3）：13－16.張細(xì)兵，金琨，林木松.潮流河段溫排水影響的平面二維數(shù)值模擬［J］.長江科學(xué)院院報(bào)，2006，23（3）：13－16.

［3］ CAO Y，ZHU J Z.Numerical simulation of 3Dcooling water based on FVCOM［J］.Chinese Journal of Hydrodynamic，2009，24（4）：432－439.曹穎，朱軍政，基于FVCOM 模式的溫排水三維數(shù)值模擬研究［J］.水動(dòng)力學(xué)研究與進(jìn)展，2009，24（4）：432－439.

［4］ JIANG S，DUAN M L，WANG S G.Advances in seawater thermal diffusion and new issues analysis［J］.Energy Environmental Protection，2006，20（5）：5－9.蔣爽，端木琳，王樹剛.海水熱擴(kuò)散研究進(jìn)展與新問題分析［J］.能源環(huán)境保護(hù)，2006，20（5）：5－9.

［5］ FANG G H.Quasi－h(huán)armonic constituent method for analysis and prediction of tides III［C］.Beijing：Science Press，1981，22：19－40.方國洪.潮汐分析和預(yù)報(bào)的準(zhǔn)調(diào)和分潮方法III［C］.北京：科學(xué)出版社，1981，22：19－40.

［6］ HE S Q，HUAI W X，ZENG Y H.Numerical simulation on 2－D thermal discharge into continuous curved channel［J］.Resources and Environment in the Yangtze Basin，2007，16（1）：17－20.何書琴，槐文信，曾玉紅.連續(xù)彎道平面二維溫排水的數(shù)值模擬［J］.長江流域資源與環(huán)境，2007，16（1）：17－20.

［7］ CHANG X J.Numerical simulation of cooling water in Jiaozhou Bay［D］.Qingdao：Ocean University of China，2011.常小軍.膠州灣溫排水?dāng)?shù)值模擬［D］.青島：中國海洋大學(xué)，2011.