渤海灣溫帶風(fēng)暴潮數(shù)值預(yù)報(bào)模型

李大鳴，徐亞男，白 玲，解以揚(yáng)，吳丹朱，何乃光

(1. 天津大學(xué)建筑工程學(xué)院暨港口與海洋工程教育部重點(diǎn)實(shí)驗(yàn)室，天津 300072；2. 天津市氣象科學(xué)研究所，天津 300074)

渤海灣海域是中國(guó)渤海三大海灣之一，位于渤海西部，海底地形大致自南向北，自岸向海傾斜，屬風(fēng)暴潮災(zāi)害的多發(fā)區(qū)和嚴(yán)重區(qū)．其災(zāi)害多發(fā)生在盛夏臺(tái)風(fēng)活動(dòng)季節(jié)和春秋過(guò)渡季節(jié)[1]．渤海灣中有豐富的石油儲(chǔ)藏(大港油田、冀東南堡油田)，灣內(nèi)有天津新港，每年因?yàn)轱L(fēng)暴潮災(zāi)害經(jīng)濟(jì)損失達(dá)億元以上，因此建立適用于渤海灣海域的風(fēng)暴潮數(shù)值預(yù)報(bào)系統(tǒng)，對(duì)減少港區(qū)經(jīng)濟(jì)損失，降低由溢油擴(kuò)散及污染源引起的環(huán)境污染有重要的經(jīng)濟(jì)效益和環(huán)境效益．研究表明在溫帶風(fēng)暴潮過(guò)程中，氣壓場(chǎng)與波浪場(chǎng)較風(fēng)應(yīng)力場(chǎng)在渤海灣海域的影響有限[2-3]，可以在此區(qū)域忽略；但由于渤海灣為 3面環(huán)陸的半封閉性海灣，岸線較長(zhǎng)且地形復(fù)雜，溫帶風(fēng)暴潮數(shù)值預(yù)報(bào)模型的建立仍具有一定難度．筆者所建立的嵌套模型采用多分潮調(diào)和分析方法解決了海域開(kāi)邊界處理問(wèn)題，提高了渤海灣海域模型空間網(wǎng)格的分辨率，并應(yīng)用顯隱交替的有限差分格式(alternating direction implicit，ADI)對(duì)風(fēng)暴潮控制方程離散求解，采用局部深槽、縮小水域的計(jì)算模式處理了動(dòng)態(tài)淺水岸邊界，提高了模型計(jì)算的穩(wěn)定性．

國(guó)內(nèi)外學(xué)者已經(jīng)發(fā)展了多種風(fēng)暴潮預(yù)報(bào)模式[4]，早期美國(guó)的Jelesnianski等[5]發(fā)展了 SPLASH(special program to list amplitudes of surges from hurricane)模式，后在該模式的基礎(chǔ)上美國(guó)發(fā)展了新的 SLOSH(sea，lake and overland surges from hurricanes)模式來(lái)預(yù)報(bào)海上、陸上、湖上的臺(tái)風(fēng)風(fēng)暴潮；Blumberg[6]發(fā)展了POM(Princeton ocean model)模式，該模式可以計(jì)算小尺度河川的水理運(yùn)動(dòng)，也可以模擬大尺度海洋、海岸的水位與流場(chǎng)的變化．在國(guó)內(nèi)，李艷蕓等[7]在風(fēng)暴潮預(yù)報(bào)模式于渤海海域中的應(yīng)用研究中，采用COHERENS(coupled hydrodynamic-ecological model for regional and shelf seas)三維多功能大陸架水動(dòng)力數(shù)學(xué)模型模擬了熱帶風(fēng)暴下渤海的增水過(guò)程；于福江等[8]建立了球坐標(biāo)系下的溫帶風(fēng)暴潮模式，對(duì)渤海灣一次特大風(fēng)暴潮過(guò)程進(jìn)行了數(shù)值模擬．但連續(xù)、快速的溫帶風(fēng)暴潮預(yù)報(bào)模型的研究，特別是對(duì)多個(gè)例風(fēng)暴潮過(guò)程預(yù)報(bào)的數(shù)學(xué)模型研究還需要做大量的工作．筆者自主創(chuàng)建溫帶風(fēng)暴潮預(yù)報(bào)模型，與氣象部門合作，已應(yīng)用于對(duì)渤海灣海域長(zhǎng)系列天文潮、風(fēng)暴潮、風(fēng)暴潮增水的長(zhǎng)期、連續(xù)、快速預(yù)報(bào)計(jì)算，在大模型網(wǎng)格 10 km×10 km，小模型網(wǎng)格 1 km×1 km，模擬3,d的72 h風(fēng)暴潮過(guò)程，在現(xiàn)行一般配置的個(gè)人計(jì)算機(jī)上只需運(yùn)行5 min．筆者選擇了近期預(yù)報(bào)的3次渤海灣風(fēng)暴潮過(guò)程，模擬結(jié)果與實(shí)際觀測(cè)進(jìn)行比較，表明該模型對(duì)渤海灣風(fēng)暴潮過(guò)程預(yù)報(bào)具有一定的實(shí)際應(yīng)用價(jià)值．

1 數(shù)值模型建立的理論基礎(chǔ)

1.1 水動(dòng)力數(shù)學(xué)模型的控制方程

溫帶風(fēng)暴潮是由大氣強(qiáng)迫力(風(fēng)場(chǎng)、氣壓)作用于海面造成的海水水位與潮流的劇烈變化．因?yàn)椴澈呈瞧骄钪挥?8 m的3面環(huán)陸的淺海，地形變化復(fù)雜，氣壓對(duì)風(fēng)暴潮過(guò)程的影響較小[9]，只考慮風(fēng)應(yīng)力場(chǎng)和潮汐作用．模型采用直角坐標(biāo)系，假定沿水深方向的動(dòng)水壓強(qiáng)分布符合靜水壓強(qiáng)分布，將三維流動(dòng)的基本方程沿水深平均積分即可得到沿水深平均的平面二維流動(dòng)的基本方程[10]．以風(fēng)應(yīng)力、水流與海底的摩擦應(yīng)力[11]為主要影響因子，水動(dòng)力數(shù)學(xué)模型控制方程為

式中：ξ為增水位；h為平均水深；H為全水深，H＝ξ+h ；u、v分別為x、y方向上流速分量；τx,s、τy,s分別為x和y方向的海面風(fēng)應(yīng)力；f為柯氏系數(shù)．

模型岸邊界條件為：vn= 0 (n為邊界法線方向)．

1.2 方程離散求解

采用 ADI方法對(duì)方程進(jìn)行離散求解，差分的交錯(cuò)網(wǎng)格為正方形網(wǎng)格，網(wǎng)格線分別平行于x軸和y軸，間距為 Δ x = Δy =Δs ．

1.3 風(fēng)場(chǎng)的預(yù)報(bào)模式

模型驅(qū)動(dòng)力風(fēng)場(chǎng)采用 MM,5模式[12-14]計(jì)算成果．MM,5以 NECP資料中的 0.5×0.5格距的 GFS資料做背景場(chǎng)，中心點(diǎn)位置選在 41.1,N，118.2,E處．模型網(wǎng)格分辨率9,km，利用探空和地面觀測(cè)資料對(duì)背景場(chǎng)進(jìn)行修正，從而得到 10,m 風(fēng)場(chǎng)作為風(fēng)暴潮的主要驅(qū)動(dòng)力；模型能夠預(yù)報(bào)連續(xù) 72 h的逐時(shí)風(fēng)場(chǎng)．將計(jì)算得到的經(jīng)緯度坐標(biāo)下的風(fēng)場(chǎng)數(shù)據(jù)插值到直角坐標(biāo)系下的模型模擬區(qū)域的格點(diǎn)內(nèi)．

模型中風(fēng)應(yīng)力計(jì)算則采用應(yīng)用較廣泛的公式

式中：W為海面10,m風(fēng)速；aρ為空氣密度，取為1.226,kg/m3；DC 按經(jīng)驗(yàn)取為2.6×10-3．

1.4 調(diào)和分析方法

模型計(jì)算區(qū)域是整個(gè)渤海與部分黃海，開(kāi)邊界為青島港(35.43,N，119.58,E)到韓國(guó)西岸港口HAMPYEONG MAN(35.15,N，126.35,E)的連線，對(duì)此邊界進(jìn)行了多分潮調(diào)和分析計(jì)算．根據(jù)驗(yàn)潮站 1 a的潮位實(shí)測(cè)資料，采用Sa、SSa、Mm、Mf等 35個(gè)分潮[14]，見(jiàn)表 1(部分分潮)．最后 1列為調(diào)和分析計(jì)算成果．

表1 調(diào)和分析中采用的部分分潮參數(shù)Tab.1 Chon-tide schedule in harmonic analysis

各分潮的調(diào)和常數(shù)即振幅與遲角的展開(kāi)公式[15]為

式中：0a為基準(zhǔn)面的平均海平面高度；0V u+ 為分潮的天文初相角；jH、jg為分潮振幅與遲角；j、m為分潮序號(hào)與總數(shù)；jσ為分潮角頻率；jf為平均振幅的訂正系數(shù)；t為時(shí)間．

1.5 動(dòng)態(tài)淺水岸邊界處理

ADI差分計(jì)算模式要求整個(gè)計(jì)算域應(yīng)保持在水深以下，對(duì)淺水岸邊界的露灘、淹沒(méi)變化應(yīng)是連續(xù)、穩(wěn)定過(guò)程．本文采用局部深槽、縮小水域的活動(dòng)邊界處理方法，當(dāng)全水深接近0.1 m時(shí)，在淺水網(wǎng)格區(qū)全水深保持為 H10= 0 .1m ，保持流量、流速不變，流量為

變化后的水域?qū)挾萐B為

2 風(fēng)暴潮數(shù)值預(yù)報(bào)模型的應(yīng)用及結(jié)果分析

圖1 模型計(jì)算區(qū)域范圍及嵌套模式示意Fig.1 Computed domain and setup of nesting domain

模型采用嵌套網(wǎng)格模式如圖1所示，網(wǎng)格剖分為正方形網(wǎng)格[16]，其中第 1套網(wǎng)格構(gòu)成的大模型空間步長(zhǎng)為10 km，模型范圍涵蓋整個(gè)渤海及部分黃海海域，模型計(jì)算域跨度為 117°38'47"E至 126°32'38"E，35°18'5"N至40°50'43"N，計(jì)算結(jié)果包括部分黃海及渤海的潮流與潮位過(guò)程，為第2套網(wǎng)格提供海洋開(kāi)邊界條件；第 2套網(wǎng)格構(gòu)成的小模型空間步長(zhǎng)為 1,km覆蓋整個(gè)渤海灣，岸邊界網(wǎng)格更加精細(xì)化，為實(shí)現(xiàn)岸邊界的淺水動(dòng)態(tài)處理提供了條件，模型網(wǎng)格見(jiàn)圖2．大模型網(wǎng)格總數(shù)為 78×61，計(jì)算時(shí)間步長(zhǎng) 60,s；小模型網(wǎng)格總數(shù) 121×161，計(jì)算時(shí)間步長(zhǎng) 10,s，其中圖2(b)的標(biāo)記點(diǎn)為渤海灣內(nèi)的主要港口和地區(qū)，包括塘沽、曹妃甸、黃驊港等，模型以塘沽為驗(yàn)證點(diǎn)，計(jì)算了3次典型的風(fēng)暴潮過(guò)程，驗(yàn)證結(jié)果表明，所建風(fēng)暴潮數(shù)值預(yù)報(bào)模型可應(yīng)用于渤海灣地區(qū)的風(fēng)暴潮預(yù)警預(yù)報(bào)系統(tǒng).

圖2 黃渤海區(qū)域與渤海灣網(wǎng)格劃分示意Fig.2 Mesh grid of Yellow-Bohai Sea and Bohai Bay

2.1 調(diào)和分析結(jié)果驗(yàn)證

由 2002年青島港實(shí)測(cè)潮位資料計(jì)算出 35個(gè)分潮的調(diào)和常數(shù)，以此分別計(jì)算選定的3個(gè)不同時(shí)段的青島港潮位過(guò)程，將計(jì)算結(jié)果與實(shí)測(cè)資料對(duì)比驗(yàn)證如圖 3所示．對(duì)比曲線表明，應(yīng)用調(diào)和分析方法具有可行性．韓國(guó)西岸港口HAMPYEONG MAN的潮位資料采用潮位軟件提供，模型以青島港潮位過(guò)程為主，考慮韓國(guó)西岸港口潮位過(guò)程修正，確定大模型計(jì)算域的開(kāi)邊界條件．

圖3 青島港潮位調(diào)和分析驗(yàn)證Fig.3 Comparison and analysis of tidal level in Qingdao Harbor

2.2 風(fēng)場(chǎng)模式計(jì)算結(jié)果分析

圖4 2007-08-12 14:00黃渤海海域與渤海灣風(fēng)速分布Fig.4 Distribution of wind speed of Yellow-Bohai Sea and Bohai Bay at 2 pm Aug 12th, 2007

渤海灣為3面環(huán)陸的內(nèi)海，風(fēng)暴潮過(guò)程受風(fēng)場(chǎng)影響最為顯著，尤其受到東風(fēng)與東北風(fēng)的作用后海面水位抬高顯著．模型預(yù)報(bào)的渤海灣 3次典型風(fēng)暴潮過(guò)程分別為070812次、090415次和090509次，風(fēng)場(chǎng)在渤海灣的分布主要呈東北偏東風(fēng)形勢(shì)．如圖 4所示，渤海灣 070812次風(fēng)暴潮受自身海域氣旋影響較小，而受黃海海域氣旋影響顯著，從而導(dǎo)致潮位超過(guò)警戒水位；影響 090415次風(fēng)暴潮過(guò)程的風(fēng)場(chǎng)是東北大風(fēng)，風(fēng)速最大達(dá)到 21.29,m/s；090509次風(fēng)暴潮過(guò)程直接受當(dāng)日的7級(jí)偏東大風(fēng)影響，導(dǎo)致塘沽驗(yàn)潮站的風(fēng)暴潮潮位超過(guò)警戒水位．警戒水位預(yù)設(shè)為4.7,m．

2.3 數(shù)值模擬結(jié)果

圖5 風(fēng)暴潮潮位過(guò)程曲線Fig.5 Tidal level process curves of storm surge

將已建立的風(fēng)暴潮數(shù)值預(yù)報(bào)模型應(yīng)用到渤海灣，調(diào)試驗(yàn)證了渤海灣 070812次風(fēng)暴潮、同時(shí)對(duì)風(fēng)暴潮090415次、090509次進(jìn)行了預(yù)報(bào)檢驗(yàn)．圖 5與圖 6給出了計(jì)算潮位與實(shí)測(cè)潮位過(guò)程和計(jì)算增水值與實(shí)測(cè)增水值的曲線對(duì)比驗(yàn)證．從潮位的變化趨勢(shì)看，070812次風(fēng)暴潮計(jì)算潮位與實(shí)測(cè)潮位擬合一般相關(guān)系數(shù)為 0.93，同時(shí) 090415次、090509次預(yù)報(bào)的風(fēng)暴潮潮位與實(shí)測(cè)潮位符合程度較高相關(guān)系數(shù)分別為0.96和 0.97(見(jiàn)表 2)；從增水過(guò)程變化趨勢(shì)看，對(duì) 3次風(fēng)暴增水的峰值能夠較好地模擬出來(lái)．070812次風(fēng)暴潮實(shí)測(cè)最高潮位 4.83,m，增水 0.86,m，而模型預(yù)報(bào)結(jié)果顯示塘沽最高潮位達(dá) 4.96,m，增水 1.38,m，與實(shí)測(cè)值較為接近；090415次風(fēng)暴潮受冷空氣與氣旋影響，塘沽驗(yàn)潮站測(cè)得該日6∶00發(fā)生4.94,m大潮，增水為 1.84,m，模型提前 2日計(jì)算得到 4月 15日6∶00風(fēng)暴潮最高潮位為 5.08,m，計(jì)算增水達(dá)到1.94,m；090509次風(fēng)暴潮受南下冷空氣影響，風(fēng)力最高達(dá) 7級(jí)，實(shí)際觀測(cè)塘沽站最高潮位 4.95,m，增水1.27,m，模型提前2日計(jì)算結(jié)果及誤差分析如表2所示，模型模擬潮位峰值的平均相對(duì)誤差為 0.024，增水峰值的平均相對(duì)誤差也在允許的范圍內(nèi)．上述結(jié)果表明，建立的風(fēng)暴潮數(shù)值預(yù)報(bào)模型能夠?qū)Σ澈澈^(qū)風(fēng)暴潮進(jìn)行較為準(zhǔn)確的預(yù)報(bào)，可以用來(lái)對(duì)該地區(qū)風(fēng)暴潮進(jìn)行預(yù)報(bào)．

但是，該模型在模擬增水峰值出現(xiàn)的時(shí)間上與實(shí)測(cè)值存在誤差，其原因是多方面的：在對(duì)天文潮進(jìn)行非線性模擬時(shí)，由于缺失水邊界實(shí)測(cè)資料采用調(diào)和分析方法，雖然考慮了 35個(gè)分潮的作用，但以青島港的潮位值來(lái)代替整個(gè)水邊界的潮位值，可能導(dǎo)致峰值與相位的偏差；在本文所建立的風(fēng)暴潮預(yù)報(bào)模型中，近岸地形采用的是20世紀(jì)90年代資料，這與實(shí)際近岸地形分布可能存在一定偏差，從而造成計(jì)算模擬誤差，但是隨著預(yù)報(bào)經(jīng)驗(yàn)的增加和資料的不斷收集完善，可望減小預(yù)報(bào)誤差．

圖6 風(fēng)暴潮增水驗(yàn)證曲線Fig.6 Tidal rising process curves of storm surge

表2 塘沽站風(fēng)暴潮潮位及增水值誤差分析Tab.2 Error analysis on tidal level and tidal rising of storm surge at Tanggu station

3 結(jié) 論

(1)渤海灣溫帶風(fēng)暴潮數(shù)值預(yù)報(bào)模型，采用大小嵌套模型模式，以多分潮調(diào)和分析提供海洋邊界條件，以MM5風(fēng)場(chǎng)計(jì)算成果形成時(shí)空分布風(fēng)應(yīng)力場(chǎng)過(guò)程，模型主體網(wǎng)格采用 ADI差分格式進(jìn)行顯隱交替計(jì)算，實(shí)現(xiàn)了連續(xù)、快速達(dá)到一定精度要求的風(fēng)暴潮預(yù)報(bào)運(yùn)行模式，為渤海灣溫帶風(fēng)暴潮數(shù)值預(yù)報(bào)提供了可靠的研究手段．

(2)渤海灣屬淺海類型海灣，岸灘地形變化復(fù)雜，本文中提出淺水變動(dòng)岸邊界的計(jì)算模式，在 ADI差分格式連續(xù)計(jì)算中，采用局部深槽、縮小水域的活動(dòng)邊界處理方法，增加了模型計(jì)算的穩(wěn)定性，提高了模型模擬風(fēng)暴潮在岸灘附近增水計(jì)算的能力．

(3)本研究建立的海洋潮波動(dòng)力和風(fēng)應(yīng)力場(chǎng)聯(lián)合作用的溫帶風(fēng)暴潮數(shù)值預(yù)報(bào)模型，預(yù)報(bào)并驗(yàn)證了渤海灣海域近期 3次風(fēng)暴潮過(guò)程，預(yù)報(bào)潮位過(guò)程、增水過(guò)程與實(shí)測(cè)值進(jìn)行比較，潮位過(guò)程吻合較好；增水過(guò)程在峰值處比較接近，平均相對(duì)誤差較小，其他各處趨勢(shì)基本一致，表明該模型在渤海灣海域溫帶風(fēng)暴潮預(yù)報(bào)模擬中具有一定的應(yīng)用價(jià)值，可以用來(lái)預(yù)報(bào)該海區(qū)的風(fēng)暴潮過(guò)程．

［1］ 李大鳴，徐亞男，宋雙霞，等. 波浪輻射應(yīng)力在渤海灣海域?qū)︼L(fēng)暴潮影響的研究[J]. 水動(dòng)力學(xué)研究與進(jìn)展A輯，2010，25(3)：374-382.Li Daming，Xu Yanan，Song Shuangxia，et al. The study on the effect of wave-radiation stress on storm surge in Bohai Bay[J]. Chinese Journal of Hydrodynamics，Ser A，2010，25(3)：374-382(in Chinese).

［2］ Andrey P，Eppel D P，Kapitza H. Interaction of waves，currents and tides，and wave-energy impact on the beach area of Sylt Island[J]. Ocean Dynamics，2009，59(3)：451-461.

［3］ 林 祥，尹寶樹(shù)，侯一筠，等. 輻射應(yīng)力在黃河三角洲近岸波浪和潮汐風(fēng)暴潮相互作用中的影響[J]. 海洋與湖沼，2002，33(6)：615-621.Lin Xiang，Yin Baoshu，Hou Yijun，et al. Effects of radiation stress in the interaction of coupled wave-tidesurge in the coastal area of Huanghe Delta[J]. Oceanologia Etlimnologia Sinica，2002，33(6)：615-621(in Chinese).

［4］ Kim K O，Yamashita T. Storm surge simulation using wind-wave-surge coupling model[J]. Journal of Oceanography，2008，64：621-630.

［5］ Jelesnianski C P，Chen J，Shaffer W A. SLOSH：Sea，Lake，and Overland Surges from Hurricanes[R].National Weather Service，USA，1992.

［6］ Blumberg A F，Mellor G L. A description of a threedimensional coastal ocean circulation model[C]//Three-Dimension Coastal Ocean Models. Washington，USA，1987：1-16.

［7］ 李艷蕓，李紹武. 風(fēng)暴潮預(yù)報(bào)模式在渤海海域中的應(yīng)用研究[J]. 海洋技術(shù)，2006，25(1)：101-106.Li Yanyun，Li Shaowu. Application research of a storm surge prediction model in Bohai Sea[J]. Ocean Technology，2006，25(1)：101-106(in Chinese).

［8］ 于福江，王喜年，宋 珊，等. 渤?！?216”特大風(fēng)暴潮過(guò)程的數(shù)值模擬[J]. 海洋預(yù)報(bào)，2000，17(4)：9-15.Yu Fujiang，Wang Xinian，Song Shan，et al. The numerical simulation of storm surge in Bohai Sea caused by tropical storm POLLY[J]. Marine Forecast，2000，17(4)：9-15(in Chinese).

［9］ 陳士蔭，顧家龍. 海岸動(dòng)力學(xué)[M]. 北京：人民交通出版社，1988.Chen Shiyin，Gu Jialong. Coastal Dynamics[M].Beijing：People’s Communication Press，1988(in Chinese).

［10］ Jones J E，Davies A M. Storm surge computations for the west coast of Britain using a finite element model(TELEMAC)[J]. Ocean Dynamics，2008，58(5/6)：337-363.

［11］ Lee D. Bottom shear stress under wave-current interaction[J]. Journal of Hydrodynamics，Ser B，2008，20(1)：88-95.

［12］ Lee S M，Princevac M，Mitsutomi S，et al. MM5 simulations for air quality modeling：An application to a coastal area with complex terrain[J]. Atmospheric Environment，2009，43(2)：447-457.

［13］ Heo K Y，Lee J W，Ha K J，et al. Simulation of atmospheric states for a storm surge on the west coast of Korea：Model comparison between MM5，WRF and COAMPS[J]. Nat Harzards，2009，51(1)：151-162.

［14］ Ivanov S，Palamarchuk J，Pyshniak D. Upscale feedbacks through microphysics fields at nesting domains of the MM5 model[J]. Atmospheric Research，2009，94(1)：726-735.

［15］ 陳宗鏞. 潮汐學(xué)[M]. 北京：科學(xué)出版社，1980.Chen Zongyong. Tidal Science[M]. Beijing：Science Press，1980(in Chinese).

［16］ 于福江，張占海. 一個(gè)東海嵌套網(wǎng)格臺(tái)風(fēng)暴潮數(shù)值預(yù)報(bào)模式的研制與應(yīng)用[J]. 海洋學(xué)報(bào)，2002，24(4)：23-33.Yu Fujiang，Zhang Zhanhai. Implementation and application of a nested numerical typhoon storm surge forecast model in the East China Sea[J]. Acta Oceanologica Sinica，2002，24(4)：23-33(in Chinese).