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

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

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

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

1 模型的建立

1.1 臺(tái)風(fēng)場(chǎng)的構(gòu)建

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

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

(1)

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

(2)

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

(3)

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

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

(4)

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

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

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

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

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

2 模型驗(yàn)證

2.1 臺(tái)風(fēng)場(chǎng)驗(yàn)證

采用澳門(mén)氣象站(站位分布見(jiàn)圖2)“山竹”臺(tái)風(fēng)實(shí)測(cè)風(fēng)速風(fēng)向資料對(duì)臺(tái)風(fēng)模型計(jì)算的風(fēng)速風(fēng)向結(jié)果進(jìn)行驗(yàn)證，結(jié)果如圖3所示。

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

圖3 澳門(mén)氣象站“山竹”臺(tái)風(fēng)過(guò)程風(fēng)速風(fēng)向驗(yàn)證結(jié)果示意

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

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

2.2 天文潮驗(yàn)證

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

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

圖5 南海模型天文潮驗(yàn)證示意

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

2.3 風(fēng)暴潮水位過(guò)程驗(yàn)證

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

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

圖6 “山竹”臺(tái)風(fēng)過(guò)程水位驗(yàn)證示意

表1 誤差統(tǒng)計(jì)

3 結(jié)果分析

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

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

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

臺(tái)風(fēng)登陸后(見(jiàn)圖7b)，珠江口海域普遍為東南偏南風(fēng)，珠江河口東南向的開(kāi)口方向有利于水體向岸堆積并沿河道向上游輸運(yùn)，此時(shí)珠江口海域普遍達(dá)到增水的最大值。伶仃洋河口灣頂增水值較大約3.1 m，是由于伶仃洋河口灣喇叭狀的形態(tài)有利于水體向?yàn)稠斁奂Ｆ溆嗫陂T(mén)區(qū)域最大增水值普遍為2.6～3.1 m。

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

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

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

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

圖8 各代表站不同水位過(guò)程

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

4 結(jié)語(yǔ)

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

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