羅智豐，陳 剛，道付海

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

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

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

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

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

1 ADCIRC模型介紹

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

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

1.1 控制方程

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

(1)

二維動量方程為：

(2)

(3)

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

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

(4)

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

(5)

(6)

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

(7)

(8)

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

1.2 底摩擦項

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

τbx≡Uτ*

(9)

τbx≡Uτ*

(10)

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

τ*=Cf

(11)

(12)

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

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

(13)

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

(14)

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

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

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

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

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

(15)

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

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

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

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

2.1 計算網(wǎng)格

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

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

2.2 臺風風場

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

(16)

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

Jelesnianski臺風風場模型是一種經(jīng)驗?zāi)Ｐ停苯訌娘L場外觀相似性出發(fā)建立模型風場，可表示為圓形風場與移動風場的和：

V=Vr+Vs

(17)

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

(18)

(19)

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

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

(20)

3 天文潮驗證

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

a)香港

c)南沙

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

4 歷史臺風暴潮驗證

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

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

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

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

a)南沙

a)南沙

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

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

5 結(jié)語

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

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

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