粵港澳大灣區(qū)風(fēng)暴潮數(shù)值模型的建立與應(yīng)用

羅志發(fā)，黃本勝，譚 超，黃廣靈

(1.廣東省水利水電科學(xué)研究院，廣東 廣州 510635；2.廣東省水動力學(xué)應(yīng)用研究重點實驗室，廣東 廣州 510635；3. 廣東省流域水環(huán)境治理與水生態(tài)修復(fù)重點實驗室，廣東 廣州 510635；4. 中山大學(xué)，廣東 廣州 510275)

粵港澳大灣區(qū)是我國開放程度最高、經(jīng)濟活力最強的區(qū)域之一，以不足全國1%的土地面積和不足全國5%的人口，貢獻了全國經(jīng)濟總量的17%，在國家發(fā)展大局中具有重要戰(zhàn)略地位?；浉郯拇鬄硡^(qū)位于中國大陸南端，瀕臨南海，地處珠江流域下游，河網(wǎng)水系發(fā)達。由于特殊的地理位置和氣候，導(dǎo)致大灣區(qū)臺風(fēng)暴潮災(zāi)害易發(fā)頻發(fā)，由此造成的人員傷亡、經(jīng)濟損失相當(dāng)巨大，已經(jīng)成為影響人民生活質(zhì)量、制約國民經(jīng)濟高質(zhì)量發(fā)展的重要因素。開展粵港澳大灣區(qū)風(fēng)暴潮研究既是重要科學(xué)問題也是風(fēng)暴潮預(yù)警預(yù)報關(guān)鍵技術(shù)研發(fā)的基礎(chǔ)，對提升大灣區(qū)水安全保障能力具有重要的意義。

粵港澳大灣區(qū)所在的珠江河口是多種動力因子協(xié)同作用的復(fù)雜系統(tǒng)，呈“三江匯流，八口入?！钡男蝿?。河網(wǎng)區(qū)水網(wǎng)密布，橫向支汊發(fā)育，其間多種動力因子相互耦合，動力復(fù)雜多變。珠江河口風(fēng)暴潮數(shù)值模擬已開展較多的研究，如二維風(fēng)暴潮模型的建立[1]，路徑、風(fēng)速對風(fēng)暴潮增水的影響[2]，地形對局部增水的影響等[3]。針對珠江河口復(fù)雜的動力系統(tǒng)，本文基于無結(jié)構(gòu)三角網(wǎng)格的數(shù)值模式，優(yōu)化擬合復(fù)雜岸形，通過二重嵌套的方式合理提供外海余水位邊界，構(gòu)造臺風(fēng)風(fēng)場和氣壓場作為模型的海表面邊界條件，構(gòu)建珠江河口三維風(fēng)暴潮數(shù)值模型。對多組風(fēng)暴潮增水過程進行數(shù)值模擬，驗證結(jié)果較好。選取1822號“山竹”臺風(fēng)的風(fēng)暴潮增水過程進行數(shù)值模擬，研究天文潮與風(fēng)暴潮非線性作用對風(fēng)暴潮增水的影響。

1 模型的建立

1.1 臺風(fēng)場的構(gòu)建

臺風(fēng)是風(fēng)暴潮模擬的重要因子，臺風(fēng)風(fēng)場的質(zhì)量決定了風(fēng)暴潮水位模擬的準(zhǔn)確性。利用臺風(fēng)風(fēng)場模型可以較好地模擬臺風(fēng)中心區(qū)域的氣旋風(fēng)，但是離臺風(fēng)中心較遠的區(qū)域，由于氣旋風(fēng)場的減小，背景風(fēng)場的影響逐漸增大，需要利用再分析風(fēng)場予以彌補。再分析風(fēng)場能很好的反映大范圍的風(fēng)場結(jié)構(gòu)特征，但由于其空間分辨率較低，難以刻畫臺風(fēng)中心的風(fēng)場結(jié)構(gòu)，其風(fēng)速值顯著小于實際風(fēng)速。針對臺風(fēng)風(fēng)場模型和再分析風(fēng)場資料存在的不足，本文結(jié)合兩者優(yōu)勢進行臺風(fēng)風(fēng)場的構(gòu)造，即克服了經(jīng)驗臺風(fēng)公式外圍風(fēng)場計算結(jié)果較小的問題，也彌補了再分析風(fēng)場臺風(fēng)中心風(fēng)強不足的缺陷。

臺風(fēng)風(fēng)場模型根據(jù)梯度風(fēng)原理，由臺風(fēng)氣壓場計算出風(fēng)場。本文采用國內(nèi)外應(yīng)用較廣泛的Fujita[4]、Ueno[5]公式構(gòu)建氣壓場和臺風(fēng)風(fēng)場，其式如下：

(1)

·[(x-x0)sinθ+(y-y0)cosθ]

(2)

·[(x-x0)cosθ-(y-y0)sinθ]

(3)

式中P(r)為距離臺風(fēng)中心r處的氣壓值，hPa；P∞為臺風(fēng)外圍無窮遠處的大氣壓，取1 010 hPa；P0為臺風(fēng)中心氣壓，hPa；R為臺風(fēng)最大風(fēng)速半徑，km，根據(jù)經(jīng)驗公式計算：R=Rk-0.4×(P0-900)+0.01×(P0-900)2，Rk為經(jīng)驗常數(shù)，介于30～60，本文取40；r為計算點離臺風(fēng)中心的距離，km；Vdx、Vdy分別為臺風(fēng)移動速度在x、y的分量，m/s；f為科氏力參數(shù)；ρα為空氣密度，取值1.292 9 g/m3；ΔP=P∞-P0為臺風(fēng)中心氣壓示度，hPa；x0、y0為臺風(fēng)中心坐標(biāo)；θ為臺風(fēng)流入角，取20°；C1、C2為訂正系數(shù)，本文取0.8。文中構(gòu)建臺風(fēng)風(fēng)場模型需要的臺風(fēng)相關(guān)數(shù)據(jù)來自于中央氣象臺臺風(fēng)網(wǎng)發(fā)布的信息。

外圍背景風(fēng)場采用歐洲中期數(shù)值預(yù)報中心的ERA50再分析風(fēng)場，分辨率為0.125°×0.125°。臺風(fēng)模型風(fēng)場與ERA50再分析風(fēng)場合成方法如下：

(4)

1.2 風(fēng)暴潮模型的建立

風(fēng)暴潮增水是大尺度的動力過程，模擬好風(fēng)暴潮需要建立大范圍的模型，同時，珠江口區(qū)域精細化模擬需要構(gòu)建高分辨率的計算網(wǎng)格，為了解決模擬范圍、網(wǎng)格分辨率、計算效率的問題，本文采用兩重嵌套的方法進行模擬計算(見圖1所示)。選用海洋環(huán)流模式SELFE建立風(fēng)暴潮模型，該模式基于非結(jié)構(gòu)化網(wǎng)格，可精細化擬合復(fù)雜岸線和地形，采用半隱式的歐拉—拉格朗日有限元算法求解N-S方程組。模式最大的特點是減小CFL條件的限制，在保證計算結(jié)果準(zhǔn)確的前提下，可以適當(dāng)放大時間步長，以達到計算精度和計算效率的雙贏。

圖1 模型嵌套計算網(wǎng)格示意

大范圍的南海模型計算范圍為98°E～126°E，0°N～30°N，涵蓋整個南海及西北太平洋海域，南邊界至卡里馬塔海峽、北邊界至浙江省沿岸海域、東邊界至48 h警戒線。網(wǎng)格分辨率從近岸的1 km逐漸過渡到外海20 km，水深數(shù)據(jù)采用ETOP01全球1′×1′分辨率的地形資料。珠江河口計算范圍為111°E～ 116°E， 21°N～ 23.7°N，模型上邊界為西江高要，北江石角，東江博羅，流溪河及潭江上游。外海下邊界取約100 m等深線處。模型采用無結(jié)構(gòu)網(wǎng)格，擬合復(fù)雜河岸邊界，對局部進行加密，提高網(wǎng)格分辨率。模型網(wǎng)格共有101 752個節(jié)點，173 045個網(wǎng)格單元，網(wǎng)格大小為從河網(wǎng)區(qū)10 m逐漸過渡到外海的20 km。模型在垂向上采用Sigma坐標(biāo)，均勻分為10層。模型中珠江三角洲網(wǎng)河區(qū)采用2005—2008年的大范圍實測地形，河口區(qū)及近岸海區(qū)采用2000—2008年海圖地形，外海海域采用ETOP01全球海洋地形。

南海模型由風(fēng)場、氣壓場及8個主要分潮(M2、S2、N2、K2、K1、O1、P1、Q1)進行驅(qū)動，計算得到余水位，將余水位以及風(fēng)場、氣壓場、8個分潮作為珠江口模型的驅(qū)動條件，對珠江河口風(fēng)暴潮增水進行數(shù)值模擬。

2 模型驗證

2.1 臺風(fēng)場驗證

采用澳門氣象站(站位分布見圖2)“山竹”臺風(fēng)實測風(fēng)速風(fēng)向資料對臺風(fēng)模型計算的風(fēng)速風(fēng)向結(jié)果進行驗證，結(jié)果如圖3所示。

圖2 水文氣象站位分布示意

圖3 澳門氣象站“山竹”臺風(fēng)過程風(fēng)速風(fēng)向驗證結(jié)果示意

計算的臺風(fēng)風(fēng)速、風(fēng)向其變化過程與實測數(shù)據(jù)相符合，最大風(fēng)速絕對誤差僅1 m/s。說明計算的風(fēng)場較好地反映了臺風(fēng)過程及最大風(fēng)速。圖4為采用本文臺風(fēng)場構(gòu)建方法計算的臺風(fēng)場圖，由圖4可知，在經(jīng)驗?zāi)Ｐ偷幕A(chǔ)上融合了ERA50的再分析數(shù)據(jù)，對臺風(fēng)場的計算避免了經(jīng)驗公式關(guān)于臺風(fēng)外圍風(fēng)速計算明顯偏小的問題，構(gòu)建的風(fēng)場結(jié)構(gòu)較為合理，可用于風(fēng)暴潮數(shù)值模擬。

圖4 2018年9月16日6：00“山竹”臺風(fēng)風(fēng)場示意

2.2 天文潮驗證

模擬風(fēng)暴潮和天文潮耦合的水位過程時，天文潮的模擬計算是基礎(chǔ)。在不考慮風(fēng)應(yīng)力的情況下，對南海天文潮進行數(shù)值模擬，以檢驗?zāi)Ｐ蛯μ煳某庇嬎愕目煽啃?。計算時間為2018年7月1—31日，選用沿海大萬山、珠海港、上川島、馬鞭洲等潮位站(站位分布見圖2)實測資料對模型計算結(jié)果進行驗證，比較結(jié)果見圖5所示。

由圖5可知，各潮位站計算的天文潮水位過程線，無論是高低潮位還是相位與實測數(shù)據(jù)基本一致，潮位計算絕對誤差介于8.8～18.4 cm，驗證結(jié)果良好。表明模型開邊界條件基本正確，底摩擦參數(shù)設(shè)置合理，可在此基礎(chǔ)上進一步開發(fā)風(fēng)暴潮模式，為珠江口風(fēng)暴潮模型提供合理的余水位開邊界條件。

圖5 南海模型天文潮驗證示意

珠江口模型已經(jīng)過2001年、2005年、2008年、2009年等多組水文觀測資料的驗證，并成功應(yīng)用于珠江河口三維水動力及鹽度輸移的數(shù)值模擬，模式簡介及天文潮、水動力驗證結(jié)果詳見文獻[6]。

2.3 風(fēng)暴潮水位過程驗證

以“山竹”臺風(fēng)為例，對珠江河口風(fēng)暴潮模型進行驗證。模型上游高要、石角、博羅、人和給定實測徑流數(shù)據(jù)，石咀因缺乏實測流量數(shù)據(jù)而給定當(dāng)月份的多年平均值。外海開邊界水位通過天文潮水位+余水位的形式給定，天文潮水位由調(diào)和常數(shù)計算得到，余水位由南海模型提供。海表面臺風(fēng)場及氣壓場通過臺風(fēng)模型及氣壓模型給出。模型時間步長設(shè)置為200 s，模擬時間為2018年9月10—20日，前5 d用于天文潮計算的穩(wěn)定，后5 d對風(fēng)暴潮過程進行模擬。

選用珠江口門及河網(wǎng)區(qū)水文站位(泗盛圍、南沙、萬頃沙、橫門、燈籠山、黃金、西炮臺，站位分布見圖2)實測水位資料對風(fēng)暴潮計算結(jié)果進行驗證。由于篇幅有限，僅給出部分驗證結(jié)果圖(見圖6)。由圖6可知，計算結(jié)果水位變化過程與實測數(shù)據(jù)水位變化趨勢一致，計算結(jié)果與實測結(jié)果吻合較好。統(tǒng)計了模型計算最高潮位的絕對誤差以及相位誤差(見表1)，由表1可知，“山竹”臺風(fēng)計算結(jié)果的最高潮位絕對誤差介于0.15～0.49 m，相對誤差控制在25%以內(nèi)，最高潮位的相位誤差均在1 h以內(nèi)。綜上分析，模型計算結(jié)果較為合理，所建立的模型能夠較好地反映珠江河口徑流、潮流、風(fēng)等多種動力因子相互作用的結(jié)果，能夠較為準(zhǔn)確地模擬珠江河口的風(fēng)暴潮過程。

圖6 “山竹”臺風(fēng)過程水位驗證示意

表1 誤差統(tǒng)計

3 結(jié)果分析

本節(jié)利用上文已建立的風(fēng)暴潮模型，模擬“山竹”臺風(fēng)風(fēng)暴潮增水過程，討論天文潮與風(fēng)暴潮非線性相互作用的特征。

3.1 風(fēng)暴潮增減水特征分析

2018年9月7日20:00，臺風(fēng)“山竹”在西北太平洋洋面上生成；9月15日，臺風(fēng)“山竹”從菲律賓北部登陸；16日17:00在廣東臺山海宴鎮(zhèn)登陸，登陸時中心附近最大風(fēng)力14級，中心最低氣壓為955 hPa。9月16日13:00，臺風(fēng)中心位于珠江口南側(cè)(見圖7a)，伶仃洋海域為東北偏東風(fēng)，伶仃洋東岸在離岸風(fēng)的驅(qū)動下，產(chǎn)生減水，深圳灣區(qū)域最大減水值約-1.5 m。伶仃洋西岸由于水體的橫向堆積，產(chǎn)生0.5～1.0 m的增水。磨刀門、黃茅海區(qū)域為偏北風(fēng)，有利于水體離岸輸運，形成-0.5～-1.0 m的減水。

臺風(fēng)登陸后(見圖7b)，珠江口海域普遍為東南偏南風(fēng)，珠江河口東南向的開口方向有利于水體向岸堆積并沿河道向上游輸運，此時珠江口海域普遍達到增水的最大值。伶仃洋河口灣頂增水值較大約3.1 m，是由于伶仃洋河口灣喇叭狀的形態(tài)有利于水體向灣頂聚集。其余口門區(qū)域最大增水值普遍為2.6～3.1 m。

注：a、b分別為2018年9月16日13時以及2018年9月16日17時珠江口風(fēng)暴潮增減水分布，紅色線條為臺風(fēng)路徑，紅色圓點為臺風(fēng)中心。

3.2 天文潮風(fēng)暴潮非線性相互作用分析

對于風(fēng)暴潮增水而言，除了氣象因子之外，天文潮與風(fēng)的非線性作用對風(fēng)暴潮增水也是重要的影響因子。設(shè)置3個數(shù)值試驗：考慮天文潮風(fēng)的耦合作用、僅考慮潮的作用、僅考慮風(fēng)的作用。利用天文潮與風(fēng)耦合計算的總水位減去天文潮位可得到風(fēng)暴潮增水，風(fēng)暴潮增水包含了純風(fēng)生增水(僅考慮風(fēng)的作用)和非線性增水。利用風(fēng)暴潮增水減去純風(fēng)生增水后可得到非線性增水。在珠江口區(qū)域設(shè)置4個計算點(見圖7a)，分別位于虎門(P1)、伶仃洋西岸(P2)、磨刀門(P3)、黃茅海河口灣頂(P4)，分析各計算點不同水位的變化特征。

由風(fēng)暴潮水位的時間序列(見圖8)可以看出，由于風(fēng)暴潮最大增水發(fā)生在小潮的低水位階段，因此總水位最大值主要是由風(fēng)暴潮最大增水決定。除了P2點外，其余各計算點均呈現(xiàn)先減水后增水的特征，其中最大增水值位于P1點約3.1 m，與伶仃洋河口灣喇叭狀的形態(tài)有利于水體聚集有關(guān)；最小減水值位于P3點約-2 m，與黃茅海河口灣的開口方向和臺風(fēng)登陸前的風(fēng)向有關(guān)。最大減水發(fā)生在臺風(fēng)登陸前1～2 h，最大增水值發(fā)生在臺風(fēng)登陸后1～4 h。由于P2點位于伶仃洋西岸，臺風(fēng)登陸前后的風(fēng)向均有利于水體的向岸堆積，因而在整個臺風(fēng)過程中均處于增水的狀態(tài)，也是最早出現(xiàn)最大增水的位置。

圖8 各代表站不同水位過程

各計算點非線性增水約1.19～1.39 m 對于風(fēng)暴潮增水的貢獻介于37.5%～50.0 %之間。由天文潮水位曲線與非線性增水曲線可以看出，非線性增水極小值一般出現(xiàn)在高潮位階段不利于增水，極大值一般出現(xiàn)在低潮位階段有利于增水。計算結(jié)果與前人研究結(jié)果基本一致[7]，從動力學(xué)上的解釋認為：耦合作用中風(fēng)與水位效應(yīng)的非線性項作用影響的強弱與潮汐、純風(fēng)暴潮和水深都有關(guān)系，其中水深的影響最大，非線性項的作用影響與水深成反比；潮汐高潮削弱了風(fēng)應(yīng)力作用，從而減弱了風(fēng)暴潮，低潮時相反[8]。珠江河口天文潮與風(fēng)耦合作用較為復(fù)雜，與理論推導(dǎo)的結(jié)果有所差異，非線性增水的極值與高低潮位存在一定的相位差。

4 結(jié)語

1) 本文基于SELFE模式和圓形臺風(fēng)風(fēng)場，考慮了徑流、天文潮、臺風(fēng)的耦合作用，建立了粵港澳大灣區(qū)風(fēng)暴潮數(shù)值模型。采用實測水文氣象資料對“山竹”進行模擬驗證，計算最大水位相對誤差控制在25%以內(nèi)，模型計算結(jié)果較好地反映了風(fēng)暴潮增水的過程，可用于粵港澳大灣區(qū)風(fēng)暴潮數(shù)值模擬研究。

2) 基于本文建立的風(fēng)暴潮數(shù)值模型，以1822號臺風(fēng)“山竹”為例，討論了非線性增水特征，結(jié)果表明非線性增水對風(fēng)暴潮增水起到重要的作用，非線性增水極小值出現(xiàn)在高潮位階段不利于風(fēng)暴潮增水，極大值出現(xiàn)在低潮位階段有利于風(fēng)暴潮增水。