珠江口伶仃洋風暴潮傳播特性及影響因素數值模擬研究
——以臺風“山竹”為例

陳鵬，張卓,，宋志堯，章衛勝，葉榮輝，李玉婷

（1.南京師范大學 虛擬地理環境教育部重點實驗室，江蘇 南京 210023；2.江蘇省地理信息資源開發與利用協同創新中心，江蘇 南京 210023；3.水利部水旱災害防御重點實驗室，江蘇 南京 210029；4.珠江水利委員會珠江水利科學研究院，廣東 廣州 510611；5.水利部珠江河口動力學及伴生過程調控重點實驗室，廣東 廣州 510611）

風暴潮是一種重要的天氣現象，會導致沿海地區洪水泛濫、屋舍沖毀、農田淹沒等災害，威脅人民的生命安全且造成巨大的財產損失[1]。珠江口面臨南海，是臺風的多發區，珠江三角洲又是中國重要的經濟區域，風暴潮災害受到普遍的關注和重視。據統計，僅2018 年臺風“山竹”就造成了廣東、廣西、海南、湖南、貴州5 省（自治區）近300 萬人受災，直接經濟損失超50 億元[2]。因此，風暴潮傳播規律的研究對于災害區域的防災減災及應急響應措施的制定是非常必要的[3]。

通過理論推導，Proudman 認為風暴潮增水最大值應出現在天文潮潮位較低時刻[4]；Prandle 等在研究泰晤士河風暴潮增水時發現，風暴潮增水在天文潮漲潮階段較落潮時偏大[5]。國內，很多學者對我國沿海地區的風暴潮開展了大量研究。在觀測研究方面，王永信等利用40 年的風暴潮增水資料，分析了風暴潮沿珠江河道上溯運動問題[6]；史鍵輝等根據赤灣等站點的實測資料，研究了風暴潮和天文潮之間的相互作用[7]；陳杏文等利用浮標觀測數據分析了2021 年臺風期間海洋動力學和熱力學的響應特征[8]；羅志發等利用1970－2018年共49 年的珠江口實測數據，結合數值模擬結果分析得到各測點歷史最高潮位以0.02～0.03 m/a 的平均速率增加[9]；劉媛媛等將風速、風向、氣壓以及前時序的潮位數據與LSTM 模型結合，建立了風暴潮預報模型[10]。在數值模擬方面，于斌等對風暴潮沿珠江河道上溯進行模擬，發現風暴潮位和風暴強度和路徑有很大關系[11]；李杰等發現當臺風中心距離驗潮站最短時該站通常會出現最大增水[12]；Du 等對珠江口不同臺風路徑及臺風中心移動速度對風暴潮增水和海岸淹沒的影響進行比較，發現當臺風移動速度越慢時引發的海岸淹沒越嚴重[13]；嚴楓等對雙臺風作用下的增水、流場的相互影響及影響區域進行了探討[14]。近十年來，珠江口風暴潮數值也在逐步完善中，這方面有傅賜福等對風場適應性研究[15]；葉榮輝等的珠江三角洲風暴潮模型[16]；李心雨等研究了WRF 與臺風經驗模型在臺風“山竹”匯總的應用對比，認為臺風經驗模型在臺風模擬中仍是可接受且便捷的方法[17]。在未來氣候變化影響研究方面，Yin 等研究了未來海平面上升和臺風強度增強對珠江口風暴潮的影響[18]。

總體來說，目前從臺風登陸位置和移動方向對珠江口風暴增水時空分布影響的研究不少，但多是實測點的風暴潮增水分析，而采樣點數和臺風案例有限，無法對影響因素實行單因子定量分析。數值模擬多針對整個珠江口，而針對西進型臺風在廣東近海改變登陸地點、登陸角度以及登陸時間條件下，對伶仃洋內風暴潮增水的時空分布特征變化進行影響因子量化分析的相關研究較少，更缺乏伶仃洋內天文潮與風暴潮的相互作用對風暴潮增水影響的典型案例。因此，研究選取近年來對廣東地區影響較大的一次西進型臺風——“山竹”作為典型案例進行數值模擬，分析西進型臺風在伶仃洋內的風暴潮分布及變化特征，探討西進型臺風改變條件（如登陸時間、中心氣壓、登陸位置和移動方向）對伶仃洋風暴潮增水的響應規律，為有關部門的決策制定提供依據。

1 模型介紹

1.1 研究區域和計算范圍

本文中，為了獲取熱帶氣旋引起風暴潮的整個過程，模型計算采用大、小區雙層嵌套的方法。大區范圍包括南海和部分西北太平洋區域，模型邊界采用主要4 個半日潮（M2，S2，K2，N2）和4 個全日潮（K1，O1，Q1，P1）來估算外海邊界上的天文潮。小區范圍（圖1）包括珠江三角洲河網、八大入海口門以及伶仃洋和黃茅海，簡稱珠江口模型，模型采用無結構三角網格，網格分辨率在外海為5 km×5 km，河口海灣內加密為100 m×100 m。該模型在北江、西江和東江的上游設置流量邊界，外海開邊界由大范圍計算的水位結果提供。研究的主要區域位于伶仃洋內，該區域口門外受到香港離島及珠海東澳島等眾多島嶼的掩護，風浪對于潮位的影響相對較小[19]，因此本文不考慮風浪對風暴潮的影響。

圖1 珠江口模型網格和地形

1.2 臺風模型

選取登陸華南地區的最強臺風—“山竹”作為典型案例進行模擬（圖2），臺風“山竹”在西北太平洋洋面生成，并于2018年9月16日17時在廣東臺山登陸，登陸時仍有14 級風力。“山竹”引發的大風、風暴潮以及暴雨洪水等對廣東、香港和澳門等地造成了相當大的影響，據媒體報道，至少有70人死亡，超過252萬人被轉移安置。

圖2 臺風“山竹”行進路徑及在登陸時刻圖

本文采用Holland[20]模型作為風場驅動，與其他臺風氣壓參數模型相比，Holland 模型引入氣壓剖面特征參數B 來描述相同的強度和半徑上的風速分布。其值可由Holland 模型和梯度風公式結合日本氣象廳公布的50 節（25.722 m/s）風速圈半徑估算得到。臺風“山竹”在移動和發展過程中，估算B的變化范圍在0.9～1.0之間，本文取B=1.0。

1.3 風暴潮模型

FVCOM 是一個非結構網格、有限體積的海岸和海洋模型，最初由陳長勝團隊開發[21]，并由UMASDSD 和WHOI 研 究 團 隊 改 進[22-23]。基 于FVCOM 主程序構建的風暴潮模型，采用干濕處理和通用湍流模塊，其中驅動力有外海潮邊界，上游徑流以及海表面風應力。本文構建中國海-珠江口風暴增水計算模型，該模型既能模擬天文潮和風暴潮的相互作用，又可以模擬波流相互作用，能夠得到比前人研究更合理的結果[24]。珠江口模型的外海邊界可由中國海模型的結果提供，上游徑流洪水為西江的高要、北江的石角和東江的博羅等3個水文站當日的實測徑流量。

風應力是海氣能量交互的關鍵，在數值模擬中風應力一般按海表面10 m風速的2次方來計算：

式中：τ為海表面切應力；ρa為空氣密度；Cd為拖曳力系數；U10為海面10 m風速。其中風拖曳力系數采用Wu[25]基于實測資料得到的經驗公式計算得到，公式如下：

2 模型驗證與結果分析

2.1 風場驗證

對風場的驗證采用收集到的臺風登陸前后的3個實測站位（W1、W2和W3）的風速風向資料，站位分布見圖3。其中W1 站點位于澳門機場東南側，W2 號站點位于磨刀門口門處，W3 號站點位于黃茅海灣口處。由驗證結果可知，“山竹”臺風期間風速風向的模擬結果較好（圖4），說明本文選取的臺風風場計算方案合理。

圖3 珠江口模型風速風向及水位測站分布圖

圖4 臺風登陸前后風速風向計算值（線）與實測值（點）驗證結果

2.2 水位驗證

風暴潮模擬計算從2018 年9 月14 日0 時至2018 年9 月18 日0 時，歷時4 天，涵蓋了臺風“山竹”從臨近、登陸和消亡的全過程。采用6 個潮位測站（圖3）的實測資料對珠江口風暴潮模型進行驗證，通過對糙率和臺風參數的率定，除西側萬頃沙站的風暴增水退潮階段相對誤差較大外，其他站位模擬結果與實測值吻合較好（圖5）。萬頃沙站退潮階段誤差較大的原因可能是沒有考慮珠江三角區內臺風暴雨引發的匯流洪水，也可能是該站位上游的河道地形有誤所致。

圖5 臺風期間水位過程計算值（線）與實測值（點）驗證結果

2.3 臺風期間的風暴增水分析

伶仃洋南部灣口寬約35 km，北部灣頂寬僅有5 km，呈口外寬里灣窄的喇叭狀。在喇叭狀海灣中，風暴潮從口外深海向淺水灣頂傳播時，受到河口斷面收縮和底摩擦的影響，潮波能力向灣頂聚集，風暴增水被放大，從而在灣頂形成極端高水[26]。珠江口模型計算的風暴潮總水位減去天文潮可以得到風暴增水幾個關鍵時刻的空間分布（圖6）。臺風登陸前8小時，東北風掃過海灣，吹向西南，導致伶仃洋灣頂和虎門的風暴增水為負。隨著臺風中心向珠江口靠近，風向逐漸沿順時針旋轉，16 日13－16 時，在風力驅動下水體向伶仃洋西海岸堆積，伶仃洋出現西高東低的增水分布。登陸時風力達到最大，西北風向幾乎與海灣走向平行，海岸及河口附近的多個觀測站的增水均達到了最大值。登陸后臺風強度明顯減弱，風暴增水也隨之減弱。

圖6 臺風“山竹”登陸期間的增水過程

2.4 風暴潮水位過程線分析

為研究臺風登陸前后風暴增水過程的變化，在伶仃洋的灣口、中部、灣頂及虎門河口選取12個特征點（圖7（a）），記錄特征點逐時的風暴潮位變化。從特征點的風暴潮水位過程線圖（7（b））可以看出，12 個特征點的風暴潮過程曲線相似，可分為初振段、主振段和余振段。特征點所在位置的不同，各曲線顯示出不同的特征。

圖7 伶仃洋和虎門口12個特征點在各階段增水過程（臺風在第65 h登陸）

初振段之前，一系列較弱增水波傳入伶仃洋，持續時間約30 h（第17～48 h），振幅約為0.2 m（圖7（c）、（f））。這種先兆增水是由大洋或外海遠程風引發的邊緣波形成的，在靠近海灣時會略有增強。在第48～60 h 的增水上升段，灣口附近的1～5點（圖7（c））和灣中至灣頂6～12點（圖7（f））的水位線特征不同，主要區別在于后者上升段疊加了更強的短期振蕩，表明在灣內天文潮對風暴增水的影響更大。在增水的峰值階段（圖7（d）、（g）），由于伶仃洋為喇叭狀，風暴引起的風暴增水可以無阻擋地向灣內堆積，從灣口至灣頂的峰值增水逐漸加大。同時，灣內不同點之間的增水峰值出現的時間也有明顯差異。因伶仃洋是開敞式灣口，增水水位沿橫向坡降較小，如1 和2、6和8、9和10點，而4點和5點由于受到島嶼遮擋，此兩點的峰值增水小于3 點。主增水階段過后，灣口附近（1～5）和灣內（6～12）之間的風暴增水下降段有所不同（圖7（e）、（h））。與上升段相似，特征點在灣內的短期波動比海灣口附近的更強烈。在峰值增水出現約10 h后，第75～77 h，灣內各點出現明顯的增水次高峰，而灣口附近各點則不明顯。

從對圖6、圖7 的分析中，可以得到珠江口風暴潮的時空分布受臺風風場和海岸幾何形態及水深地形的共同影響。由于海灣開口寬闊，灣內外水量交換通暢，故12 個特征點的水位過程線形態趨于一致，而伶仃洋灣內比灣口附近點的水位過程線的波動現象更強烈，這是地形淺化效應和天文潮與風暴潮相互作用的結果。

3 珠江口伶仃洋風暴潮影響因子分析

3.1 臺風登陸時間對風暴增水的影響

臺風“山竹”引發的風暴增水在9月16日17時（臺風登陸時刻）到達峰值，此時的天文潮處于落潮階段。臺風的不同登陸時間會對風暴增水造成的影響與天文潮和風暴增水的相互作用有關，也是風暴潮研究的焦點問題。Proudman 對天文潮與風暴潮相互作用對風暴增水的影響進行了理論研究，發現在高潮位附近的風暴增水總是小于低潮附近的風暴增水[4]，為驗證珠江口風暴潮是否遵循該規律，設計如下試驗：將臺風“山竹”的時間序列分別提前1、2、3 和4 h，以+1、+2、+3、+4來表示，和滯后1 和2 h，以-1、-2 來表示，利用修正后的風場進行數值試驗。取代表灣口的4點和代表灣內的11 點的總水位線過程線來說明天文潮和風暴潮相互作用的影響。如圖8 所示，總水位最大值與潮位成正比，這意味著當風暴潮與潮位同步時可能會出現極大的總水位。而去除掉天文潮后的風暴潮增水卻與潮位并非成正比（圖8（c）、（d）），而與之相反，在高潮位（+2）登陸的風暴增水小于在低潮（0、-1 和-2）登陸的風暴增水。通過對4 點和11 點的比較，表明灣內的風暴增水符合Proudman 的分析結論[4]，即在高潮位出現風暴增水最小，并且灣頂附近的增水受天文潮的影響更加顯著。經過分析認為，這種由于風暴潮和天文潮之間相位差造成的增水差異，主要源于非線性阻力和淺水效應。因此，在對風暴潮增水進行預報時，尤其在灣內淺海區域，必須考慮臺風登陸期間天文潮相位的影響。

圖8 4點和11點風暴潮總水位（上）和增水（下）的時間變化

3.2 臺風中心最低氣壓的影響

臺風強度與中心氣壓差成正比，為了定量分析不同的中心氣壓差對增水大小的影響，將“山竹”臺風的中心氣壓原始數據進行了±10 和±20 hPa 的偏移（圖9（a）），對四種設計情景的風暴增水過程進行模擬。從4 點（圖9（b））和11 點（圖9（c））設計試驗情景和原始臺風之間的總水位對比表明，增大壓差（-10，-20）既能增大正、負風暴潮位的峰值，又能延長高水位持續時間。反之，則峰值水位降低，且高水位持續時間縮短。對比4點和11點的水位過程線，中心氣壓差變化對灣內11點的影響比灣口的4點更為顯著，說明風暴潮增水過程對同一臺風氣壓差變化在空間上的分布是不均勻的。一般而言，水深較淺的沿岸和灣內的風暴增水對臺風強度（或氣壓差）的變化更為敏感，這與淺水效應和灣頂潮能放大效應有關。

圖9 4點和11點在不同臺風中心氣壓時間序列下總水位對比

3.3 登陸地點的影響

有學者探討了臺風風速強弱、尺度大小以及海岸線的幾何形態等因素對伶仃洋風暴潮的影響，但對于西行臺風登陸位置的影響探討較少。下面以臺風“山竹”為例，研究登陸位置的變化對伶仃洋風暴潮的影響。實驗設計是將“山竹”臺風的原路徑向左移動4 次、右移動1 次，形成5 條新路徑，每條臺風路徑相差0.5°（圖10）。

圖10 臺風設計路徑與原始路徑及登陸位置圖

利用中國海-珠江風暴潮模型，模擬5 條設計路徑的臺風登陸時伶仃洋風暴增水的空間分布，并與原始路徑進行比較。模擬結果表明，臺風路徑向西移動0.5°和1°時，風暴潮增水的強度和范圍均較大。繼續向西移動幅度達1.5°和2°時，風暴潮增水轉向逐漸減弱。從圖11（b）看出，路徑西移0.5°時伶仃洋內產生的增水最大，超過原路徑，因為此時臺風最大風速半徑正好位于伶仃洋口門中心區域，大量海水在風的驅動涌入灣內。相對地，原路徑向東移動0.5°時，伶仃洋內的風暴增水強度則急劇減弱，這是因為最大風速半徑已經位于伶仃洋灣口的東部，東移后灣口處的風力比原路徑更弱。圖12 顯示了臺風路徑變化對風暴潮增水過程的影響，可以看出登陸位置的不同對增水過程有很大的影響。具體而言，路徑西移使增水上升段提前到來，路徑東移則推遲增水上升段的時間，且增水下降段路徑東移和西移存在較大差別，分析原因是路徑西移使最大風速半徑圈提前經過伶仃洋灣口，因此提前引發增水，東移則相反。

圖11 伶仃洋灣最大風暴增水分布

圖12 4點和11點在不同臺風路徑下的風暴潮增水過程線

3.4 登陸角度的影響

根據1968－2017 年的統計資料，對珠三角區域影響較大的臺風主要是西進型臺風，即以相對平直的西北向移動路徑靠近海岸并登陸的臺風類型。這種臺風的特點是路徑長且行進速度快，誕生于北太平洋中西部，較長的生命周期使其在海面充分發展，等到臨近登陸時往往形成超強臺風，2018 年的超強臺風“山竹”就是其中一個典型案例。西進型臺風的路線往往比較平直，因此可以用一個登陸角度（登陸時臺風路徑跟海岸線切線方向的夾角）來代表臺風的移動方向。為了研究登陸角度變化對于伶仃洋內風暴潮增水的影響，在臺風“山竹”原路徑的基礎上，將整個路徑分別沿逆時針（減小登陸角度）和順時針（增大登陸角度）方向旋轉5°，得到兩個設計臺風路徑（圖13）并保持其他臺風參數不變。從圖14 的模擬結果可知，減小登陸角度使伶仃洋風暴增水增加，而增加登陸角度使伶仃洋風暴潮增水明顯減小。分析原因，認為廣東沿海海岸線基本沿東西向，當臺風登陸時與海岸的夾角減小時，臺風最大風圈跟伶仃洋之間的距離更小，從而會增大伶仃洋的最大風暴潮影響的面積。從圖15 可以看到，當登陸角度變小，使風暴潮增水前期（60 h以前）增水變小，而在峰值及峰值之后增水增加。增大登陸角度，風暴潮前期增水變大，但在峰值之后增水明顯變小，這都是由于登陸角度改變伶仃洋風場過程造成的。

圖13 原始路徑（方案1）與設計路徑（方案2、3）分布圖

圖14 原始路徑（a）與設計路徑（b）、（c）風暴潮最大增水空間分布

圖15 4點和11點不同路徑下風暴增水過程線

4 結論

本文采用FVCOM 數值模型，以近十年來登陸廣東省的超強臺風“山竹”為例，對伶仃洋及其附近海域的風暴潮過程進行模擬和分析。模型計算結果與珠江口6 個臺站的實測值吻合較好。伶仃洋風暴潮增水的時間過程及空間分布表明，伶仃洋的風暴增水過程受風向、岸線幾何形態和局部水深地形的影響。伶仃洋的12 個特征點的水位過程線表明，風暴潮的時間過程曲線可分為初振段、主振段和余振段，灣內各點的風暴增水表現出比灣口附近和灣外更強的短期振動，表明灣內天文潮與風暴潮相互作用比灣口和灣外更強。

文章研究了臺風的登陸時間、中心氣壓差、登陸位置和移動方向等特征參數對風暴潮增水的影響。結果表明，不同登陸時間變化通過天文潮與風暴潮的相互作用導致了風暴增水的變化，且在灣內變化更為明顯，符合低潮位時風暴增水更大的規律。中心最低氣壓的變化影響增水的峰值大小，當氣壓差增大時，可以同時增大增水和減水的峰值，但對風暴潮增水的過程線影響較小。臺風登陸位置對伶仃洋的風暴增水的時空過程分布有很大影響，當臺風在伶仃洋西側海岸登陸，且風向近似垂直于海岸時，可能會出現極端的風暴增水。對于西北型登陸的臺風，當臺風登陸方向與海岸夾角減小時，會增大風暴增水的峰值，反之，增大角度則會削弱增水峰值。本文的研究結果有助于對珠江口伶仃洋及其周邊海區典型的西北向移動臺風風暴潮過程的認識，可為沿海防臺、防潮減災政策制定提供參考。