基于ADCIRC的廣州市風暴潮精細化預報模型的建立與驗證

羅智豐，陳 剛，道付海

(1.廣東省水文局廣州水文分局，廣東 廣州 510150；2.廣東海啟星海洋科技有限公司，廣東 廣州 511400)

廣州市網河區大小水道、河涌縱橫交錯，水網密布，珠江三角洲出海八大口門占其三：虎門、蕉門、洪奇瀝。城區珠江主河段均為感潮河段，據統計，影響廣州的臺風平均每年2.5個，均給廣州帶來不同程度的風暴潮災害[1]。

臺風具有發生頻率高、受災范圍廣、災害強度大以及次生災害多等特點[2]，近年來廣州風暴潮災害呈愈發嚴重的趨勢，以廣州中心城區代表站——位于海珠區中山大學北門珠江側的中大站為例，該站歷史最高水位排名前6名中，除1次受洪水影響，其余均為風暴潮增水，且均出自2000年之后，0104號“尤特”、0814號“黑格比”、1713號“天鴿”、1822號“山竹”等臺風均刷新歷史紀錄(表1)，尤其“天鴿”“山竹”給廣州市帶來嚴重的災害損失。據統計，“山竹”導致廣州南沙、番禺、黃埔、海珠、天河、越秀、荔灣、白云等區共68處堤防浸水，黃埔區江水越堤水浸達15 km，越秀區珠江前航道二沙島區域江水越堤水浸達7.6 km，受淹車庫54個(數據來源：廣州市三防辦)。

表1 中大站實測最高水位排名(珠江基面) 單位：m

因此，構建廣州市風暴潮精細化預報模型，形成面對對象的風暴潮精細化預報，提高風暴潮實時觀測與業務化預警報能力，對摸清廣州市風暴潮特性、重點保障目標排查，大力提升風暴潮災害綜合防御能力具有重要意義。

1 ADCIRC模型介紹

ADCIRC(An Advanced Circulation Model For Oceanic,Coastal And Estuarine Waters)是由美國北卡羅來納大學的 R.A.LUETTICH和美國圣母大學的J.JWESTERINK 于1992年研制的基于有限元方法，可用于海洋、海岸、河口區域的水動力模式[3]。被美國工程兵部隊(ACE)和美國海軍研究所(NRL)廣泛應用于各個軍港的潮汐、海流和風暴潮的預報中，ADCIRC也被NOAA應用于美國東海岸的風暴潮預報及路易斯安那州的風暴潮風險評估中。近年來，國內許多學者也將其應用于中國沿海的風暴潮數值模擬，劉克強等[4]構建了適用于西北太平洋海域的風暴潮預報模型，李歡等[5]將其應用于寧波沿海風暴潮預報，陳海軍等[6]建立了渤海高分辨率的二維潮汐潮流模型，羅鋒等[7]建立了江蘇海域的精細化風暴潮數值預報模型。

ADCIRC有以下特征：采用廣義波動連續方程(GWCE)與動量方程結合，基于伽遼金有限元方法求解；可應用笛卡爾坐標或球坐標、進行二維或三維的計算，并采用并行計算提高計算效率；采用不規則的三角形網格，可精確的模擬河口海岸地區的地形，且可方便地對所關注的區域進行加密，在非關注區域則可以設置較大間距的網格，以節省計算量；采用干濕邊界處理技術；可用于計算河堤、海堤等障礙物的溢流計算。

1.1 控制方程

采用垂向平均的二維方式計算[8]，通過笛卡爾坐標下的連續方程和動量方程求解自由表面起伏、二維流速等變量，其二維連續方程為：

(1)

二維動量方程為：

(2)

(3)

式中ζ——從海平面起算的水位高度；U、V——2個方向的垂向平均流速；H——總水深；f——科氏參量；Ps——水面大氣壓力；ρ0——水密度；g——重力加速度；τbx、τby——底部切應力；τsx、τsy——風應力項；Tx、Ty——波浪輻射應力項；Dx、Dy——擴散項。

為避免和減小伽留金法離散出現的虛假振蕩，該模式中連續方程的形式進行了一些處理。將式(1)對時間求偏導數，再將其與連續方程和一個隨空間變化的加權函數τ0的乘積相加，整理得到：

(4)

其中Ax、Ay的表達式如下：

(5)

(6)

根據連續方程和動量方程的表達，將Ax、Ay表達為:

(7)

(8)

運用式(7)、(8)中的表達式對式(4)當中的Ax、Ay進行替換，就得到了模式計算中所運用的GWCE(Generalized Wave Continuity Equation)方程[9]。ADCIRC通過對式(2)、(3)、(7)、(8)的聯立求解，獲得所需的水位和速度值。在球坐標系下進行計算時，利用CPP(Carte Parallelo-grammatique)圓柱法，將坐標投影到笛卡爾坐標系中，然后進行計算。

1.2 底摩擦項

在ADCIRC中，底部應力的表達式為：

τbx≡Uτ*

(9)

τbx≡Uτ*

(10)

其結果依賴于τ*的表達形式。τ*可以采用線性、二次函數或者線性和二次函數的復合公式見式(11)—(15)。摩擦系數Cf則用曼寧系數轉化而來，轉化關系見式(15)。

τ*=Cf

(11)

(12)

式(11)、(12)中U、V為深度平均速度。

復合形式公式同式(12)，但其中Cf隨Hc變化：

(13)

曼寧系數與摩擦系數轉化關系：

(14)

式中Hc——臨界水深，反映底摩擦在淺水區域增大的一個參數；θ——決定底摩擦系數與漸近線接近的快慢；λ——決定從深水到淺水底摩擦增大的快慢；n——曼寧系數；d——相對于海平面的水深。

2 廣州市風暴潮精細化模型搭建

廣州市網河區為感潮河段，水位受上游西北江、東江、本地區流溪河、增江來水影響大，因此要準確的模擬該區域河流水動力必須考慮上游來流的影響。分別選取西江的高要、北江的石角、東江的博羅、流溪河的太平場、增江的麒麟咀水位站所在位置為模型的上邊界起算點(圖 2b)，并以這5個水文站的實測流量作為模型上邊界驅動條件，把各個實測流量平均到對應河流上邊界的每個網格節點上，以本質邊界條件輸入模式，即通過在連續方程中指定對法向邊界通量積分的非零值和在動量方程中指定非零法向速度來實現的，該邊界條件滿足全局和每個邊界節點處的通量平衡。

風暴潮與潮汐并非是線性的疊加關系，它們存在非線性相互作用，且較嚴重的風暴潮災害往往發生在風暴潮遇到潮汐高潮位時，兩者的共同作用推高水位，造成危害，因此，對風暴潮與潮汐進行耦合計算，可以使計算結果更加準確。模式的潮汐計算準確與否很大程度上取決于外海開邊界上潮位的準確性，很多全球大洋潮汐模式均可通過插值的方式來提供大洋中每個點的潮汐資料比如TPXO7.2、NAO.99b、GOT00.2、DTU10等模式[10]，而鄭后俊等[11]采用293個站點的觀測資料對各個模式在南海的潮汐預報能力進行分析比較，發現DTU10模式的整體預報效果最好，高秀敏等[12]采用南海海域60個驗潮站和22個TOPEX/Poseidon衛星高度計軌道交叉點的調和常數資料，對比了4種模式4個主要分潮調和常數在南海的準確度，表明南海北部和東部區域DTU10的準確度最高。DTU10是丹麥科技大學在2010年建立的全球大洋潮汐模式，它是以FES2004為參考，基于響應法對T/P和Jason-1/2衛星高度計資料進行殘差分析而建立的，分辨率為0.125°[13-14]，本文的開邊界潮汐資料選用DTU10。南海區域的潮波主要是從太平洋經呂宋海峽傳入，本文所建立的模式，計算區域較大，涵蓋呂宋海峽和南海中、北部區域，可以更好地模擬超波的傳播過程。

外海開邊界天文潮位由8個主要的調和分潮(K1、M2、N2、O1、S2、K2、P1、Q1)計算，見式(15)，風場作為下邊界條件輸入該模型進行計算。

(15)

式中ζ——總水位；a0——平均海面；f——節點因子；H——分潮振幅；V0+μ——分潮的天文初相角；K——分潮遲角；m——分潮個數；i——各個分潮。

白新歡提出馬克思共產主義思想具有科學、現實和哲學三個基本維度，其中前兩個維度使共產主義成為科學。他同時強調了馬克思人道主義思想是共產主義思想的重要出發點。[12]

模型同時考慮上游來水、天文潮、風暴潮的共同作用，既能充分考慮多種因素的共同作用又能提高計算效率，考慮到臺風預報信息的精度時效及模型的計算效率，目前模式可實現廣州市未來3 d的精細化預報，模型每次運算時間約20 min。模型框架結構見圖1。

圖1 模型結構

2.1 計算網格

模型采用非結構的三角形網格，可精確的模擬河口海岸地區的地形變化，且可方便地對所關注的區域進行加密，網格步長從外海到近岸逐漸變小，可有效地提高計算效率。為了更好地模擬風暴潮在南海的生成與傳播過程，模式的計算范圍選擇為105.6°E～127°E和13°N～27°N，見圖4a(圖中顏色表示水深)。因為要在廣州市網河區實現風暴潮精細化預報，對于該區域內的河段提高分辨率，見圖4b，網格最小分辨率為50 m。

a)海區網格

2.2 臺風風場

風暴潮模式結果的準確性在很大程度上依賴于風場模式的質量。本文采用常用的Jelesnianski臺風氣壓模型[15]和風場模型[16]，氣壓場分布是由壓差(Pn-Pc)以及最大風速半徑Rmax決定的，其表達式為：

(16)

式中P(r)——距離臺風中心r的海表面氣壓值；Pc——臺風中心氣壓；Pn——臺風以外不受干擾的背景氣壓；Rmax——臺風最大風速半徑；r——距離臺風中心的距離。

Jelesnianski臺風風場模型是一種經驗模型，直接從風場外觀相似性出發建立模型風場，可表示為圓形風場與移動風場的和：

V=Vr+Vs

(17)

其中的圓形風場Vr和移動Vs的表達式分別為：

(18)

(19)

式中Vm——臺風最大風速；Vc——臺風中心移動速度。

風拖曳系數Cd采用Garratt公式：

(20)

3 天文潮驗證

天文潮驗證計算的時間段為2019年8月27日到9月30日23時，共計算35 d，前5 d為模式冷啟動調整時間，不進行分析。模式結果與國家海洋信息中心發布的潮汐表天文潮結果做比對。選取了珠江口附近城市的6個代表潮位站進行比對，分別為香港、赤灣(深圳)、南沙(廣州)、澳門、珠海港、臺山(江門)等潮位站。圖3為各潮位站計算與發布的天文潮對比曲線。結果顯示：各站天文潮計算結果與潮汐表所示潮位的平均絕對誤差，最大平均絕對誤差小于13 cm(表2)，模型計算結果與潮汐表的天文潮結果基本相符，表明該模型基本可以反映出珠江口沿岸天文潮狀況。

a)香港

c)南沙

表2 各站計算天文潮與潮汐表誤差比對

4 歷史臺風暴潮驗證

選取刷新廣州市多個潮位站歷史紀錄的典型臺風2017年13號臺風“天鴿”和2018年22號臺風“山竹”(圖4)，對廣州市風暴潮精細化模型進行了驗證。

圖4 “天鴿” “山竹”發展路徑

臺風“天鴿”于2017年8月23日12時50分在珠海市金灣區沿岸以強臺風級別登陸，登陸時中心最大風力14級(45 m/s)，給珠三角地區帶來嚴重的風暴潮增水，三角洲多個潮位站出現超百年一遇且破歷史極值的高潮位，中大站8月23日15 時45分出現2.76 m的高潮位，超當時歷史極值(200814“黑格比”)0.03 m，超警戒1.26 m。1822號臺風“山竹”于2018年9月16日17時在廣東江門臺山市海宴鎮登陸，登陸時為強臺風級，中心附近最大風力14級，相當于45 m/s。“山竹”給珠江三角洲地區帶來了 2.60～3.00 m的風暴潮增水，三角洲地區多個潮位站出現了突破歷史記錄極值的超百年一遇高潮位，中大站9月16日19 時35分出現3.23 m的高潮位，超歷史極值(201713“天鴿”)0.47 m，超警戒1.73 m。

選取廣州市南沙區南沙站，番禺區三沙口、三善滘站、大石站，黃埔區黃埔站及海珠區中大站等6個潮位站，對2場強臺風的風暴潮進行模擬驗證，結果顯示模型對各站風暴潮的增水過程、最高水位峰值和出現時間均模擬得較好(圖5、6)：“天鴿”臺風各站模型計算最高水位誤差在12～39 cm，最大相對誤差個別站點12.5%，大部分站在10%以內，相位最大誤差均在20 min以內；“山竹”臺風各站模型計算最高水位誤差在0～31 cm，模型相對誤差均在10%以內，相位最大誤差除三善滘站45 min附近，其余各站誤差均不超20 min。具體誤差統計見表3、4。

a)南沙

a)南沙

表3 1713號“天鴿”影響期間各站風暴潮模擬誤差統計

表4 1822號“山竹”影響期間各站風暴潮模擬誤差統計

5 結語

廣州市網河區水位影響因素復雜，建立同時考慮上游來水、天文潮、風暴潮共同影響的風暴潮精細化預報模型，并以1713號“天鴿”、1822號“山竹”臺風進行驗證，結果如下。

a)基于ADCIRC的風暴潮精細化預報模型能較好的反映廣州市網河區的天文潮狀況，可實現為無實測數據、無潮位站點的河道網格化的天文潮預報。

b)該模型能較好地模擬上游來水、天文潮、風暴潮共同影響的潮區河道總水位，并能較為準確地對廣州市網河區風暴潮過程進行模擬，可為未來廣州市風暴潮精細化預警預報工作提供參考。