粵港澳大灣區風暴潮數值模型的建立與應用

羅志發，黃本勝，譚 超，黃廣靈

(1.廣東省水利水電科學研究院，廣東 廣州 510635；2.廣東省水動力學應用研究重點實驗室，廣東 廣州 510635；3. 廣東省流域水環境治理與水生態修復重點實驗室，廣東 廣州 510635；4. 中山大學，廣東 廣州 510275)

粵港澳大灣區是我國開放程度最高、經濟活力最強的區域之一，以不足全國1%的土地面積和不足全國5%的人口，貢獻了全國經濟總量的17%，在國家發展大局中具有重要戰略地位。粵港澳大灣區位于中國大陸南端，瀕臨南海，地處珠江流域下游，河網水系發達。由于特殊的地理位置和氣候，導致大灣區臺風暴潮災害易發頻發，由此造成的人員傷亡、經濟損失相當巨大，已經成為影響人民生活質量、制約國民經濟高質量發展的重要因素。開展粵港澳大灣區風暴潮研究既是重要科學問題也是風暴潮預警預報關鍵技術研發的基礎，對提升大灣區水安全保障能力具有重要的意義。

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

1 模型的建立

1.1 臺風場的構建

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

臺風風場模型根據梯度風原理，由臺風氣壓場計算出風場。本文采用國內外應用較廣泛的Fujita[4]、Ueno[5]公式構建氣壓場和臺風風場，其式如下：

(1)

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

(2)

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

(3)

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

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

(4)

1.2 風暴潮模型的建立

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

圖1 模型嵌套計算網格示意

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

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

2 模型驗證

2.1 臺風場驗證

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

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

圖3 澳門氣象站“山竹”臺風過程風速風向驗證結果示意

計算的臺風風速、風向其變化過程與實測數據相符合，最大風速絕對誤差僅1 m/s。說明計算的風場較好地反映了臺風過程及最大風速。圖4為采用本文臺風場構建方法計算的臺風場圖，由圖4可知，在經驗模型的基礎上融合了ERA50的再分析數據，對臺風場的計算避免了經驗公式關于臺風外圍風速計算明顯偏小的問題，構建的風場結構較為合理，可用于風暴潮數值模擬。

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

2.2 天文潮驗證

模擬風暴潮和天文潮耦合的水位過程時，天文潮的模擬計算是基礎。在不考慮風應力的情況下，對南海天文潮進行數值模擬，以檢驗模型對天文潮計算的可靠性。計算時間為2018年7月1—31日，選用沿海大萬山、珠海港、上川島、馬鞭洲等潮位站(站位分布見圖2)實測資料對模型計算結果進行驗證，比較結果見圖5所示。

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

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

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

2.3 風暴潮水位過程驗證

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

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

圖6 “山竹”臺風過程水位驗證示意

表1 誤差統計

3 結果分析

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

3.1 風暴潮增減水特征分析

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

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

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

3.2 天文潮風暴潮非線性相互作用分析

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

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

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

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

4 結語

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

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