基于WRF-波浪風(fēng)暴潮耦合模型的深圳市極端風(fēng)暴潮安全分析

林澤輝

(惠州市惠陽(yáng)區(qū)永良堤圍管理所，廣東 惠州 516007)

我國(guó)東南沿海地區(qū)氣象水文條件復(fù)雜，水文方面，海水水位受天體引潮力作用而周期性漲落。氣象方面，西太平洋熱帶地區(qū)氣壓受到擾動(dòng)作用成渦，并不斷發(fā)展移動(dòng)，這些氣旋登陸后會(huì)對(duì)人類社會(huì)造成巨大的破壞。另一方面，臺(tái)風(fēng)所造成的大氣壓強(qiáng)場(chǎng)變化也會(huì)導(dǎo)致海面異常起伏，產(chǎn)生風(fēng)暴潮，當(dāng)風(fēng)暴潮恰好疊加上天文大潮，海平面會(huì)巨幅上升，這對(duì)海岸防護(hù)、河口擋潮等水利工程帶來(lái)了巨大挑戰(zhàn)。因此研究臺(tái)風(fēng)過(guò)程前后的波浪以及增水具有重要意義。文章利用…模型對(duì)深圳地區(qū)水域在臺(tái)風(fēng)作用下的水動(dòng)力特征進(jìn)行了數(shù)值模擬研究。

1 研究方法

1.1 區(qū)域概況

深圳市海岸線長(zhǎng)約250km，東起大亞灣，西至珠江口，增水>50cm的風(fēng)暴潮過(guò)程年平均出現(xiàn)2.1次。沿岸分布有機(jī)場(chǎng)、重點(diǎn)港口、大型住宅區(qū)、核電站等大型企業(yè)，人口密度極高。

深圳市年平均風(fēng)速2.6m/s，冬季各月風(fēng)速較大為3.0m/s左右，夏季各月較小約2.0m/s。大風(fēng)日數(shù)(風(fēng)速≥17m/s，即風(fēng)力相當(dāng)于8級(jí)或8級(jí)以上)年平均為7.3d。全年各月均出現(xiàn)大風(fēng)，其中以7-9月居多，共占全年大風(fēng)日數(shù)的44%，其次為6月和10月，由于受臺(tái)風(fēng)或熱帶低壓影響，常出現(xiàn)大風(fēng)，因此，大風(fēng)出現(xiàn)多的月份，也是臺(tái)風(fēng)季節(jié)。

1.2 數(shù)據(jù)資料

影響風(fēng)暴增水量級(jí)的因素有很多，包括臺(tái)風(fēng)氣壓場(chǎng)分布、臺(tái)風(fēng)風(fēng)場(chǎng)的方向與大小以及近岸地形變化因素。在從我國(guó)東南沿海登陸的臺(tái)風(fēng)里，影響到陽(yáng)江至汕尾海域的臺(tái)風(fēng)均可引起深圳海域的風(fēng)暴潮增水。目前在深圳海域，國(guó)家海洋局共設(shè)有兩個(gè)長(zhǎng)期驗(yàn)潮站，分別是赤灣和鹽田站，緊鄰深圳大亞灣岸線的有位于惠州市的惠州海洋站，各站位置見(jiàn)圖1。收集上述站點(diǎn)自建站以來(lái)臺(tái)風(fēng)期間的潮位資料統(tǒng)計(jì)歷史風(fēng)暴潮增水情況。廣東省水文局1964年在深圳市赤灣港內(nèi)建設(shè)有長(zhǎng)期水文觀測(cè)站，1986年國(guó)家海洋局南海分局在相距不到1km處建設(shè)了赤灣海洋站，并于1986年開(kāi)始投入業(yè)務(wù)化觀測(cè)。鹽田海洋站位于大鵬灣鹽田港海域，由國(guó)家海洋局南海分局建設(shè)和管理，2003年投入使用?；葜莺Ｑ笳疚挥诖髞啚碁稠斘鱾?cè)荃灣港區(qū)，由國(guó)家海洋局南海分局建設(shè)和管理，2006年建成投入使用，他們的分布位置如圖1所示。

圖1 資料站位位置示意圖

1.3 研究方法

早期人們主要通過(guò)經(jīng)驗(yàn)總結(jié)歸納的方法研究風(fēng)暴潮，觀察臺(tái)風(fēng)天時(shí)海塘外潮位變化研究風(fēng)暴潮的成因。20世紀(jì)中期，隨著計(jì)算機(jī)技術(shù)的發(fā)展，人們開(kāi)始利用數(shù)值模擬研究風(fēng)暴潮。1956年W·Hansen首次利用計(jì)算機(jī)對(duì)歐洲北海的最大風(fēng)暴潮增水進(jìn)行了研究。隨著數(shù)值模擬技術(shù)的進(jìn)步，Jelesnianski建立了研究風(fēng)暴潮的SPLASH模式以及基于二維水動(dòng)力的SLOSH模式，在不斷地摸索研究后，該模式的預(yù)報(bào)精度不斷提高。近年來(lái)，不斷發(fā)展的模型計(jì)算實(shí)現(xiàn)了對(duì)復(fù)雜地形影響的考慮，提高了計(jì)算的準(zhǔn)確性[1]。其次，有關(guān)風(fēng)暴增水與天文潮的非線性擬合機(jī)制也變得更加清晰，這也提高了風(fēng)暴潮數(shù)學(xué)模型的準(zhǔn)確性。我國(guó)自20世紀(jì)60年代也開(kāi)展了關(guān)于風(fēng)暴潮的相關(guān)研究。馮士筰所著的《風(fēng)暴潮導(dǎo)論》系統(tǒng)闡述了風(fēng)暴潮的概念，形成機(jī)制，并介紹了早期的風(fēng)暴潮數(shù)值預(yù)報(bào)方法。趙長(zhǎng)進(jìn)等利用ADCIRC潮流與SWAN波浪耦合模式對(duì)長(zhǎng)江口附近海域的風(fēng)暴潮過(guò)程進(jìn)行了模擬。黃本勝等用SELFE模型構(gòu)建了南?！榻涌陔p重嵌套模型，研究了“山竹”臺(tái)風(fēng)過(guò)程中珠江河口的風(fēng)暴增水與波浪增水分布。

在波浪模擬方面，目前工程界使用較為廣泛的數(shù)學(xué)模型有兩類，一類是波浪流體動(dòng)力學(xué)模型，另一類是動(dòng)譜密度方程。緩坡方程模型與Boussinesq方程模型是主要的波浪流體力學(xué)模型。經(jīng)典的Boussinesq僅僅適用于淺水，并不適用于水深變化大的水域和不規(guī)則波的情況。因此在工程中適用范圍不廣。動(dòng)譜密度方程的基礎(chǔ)是能量平衡理論，在模擬過(guò)程中可提高海岸波浪場(chǎng)的動(dòng)力作用過(guò)程的計(jì)算效率。該方程以動(dòng)譜密度來(lái)描述波浪變形，考慮了對(duì)波-波非線性相互作用、波流相互作用、波浪破碎和底部摩擦引起的能量耗散等過(guò)程，通過(guò)將源項(xiàng)添加到方程中，可以全面考慮動(dòng)力作用的復(fù)雜過(guò)程，使其能完成大尺度區(qū)域的波浪產(chǎn)生、成長(zhǎng)、傳播和變形的計(jì)算過(guò)程。其中對(duì)物理源項(xiàng)的處理和研究是波作用量守恒的關(guān)鍵，波作用量守恒方程發(fā)展到現(xiàn)在已經(jīng)到了第三代，逐漸開(kāi)發(fā)出WAVEWATCH、SWAN、MIKE21SW、TOMAWAC等模型。

1.4 模型設(shè)置

1.4.1 模型介紹

文章采用的是DELFT3D-FLOW與SWAN的耦合模型，DELFT3D是一款由荷蘭DELTRA開(kāi)發(fā)的開(kāi)源水動(dòng)力模擬軟件，在風(fēng)暴潮模擬方面已經(jīng)廣泛應(yīng)用。

Delft3D-Flow是基于淺水方程和Boussinesq假設(shè)，采用有限差分法求解斜壓Navier-Stokes方程的模型，模型控制方程組分別為：

質(zhì)量守恒方程：

(1)

式中：ζ為水位(總水深)；d為到參考平面的當(dāng)前水深(凈水深)；U、V-x、y為方向上的速度分量；Q為單位面積上質(zhì)量強(qiáng)度。

動(dòng)量守恒方程：

(2)

(3)

式中：f為科氏力參數(shù)；g為重力加速度；νh為動(dòng)態(tài)水平渦流黏度；τsx,τsy為x、y方向上的海平面風(fēng)壓分量；τbx,τby為x、y方向上的底部的剪應(yīng)力分量；ρ0為參考密度；ρ'為不規(guī)則密度。

SWAN模型的控制方程為：

(4)

式中：Cx，Cy，Cσ，Cθ為二維動(dòng)譜密度在x、y、σ、θ為空間上的傳播速度；σ為相對(duì)頻率；θ為波向；S為能量源項(xiàng)。

1.4.2 模型設(shè)置

1)模型范圍：

文章所采用的模型囊括整個(gè)珠江三角洲區(qū)域，見(jiàn)圖2。

圖2 模型區(qū)域示意圖

2)地形水深：

模擬范圍內(nèi)的地形水深采用的是中華人民共和國(guó)海事局2016年繪制的海圖插值數(shù)據(jù)。將地形水深數(shù)據(jù)轉(zhuǎn)化到統(tǒng)一基面(平均海平面)后，用IDW插值方法插值到網(wǎng)格點(diǎn)上，得到模型的水深地形。

3)開(kāi)邊界設(shè)置：

開(kāi)邊界采用天文潮水位，主要考慮M2、N2、S2、K2、K1、O1、P1、Q1等八個(gè)天文分潮以及MS4、M4、M6等淺水分潮共11個(gè)分潮的作用。分潮的調(diào)和分析常數(shù)由2016-2019年的潮位實(shí)測(cè)數(shù)據(jù)經(jīng)MATLAB中的T_TIDE工具包調(diào)和分析所得。模型上部采用流量邊界，輸入數(shù)據(jù)為梧州、石角、博羅三個(gè)水文站的2016-2019洪季平均流量。

4)底摩擦設(shè)置：

模型根據(jù)水深范圍設(shè)置不同大小的糙率，水深越大糙率越小。

1.4.3 風(fēng)場(chǎng)輸入

文章的風(fēng)場(chǎng)采用美國(guó)國(guó)家大氣研究中心(NCAR)等機(jī)構(gòu)共同開(kāi)發(fā)的中尺度大氣預(yù)報(bào)系統(tǒng)。WRF模式適用于臺(tái)風(fēng)模擬，可以通過(guò)嵌套網(wǎng)格運(yùn)算來(lái)達(dá)到高精度模擬的目的。本次采用了WSM6微物理方案、KainFritsch卷積云方案、Monin-Obukhov表面層方案、Noah陸面層方案、MYJ大氣邊界層方案、RRTMG輻射方案模擬臺(tái)風(fēng)。

2 模型驗(yàn)證

2.1 風(fēng)場(chǎng)驗(yàn)證

首先利用實(shí)測(cè)資料對(duì)所建立的風(fēng)暴潮-波浪耦合模型進(jìn)行率定驗(yàn)證。選取1713號(hào)臺(tái)風(fēng)(2017年第13號(hào)臺(tái)風(fēng)—天鴿，其移行路徑如圖3所示)作為模型邊界條件，為了證明WRF模式的準(zhǔn)確性其在風(fēng)暴潮模擬方面的優(yōu)勢(shì)，文章也用其他論文常用的經(jīng)驗(yàn)風(fēng)場(chǎng)模型做了驗(yàn)證對(duì)比。經(jīng)驗(yàn)風(fēng)場(chǎng)采用的是宮崎正衛(wèi)模型作為移動(dòng)風(fēng)場(chǎng)，HOLLAND模型作為梯度風(fēng)場(chǎng)，他們的公式分別為：

圖3 “天鴿”臺(tái)風(fēng)移行路徑

(5)

(5)

式中：B為HOLLAND模型氣壓參數(shù)；f為科氏力參數(shù)，與模擬區(qū)域緯度有關(guān)；R為最大風(fēng)速半徑，與中心氣壓有關(guān)。

圖4為WRF模型模擬風(fēng)場(chǎng)值以及經(jīng)驗(yàn)風(fēng)場(chǎng)模擬值與實(shí)測(cè)風(fēng)場(chǎng)的對(duì)比，可以看出，WRF模型模擬的風(fēng)場(chǎng)在峰值以及過(guò)程上都比經(jīng)驗(yàn)風(fēng)場(chǎng)模型更貼近實(shí)際，因此利用WRF模型生成的臺(tái)風(fēng)風(fēng)場(chǎng)作為風(fēng)暴潮-波浪耦合模型的邊界條件更為合理。

圖4 WRF風(fēng)場(chǎng)及模型風(fēng)場(chǎng)與實(shí)測(cè)風(fēng)場(chǎng)對(duì)比

2.2 風(fēng)暴潮及波浪率定驗(yàn)證

以1713號(hào)臺(tái)風(fēng)(2017年第13號(hào)臺(tái)風(fēng))為邊界條件分別驅(qū)動(dòng)珠江口風(fēng)暴潮模型以及珠江口波浪-風(fēng)暴潮耦合模型，得到如圖5所示的赤灣站與鹽田站水位驗(yàn)證曲線，可以看出，考慮波浪后的水位過(guò)程曲線更符合實(shí)際，最大誤差絕對(duì)值為13cm，平均相對(duì)誤差為11.3%，標(biāo)準(zhǔn)差為2.94，符合模擬精度要求。

圖5(a) 赤灣站不同模型下水位曲線

圖5(b) 鹽田站不同模型下水位曲線

3 深圳市沿岸漫堤風(fēng)險(xiǎn)評(píng)估

3.1 可能最大風(fēng)暴潮模型計(jì)算

為了尋求可能影響深圳市的最大風(fēng)暴潮對(duì)深圳市沿岸的影響，采用概率論法進(jìn)行計(jì)算可能最大臺(tái)風(fēng)中的近中心最大風(fēng)速Vmax，臺(tái)風(fēng)中心氣壓P0。

利用中國(guó)氣象局(CMA)1949-2015年熱帶氣旋最佳路徑集(Best-track)數(shù)據(jù)，選取出研究區(qū)域400km直徑范圍內(nèi)對(duì)本研究區(qū)域有重要影響的臺(tái)風(fēng)，以所選取臺(tái)風(fēng)的最小臺(tái)風(fēng)中心氣壓值作樣本。采用極值I型分布計(jì)算1000年一遇的臺(tái)風(fēng)中心氣壓值作為可能最大臺(tái)風(fēng)的中心氣壓，計(jì)算結(jié)果如圖6(a)所示，即深圳區(qū)域1000a遇的臺(tái)風(fēng)中心氣壓P0=880hPa。

圖6(a) 影響深圳的臺(tái)風(fēng)中心最低氣壓極值I型分布

圖6(b) 影響深圳的臺(tái)風(fēng)中心最大風(fēng)速極值I型分布

利用中國(guó)氣象局(CMA)1949-2015年熱帶氣旋最佳路徑集(Best-track)數(shù)據(jù)，取深圳市海灣區(qū)域400km直徑范圍內(nèi)歷年路經(jīng)本海區(qū)的臺(tái)風(fēng)最大風(fēng)速值作樣本，對(duì)于沒(méi)有臺(tái)風(fēng)經(jīng)過(guò)研究區(qū)域的年份，該年Vmax取為臺(tái)風(fēng)Vmax系列中的平均值。采用極值I型分布計(jì)算1000a一遇的最大風(fēng)速值作為可能最大臺(tái)風(fēng)的最大風(fēng)速值，計(jì)算結(jié)果如圖6(b)所示，深圳區(qū)域1000a一遇的臺(tái)風(fēng)中心最大風(fēng)速Vmax=85m/s。另外，臺(tái)風(fēng)路徑選取為東南向西北由大鵬灣登陸的路徑。

輸出堤前代表點(diǎn)(具體見(jiàn)圖7)的最高潮位與堤頂高程進(jìn)行比對(duì)，了解沿岸淹沒(méi)的情況，由于模式計(jì)算結(jié)果是相對(duì)平均海平面的，為便于比對(duì)，需要根據(jù)當(dāng)?shù)仄骄Ｆ矫媾c85高程的關(guān)系，將計(jì)算結(jié)果轉(zhuǎn)成以85黃海基面為起算面，其中西段提防采用赤灣站高程關(guān)系，即計(jì)算結(jié)果加上0.54m；東段提防采用鹽田站高程關(guān)系進(jìn)行轉(zhuǎn)換，即計(jì)算結(jié)果加上0.66m。表1列出堤前統(tǒng)計(jì)的最高潮位及堤頂高程。從代表點(diǎn)的最高潮位和海堤高程比較來(lái)看，在可能最大熱帶氣旋影響下，深圳沿岸均可能發(fā)生漫堤。

圖7 堤前代表點(diǎn)位置示意圖

表1 可能最大風(fēng)暴潮漫堤分析 (單位：m，基面：85高程)

續(xù)表1 可能最大風(fēng)暴潮漫堤分析 (單位：m，基面：85高程)

3.2 海堤防御標(biāo)準(zhǔn)

對(duì)于大部分沿海地區(qū)來(lái)說(shuō)，1000a一遇的標(biāo)準(zhǔn)是過(guò)高的。但是深圳作為特大型城市，其風(fēng)暴潮安全具有重要意義，因此考慮1000a一遇情形下的風(fēng)暴潮安全是有必要的，在此情形下得出的淹沒(méi)情況有利于應(yīng)急管理部門采取相應(yīng)的防范措施。

4 結(jié) 論

1)利用WRF中長(zhǎng)尺度氣象預(yù)報(bào)模式作為波浪-風(fēng)暴潮耦合模型的風(fēng)場(chǎng)邊界條件輸入，提高了風(fēng)暴潮模型的精確度

2)在珠江口波浪-風(fēng)暴潮模型驗(yàn)證實(shí)測(cè)資料準(zhǔn)確的基礎(chǔ)上，利用極值I型分布確定1000a一遇臺(tái)風(fēng)參數(shù)對(duì)影響深圳市海岸堤防未來(lái)發(fā)生的最大可能風(fēng)暴潮進(jìn)行計(jì)算分析。文章研究成果可為深圳市風(fēng)暴潮安全防范提供重要技術(shù)支撐，同時(shí)也為海洋水利防災(zāi)應(yīng)急部門提供決策參考。