羅智豐，陳 剛，道付海

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

廣州市網(wǎng)河區(qū)大小水道、河涌縱橫交錯(cuò)，水網(wǎng)密布，珠江三角洲出海八大口門(mén)占其三：虎門(mén)、蕉門(mén)、洪奇瀝。城區(qū)珠江主河段均為感潮河段，據(jù)統(tǒng)計(jì)，影響廣州的臺(tái)風(fēng)平均每年2.5個(gè)，均給廣州帶來(lái)不同程度的風(fēng)暴潮災(zāi)害[1]。

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

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

因此，構(gòu)建廣州市風(fēng)暴潮精細(xì)化預(yù)報(bào)模型，形成面對(duì)對(duì)象的風(fēng)暴潮精細(xì)化預(yù)報(bào)，提高風(fēng)暴潮實(shí)時(shí)觀測(cè)與業(yè)務(wù)化預(yù)警報(bào)能力，對(duì)摸清廣州市風(fēng)暴潮特性、重點(diǎn)保障目標(biāo)排查，大力提升風(fēng)暴潮災(zāi)害綜合防御能力具有重要意義。

1 ADCIRC模型介紹

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

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

1.1 控制方程

采用垂向平均的二維方式計(jì)算[8]，通過(guò)笛卡爾坐標(biāo)下的連續(xù)方程和動(dòng)量方程求解自由表面起伏、二維流速等變量，其二維連續(xù)方程為：

(1)

二維動(dòng)量方程為：

(2)

(3)

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

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

(4)

其中Ax、Ay的表達(dá)式如下：

(5)

(6)

根據(jù)連續(xù)方程和動(dòng)量方程的表達(dá)，將Ax、Ay表達(dá)為:

(7)

(8)

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

1.2 底摩擦項(xiàng)

在ADCIRC中，底部應(yīng)力的表達(dá)式為：

τbx≡Uτ*

(9)

τbx≡Uτ*

(10)

其結(jié)果依賴于τ*的表達(dá)形式。τ*可以采用線性、二次函數(shù)或者線性和二次函數(shù)的復(fù)合公式見(jiàn)式(11)—(15)。摩擦系數(shù)Cf則用曼寧系數(shù)轉(zhuǎn)化而來(lái)，轉(zhuǎn)化關(guān)系見(jiàn)式(15)。

τ*=Cf

(11)

(12)

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

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

(13)

曼寧系數(shù)與摩擦系數(shù)轉(zhuǎn)化關(guān)系：

(14)

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

2 廣州市風(fēng)暴潮精細(xì)化模型搭建

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

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

外海開(kāi)邊界天文潮位由8個(gè)主要的調(diào)和分潮(K1、M2、N2、O1、S2、K2、P1、Q1)計(jì)算，見(jiàn)式(15)，風(fēng)場(chǎng)作為下邊界條件輸入該模型進(jìn)行計(jì)算。

(15)

式中ζ——總水位；a0——平均海面；f——節(jié)點(diǎn)因子；H——分潮振幅；V0+μ——分潮的天文初相角；K——分潮遲角；m——分潮個(gè)數(shù)；i——各個(gè)分潮。

白新歡提出馬克思共產(chǎn)主義思想具有科學(xué)、現(xiàn)實(shí)和哲學(xué)三個(gè)基本維度，其中前兩個(gè)維度使共產(chǎn)主義成為科學(xué)。他同時(shí)強(qiáng)調(diào)了馬克思人道主義思想是共產(chǎn)主義思想的重要出發(fā)點(diǎn)。[12]

模型同時(shí)考慮上游來(lái)水、天文潮、風(fēng)暴潮的共同作用，既能充分考慮多種因素的共同作用又能提高計(jì)算效率，考慮到臺(tái)風(fēng)預(yù)報(bào)信息的精度時(shí)效及模型的計(jì)算效率，目前模式可實(shí)現(xiàn)廣州市未來(lái)3 d的精細(xì)化預(yù)報(bào)，模型每次運(yùn)算時(shí)間約20 min。模型框架結(jié)構(gòu)見(jiàn)圖1。

圖1 模型結(jié)構(gòu)

2.1 計(jì)算網(wǎng)格

模型采用非結(jié)構(gòu)的三角形網(wǎng)格，可精確的模擬河口海岸地區(qū)的地形變化，且可方便地對(duì)所關(guān)注的區(qū)域進(jìn)行加密，網(wǎng)格步長(zhǎng)從外海到近岸逐漸變小，可有效地提高計(jì)算效率。為了更好地模擬風(fēng)暴潮在南海的生成與傳播過(guò)程，模式的計(jì)算范圍選擇為105.6°E～127°E和13°N～27°N，見(jiàn)圖4a(圖中顏色表示水深)。因?yàn)橐趶V州市網(wǎng)河區(qū)實(shí)現(xiàn)風(fēng)暴潮精細(xì)化預(yù)報(bào)，對(duì)于該區(qū)域內(nèi)的河段提高分辨率，見(jiàn)圖4b，網(wǎng)格最小分辨率為50 m。

a)海區(qū)網(wǎng)格

2.2 臺(tái)風(fēng)風(fēng)場(chǎng)

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

(16)

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

Jelesnianski臺(tái)風(fēng)風(fēng)場(chǎng)模型是一種經(jīng)驗(yàn)?zāi)Ｐ?，直接從風(fēng)場(chǎng)外觀相似性出發(fā)建立模型風(fēng)場(chǎng)，可表示為圓形風(fēng)場(chǎng)與移動(dòng)風(fēng)場(chǎng)的和：

V=Vr+Vs

(17)

其中的圓形風(fēng)場(chǎng)Vr和移動(dòng)Vs的表達(dá)式分別為：

(18)

(19)

式中Vm——臺(tái)風(fēng)最大風(fēng)速；Vc——臺(tái)風(fēng)中心移動(dòng)速度。

風(fēng)拖曳系數(shù)Cd采用Garratt公式：

(20)

3 天文潮驗(yàn)證

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

a)香港

c)南沙

表2 各站計(jì)算天文潮與潮汐表誤差比對(duì)

4 歷史臺(tái)風(fēng)暴潮驗(yàn)證

選取刷新廣州市多個(gè)潮位站歷史紀(jì)錄的典型臺(tái)風(fēng)2017年13號(hào)臺(tái)風(fēng)“天鴿”和2018年22號(hào)臺(tái)風(fēng)“山竹”(圖4)，對(duì)廣州市風(fēng)暴潮精細(xì)化模型進(jìn)行了驗(yàn)證。

圖4 “天鴿” “山竹”發(fā)展路徑

臺(tái)風(fēng)“天鴿”于2017年8月23日12時(shí)50分在珠海市金灣區(qū)沿岸以強(qiáng)臺(tái)風(fēng)級(jí)別登陸，登陸時(shí)中心最大風(fēng)力14級(jí)(45 m/s)，給珠三角地區(qū)帶來(lái)嚴(yán)重的風(fēng)暴潮增水，三角洲多個(gè)潮位站出現(xiàn)超百年一遇且破歷史極值的高潮位，中大站8月23日15 時(shí)45分出現(xiàn)2.76 m的高潮位，超當(dāng)時(shí)歷史極值(200814“黑格比”)0.03 m，超警戒1.26 m。1822號(hào)臺(tái)風(fēng)“山竹”于2018年9月16日17時(shí)在廣東江門(mén)臺(tái)山市海宴鎮(zhèn)登陸，登陸時(shí)為強(qiáng)臺(tái)風(fēng)級(jí)，中心附近最大風(fēng)力14級(jí)，相當(dāng)于45 m/s?！吧街瘛苯o珠江三角洲地區(qū)帶來(lái)了 2.60～3.00 m的風(fēng)暴潮增水，三角洲地區(qū)多個(gè)潮位站出現(xiàn)了突破歷史記錄極值的超百年一遇高潮位，中大站9月16日19 時(shí)35分出現(xiàn)3.23 m的高潮位，超歷史極值(201713“天鴿”)0.47 m，超警戒1.73 m。

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

a)南沙

a)南沙

表3 1713號(hào)“天鴿”影響期間各站風(fēng)暴潮模擬誤差統(tǒng)計(jì)

表4 1822號(hào)“山竹”影響期間各站風(fēng)暴潮模擬誤差統(tǒng)計(jì)

5 結(jié)語(yǔ)

廣州市網(wǎng)河區(qū)水位影響因素復(fù)雜，建立同時(shí)考慮上游來(lái)水、天文潮、風(fēng)暴潮共同影響的風(fēng)暴潮精細(xì)化預(yù)報(bào)模型，并以1713號(hào)“天鴿”、1822號(hào)“山竹”臺(tái)風(fēng)進(jìn)行驗(yàn)證，結(jié)果如下。

a)基于ADCIRC的風(fēng)暴潮精細(xì)化預(yù)報(bào)模型能較好的反映廣州市網(wǎng)河區(qū)的天文潮狀況，可實(shí)現(xiàn)為無(wú)實(shí)測(cè)數(shù)據(jù)、無(wú)潮位站點(diǎn)的河道網(wǎng)格化的天文潮預(yù)報(bào)。

b)該模型能較好地模擬上游來(lái)水、天文潮、風(fēng)暴潮共同影響的潮區(qū)河道總水位，并能較為準(zhǔn)確地對(duì)廣州市網(wǎng)河區(qū)風(fēng)暴潮過(guò)程進(jìn)行模擬，可為未來(lái)廣州市風(fēng)暴潮精細(xì)化預(yù)警預(yù)報(bào)工作提供參考。