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

林澤輝

(惠州市惠陽區永良堤圍管理所，廣東 惠州 516007)

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

1 研究方法

1.1 區域概況

深圳市海岸線長約250km，東起大亞灣，西至珠江口，增水>50cm的風暴潮過程年平均出現2.1次。沿岸分布有機場、重點港口、大型住宅區、核電站等大型企業，人口密度極高。

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

1.2 數據資料

影響風暴增水量級的因素有很多，包括臺風氣壓場分布、臺風風場的方向與大小以及近岸地形變化因素。在從我國東南沿海登陸的臺風里，影響到陽江至汕尾海域的臺風均可引起深圳海域的風暴潮增水。目前在深圳海域，國家海洋局共設有兩個長期驗潮站，分別是赤灣和鹽田站，緊鄰深圳大亞灣岸線的有位于惠州市的惠州海洋站，各站位置見圖1。收集上述站點自建站以來臺風期間的潮位資料統計歷史風暴潮增水情況。廣東省水文局1964年在深圳市赤灣港內建設有長期水文觀測站，1986年國家海洋局南海分局在相距不到1km處建設了赤灣海洋站，并于1986年開始投入業務化觀測。鹽田海洋站位于大鵬灣鹽田港海域，由國家海洋局南海分局建設和管理，2003年投入使用。惠州海洋站位于大亞灣灣頂西側荃灣港區，由國家海洋局南海分局建設和管理，2006年建成投入使用，他們的分布位置如圖1所示。

圖1 資料站位位置示意圖

1.3 研究方法

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

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

1.4 模型設置

1.4.1 模型介紹

文章采用的是DELFT3D-FLOW與SWAN的耦合模型，DELFT3D是一款由荷蘭DELTRA開發的開源水動力模擬軟件，在風暴潮模擬方面已經廣泛應用。

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

質量守恒方程：

(1)

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

動量守恒方程：

(2)

(3)

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

SWAN模型的控制方程為：

(4)

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

1.4.2 模型設置

1)模型范圍：

文章所采用的模型囊括整個珠江三角洲區域，見圖2。

圖2 模型區域示意圖

2)地形水深：

模擬范圍內的地形水深采用的是中華人民共和國海事局2016年繪制的海圖插值數據。將地形水深數據轉化到統一基面(平均海平面)后，用IDW插值方法插值到網格點上，得到模型的水深地形。

3)開邊界設置：

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

4)底摩擦設置：

模型根據水深范圍設置不同大小的糙率，水深越大糙率越小。

1.4.3 風場輸入

文章的風場采用美國國家大氣研究中心(NCAR)等機構共同開發的中尺度大氣預報系統。WRF模式適用于臺風模擬，可以通過嵌套網格運算來達到高精度模擬的目的。本次采用了WSM6微物理方案、KainFritsch卷積云方案、Monin-Obukhov表面層方案、Noah陸面層方案、MYJ大氣邊界層方案、RRTMG輻射方案模擬臺風。

2 模型驗證

2.1 風場驗證

首先利用實測資料對所建立的風暴潮-波浪耦合模型進行率定驗證。選取1713號臺風(2017年第13號臺風—天鴿，其移行路徑如圖3所示)作為模型邊界條件，為了證明WRF模式的準確性其在風暴潮模擬方面的優勢，文章也用其他論文常用的經驗風場模型做了驗證對比。經驗風場采用的是宮崎正衛模型作為移動風場，HOLLAND模型作為梯度風場，他們的公式分別為：

圖3 “天鴿”臺風移行路徑

(5)

(5)

式中：B為HOLLAND模型氣壓參數；f為科氏力參數，與模擬區域緯度有關；R為最大風速半徑，與中心氣壓有關。

圖4為WRF模型模擬風場值以及經驗風場模擬值與實測風場的對比，可以看出，WRF模型模擬的風場在峰值以及過程上都比經驗風場模型更貼近實際，因此利用WRF模型生成的臺風風場作為風暴潮-波浪耦合模型的邊界條件更為合理。

圖4 WRF風場及模型風場與實測風場對比

2.2 風暴潮及波浪率定驗證

以1713號臺風(2017年第13號臺風)為邊界條件分別驅動珠江口風暴潮模型以及珠江口波浪-風暴潮耦合模型，得到如圖5所示的赤灣站與鹽田站水位驗證曲線，可以看出，考慮波浪后的水位過程曲線更符合實際，最大誤差絕對值為13cm，平均相對誤差為11.3%，標準差為2.94，符合模擬精度要求。

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

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

3 深圳市沿岸漫堤風險評估

3.1 可能最大風暴潮模型計算

為了尋求可能影響深圳市的最大風暴潮對深圳市沿岸的影響，采用概率論法進行計算可能最大臺風中的近中心最大風速Vmax，臺風中心氣壓P0。

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

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

圖6(b) 影響深圳的臺風中心最大風速極值I型分布

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

輸出堤前代表點(具體見圖7)的最高潮位與堤頂高程進行比對，了解沿岸淹沒的情況，由于模式計算結果是相對平均海平面的，為便于比對，需要根據當地平均海平面與85高程的關系，將計算結果轉成以85黃海基面為起算面，其中西段提防采用赤灣站高程關系，即計算結果加上0.54m；東段提防采用鹽田站高程關系進行轉換，即計算結果加上0.66m。表1列出堤前統計的最高潮位及堤頂高程。從代表點的最高潮位和海堤高程比較來看，在可能最大熱帶氣旋影響下，深圳沿岸均可能發生漫堤。

圖7 堤前代表點位置示意圖

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

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

3.2 海堤防御標準

對于大部分沿海地區來說，1000a一遇的標準是過高的。但是深圳作為特大型城市，其風暴潮安全具有重要意義，因此考慮1000a一遇情形下的風暴潮安全是有必要的，在此情形下得出的淹沒情況有利于應急管理部門采取相應的防范措施。

4 結 論

1)利用WRF中長尺度氣象預報模式作為波浪-風暴潮耦合模型的風場邊界條件輸入，提高了風暴潮模型的精確度

2)在珠江口波浪-風暴潮模型驗證實測資料準確的基礎上，利用極值I型分布確定1000a一遇臺風參數對影響深圳市海岸堤防未來發生的最大可能風暴潮進行計算分析。文章研究成果可為深圳市風暴潮安全防范提供重要技術支撐，同時也為海洋水利防災應急部門提供決策參考。