羅志發(fā)，黃本勝，譚 超，邱 靜，黃廣靈

(1.廣東省水利水電科學研究院，廣東 廣州 510635；2.廣東省水動力學應用研究重點實驗室，廣東 廣州 510635；3.廣東省流域水環(huán)境治理與水生態(tài)修復重點實驗室，廣東 廣州 510635)

1 概述

珠江河口位于中國大陸南端，瀕臨南海，海岸線狹長，沿岸地勢較低。特殊的地理位置及氣候條件使得珠江河口成為國內(nèi)臺風登陸最頻繁、風暴潮災害最嚴重的區(qū)域之一。近年來，在全球氣候變化背景下，風暴潮災害次數(shù)和強度呈增加的趨勢，登陸珠江河口的超強臺風引發(fā)的風暴潮災害，屢次刷新沿岸潮位站歷史最高潮位和最大增水記錄。珠江河口沿岸人口稠密、經(jīng)濟發(fā)達，每年由于風暴潮造成的經(jīng)濟損失相當巨大，已經(jīng)成為影響人民生活質(zhì)量、制約國民經(jīng)濟高質(zhì)量發(fā)展的重要因素。

為了防御風暴潮災害，珠江河口沿海陸續(xù)興建了江海堤圍。盡管如此，強臺風往往伴隨著強臺風浪，在強浪沖擊作用下，江海堤圍面臨著堤身結(jié)構(gòu)損毀的風險，可能導致海水漫過海堤入侵內(nèi)陸，引發(fā)城市內(nèi)澇。因此，準確模擬風暴潮及臺風浪過程，了解正面登陸臺風影響下珠江河口近岸海域風暴潮增水及臺風浪變化過程，對準確評估風暴潮的災害風險以及做好風暴潮災害防御與應對工作、減少風暴潮災害損失有著重要現(xiàn)實的意義。

珠江河口是多種動力因子協(xié)同作用的復雜系統(tǒng)，呈“三江匯流，八口入海”的形勢。河網(wǎng)區(qū)水網(wǎng)密布，橫向支汊發(fā)育，其間多種動力因子相互耦合，動力復雜多變。珠江河口風暴潮、波浪數(shù)值模擬已開展較多的研究，如二維風暴潮模型的建立[1]，風浪數(shù)值模擬[2-3]，路徑、風速對風暴潮增水的影響[4]，地形對局部增水的影響[5]等。波浪在近岸傳播過程中受到淺水變形的影響發(fā)生破碎，產(chǎn)生的輻射應力驅(qū)動水面抬升形成增水，在臺風期間波浪增水相對較大[6]。波浪與風暴潮耦合計算已開展了相關(guān)的研究，在不同的海域，由于岸線地形及水動力條件的不同，波浪對風暴潮增水的影響有所差異[7-9]。針對珠江河口復雜的動力系統(tǒng)，本文基于無結(jié)構(gòu)三角網(wǎng)格的數(shù)值模式，優(yōu)化擬合復雜岸形，以珠江河網(wǎng)—河口作為整體，綜合考慮風、徑流、潮流、波浪等動力因子的耦合作用，構(gòu)建珠江河口波浪—風暴潮耦合數(shù)值模型，可有效提高風暴潮數(shù)值模擬的準確度，為進一步探討風暴潮—波浪耦合作用提高有效的科學手段。以1822號“山竹”臺風為例，對臺風過程的風暴潮增水及波浪進行數(shù)值模擬，定量分析考慮波浪—風暴潮耦合作用對風暴潮數(shù)值模擬結(jié)果的影響。

2 模型的建立與驗證

2.1 模型簡介

本文選用海洋環(huán)流模式SELFE建立南海—珠江河口雙重嵌套模型(模型網(wǎng)格見圖1)，該模式基于非結(jié)構(gòu)化網(wǎng)格，可精細化擬合復雜岸線和地形，采用半隱式的歐拉—拉格朗日有限元算法求解N-S方程組，可有效提高模型計算效率。南海模型計算范圍為98°E～126°E，0°N～30°N，涵蓋整個南海及西北太平洋海域，南邊界至卡里馬塔海峽、北邊界至浙江省沿岸海域、東邊界至48小時警戒線。網(wǎng)格分辨率從近岸的 1 km 逐漸過渡到外海20 km，水深數(shù)據(jù)采用ETOP01全球1′×1′分辨率的地形資料。珠江河口模型上邊界為西江高要，北江石角，東江博羅，流溪河及潭江上游，外海下邊界取約100 m等深線處。模型采用無結(jié)構(gòu)網(wǎng)格，擬合復雜河岸邊界，對局部進行加密，提高網(wǎng)格分辨率。模型網(wǎng)格共有101 752個節(jié)點，173 045個網(wǎng)格單元，網(wǎng)格大小為從河網(wǎng)區(qū)10 m逐漸過渡到外海的20 km；在垂向上采用Sigma坐標，均勻分為10層；珠江三角洲網(wǎng)河區(qū)采用2005—2008年的大范圍實測地形，河口區(qū)及近岸海區(qū)采用2000—2008年海圖地形。南海模型由風場、氣壓場及8個主要分潮(M2、S2、N2、K2、K1、O1、P1、Q1)[10]進行驅(qū)動，計算得到余水位，將余水位以及風場、氣壓場、8個分潮作為珠江口模型的驅(qū)動條件，對珠江河口風暴潮增水進行數(shù)值模擬，模型詳細介紹見文獻[7]。

圖1 模型嵌套計算網(wǎng)格示意

在風暴潮模型的基礎(chǔ)上耦合波浪計算模塊，建立波浪—風暴潮耦合模型。波浪模塊運用Wind Wave Model-III(WWM-III)模式[11]，該模式主要用于模擬大洋—近岸波浪傳播過程，模擬的物理過程包括：波浪折射作用，淺水變形作用，風成浪作用，白浪的耗散作用，水深引起的破碎作用，海底摩擦作用等。控制方程為能量譜平衡方程：

(1)

式中N為波作用密度譜;Cgx和Cgy分別為x和y方向的波群傳播速度;u和v分別為x和y方向的水流速度;σ為相對波頻;θ為波向;Cσ和Cθ分別為σ和θ方向的波浪傳播速度;Stot為譜密度表示的源項，是各種物理現(xiàn)象的源函數(shù)的疊加形式。

風暴潮模型為波浪模型提供水位、流速、流向數(shù)據(jù)計算變化水動力條件下的波浪要素。波浪模型輸出有效波高、平均波向、有效波周期、波長和底部波軌流速，輸出結(jié)果作為邊界波浪條件輸入風暴潮模型，計算波流共同作用下的底切應力及垂向混合系數(shù)。

2.2 模型驗證

文獻[12]以“山竹”臺風為算例，對模型計算的風速、風向、天文潮位、風暴潮位等進行了率定驗證，驗證結(jié)果良好，本文在此基礎(chǔ)上對波浪計算結(jié)果進行進一步驗證。珠江口波浪模型外海波浪開邊界條件通過南海模型提供，海表面臺風場及氣壓場通過臺風模型及氣壓模型給出[12]，底摩擦系數(shù)設置為0.067[13]，波浪破波系數(shù)為0.78[11]。波浪計算時間步長設置為1 800 s，每個時間步長計算結(jié)束后與風暴潮模型進行數(shù)據(jù)反饋。模擬時間為2018年9月10日—9月20日，前5 d用于計算的穩(wěn)定，后5 d對風暴潮及波浪過程進行模擬。

圖2為外海浮標站(113.9947°E，21.4948°N)有效波高驗證示意，由圖2可知，計算結(jié)果有效波高變化過程與實測數(shù)據(jù)有效波高變化趨勢一致，相關(guān)系數(shù)為0.862，平均絕對誤差為0.53 m，二者吻合較好。計算有效波高的最大值與實測結(jié)果分別為9.4 m、9.9 m，相對誤差僅5.3%，表明模型計算結(jié)果能較好地反映有效波高的極值情況。由圖3可知，考慮波浪作用之后，潮位數(shù)值模擬結(jié)果更接近觀測值。橫門站考慮波浪作用后，最高潮位的絕對誤差由0.24 m減小為0.06 m；南沙站的模擬精度同樣有所提高。由此可見，考慮波浪—風暴潮耦合作用后可提高風暴潮數(shù)值預報的精度。綜上所述，模型計算結(jié)果能夠較好地反映“山竹”臺風期間波浪及風暴潮增水的變化過程，可為進一步分析波浪與風暴潮耦合提供基礎(chǔ)。

圖2 外海浮標站有效波高驗證示意

圖3 橫門(a)、南沙(b)站考慮和不考慮波浪作用潮位模擬結(jié)果示意

3 結(jié)果與討論

本節(jié)利用上文已建立的波浪—風暴潮耦合數(shù)值模型，模擬“山竹”臺風風暴潮增水及波浪過程，討論波浪—風暴潮耦合對風暴潮增水數(shù)值模擬的影響。

3.1 波浪及風暴潮增水計算結(jié)果分析

“山竹”臺風過程對珠江口近岸波浪要素分布影響顯著，2018年9月16日13：00，臺風中心位于珠江口南側(cè)(見圖4a)，臺風中心風力較小，有效波高較低，約為6～7 m；臺風中心右側(cè)風力較大，有效波高超過9 m。此時，伶仃洋海域為東北偏東風，珠江口有效波高相對較小，約為2～4 m。16日17:00，臺風于廣東臺山海宴鎮(zhèn)登陸，珠江口盛行東南偏南風，有效波高有顯著增大趨勢，主要高值區(qū)更偏向近岸，有效波高最大值介于4～5 m(見圖4b)。選取燈籠山、南沙、萬頃沙、橫門站位(站位位置見圖4a)分析臺風期間有效波高變化過程。由圖5可知，臺風過境期間，各站位有效波高呈現(xiàn)先增大后減小的過程。磨刀門往外延伸瀕臨南海，因此燈籠山站受波浪影響較大，最大有效波高達3.2 m。南沙、萬頃沙、橫門等站位均位于伶仃洋灣內(nèi)，外海波浪傳入伶仃洋經(jīng)底摩擦、折射、繞射等作用，波能消耗，有效波高均小于2 m。

圖5 各站位臺風期間有效波高變化

臺風登陸前4 h，伶仃洋海域為東北偏東風，伶仃洋東岸在離岸風的驅(qū)動下，產(chǎn)生減水，深圳灣區(qū)域最大減水值約為-1.5 m。伶仃洋西岸由于水體的橫向堆積，產(chǎn)生0.5～1.0 m的增水(見圖4c)。磨刀門、黃茅海區(qū)域為偏北風，有利于水體離岸輸運，形成-0.5～-1.0 m的減水。臺風登陸后(見圖4d)，珠江口海域普遍為東南偏南風，珠江河口東南向的開口方向有利于水體向岸堆積并沿河道向上游輸運，此時珠江口海域普遍達到增水的最大值。伶仃洋河口灣頂增水值較大，約為3.1 m，是由于伶仃洋河口灣喇叭狀的形態(tài)有利于水體向灣頂聚集。其余口門區(qū)域最大增水值普遍為2.6～3.1 m。

圖4 “山竹”臺風登陸前后波浪場及風暴潮增水分布(a、c：登陸前，b、d：登陸后)

3.2 波浪對風暴潮影響分析

利用風暴潮與波浪耦合計算的總水位減去未考慮波浪作用的風暴潮水位可得到波浪增水，波浪增水包含了由波浪輻射應力產(chǎn)生的增水和波流相互作用產(chǎn)生的增水。圖6為臺風登陸前后波浪增水分布示意。

臺風登陸前(見圖6a)波浪增水影響范圍較廣，波浪增水由外海向近岸逐漸增大，外海波浪增水約為10 cm，近岸海域波浪增水約為15～20 cm，與相關(guān)研究結(jié)果較為一致[9]。近岸海域中，伶仃洋、香港島附近海域波浪增水較大，黃茅海等海域波浪增水相對較小，其空間分布與有效波高空間分布大體一致。臺風登錄后(見圖6b)，波浪增水集中在近岸海域，主要分布在伶仃洋、磨刀門及黃茅海等海域近岸區(qū)域，波浪增水最大值約20 cm。

圖6 “山竹”臺風登陸前后波浪增水分布示意

選取泗盛圍、南沙、橫門、燈籠山等站位分析風暴潮增水及波浪增水過程，由圖7可知，各站位風暴潮增水呈單峰型，增水快速增大至最大值而后迅速下降，最大增水值均超過3 m。波浪增水同樣呈現(xiàn)先增大后減小的過程，最大波浪增水提前于最大風暴潮增水。各站位最大波浪增水均超過0.2 m，約占風暴潮增水的6%。

圖7 各站位風暴潮增水與波浪增水過程示意

波浪增水占風暴潮增水的比例雖然不大，但對于設計潮位具有一定的影響。表1列舉了珠江口潮位站不同頻率下的設計潮位值，由表1可知，同一測站其不同頻率的設計潮位僅相差20～30 cm，如西炮臺、南沙、燈籠山200年一遇潮位與100年一遇設計潮位值分別相差0.17、0.25、0.22 cm，與最大波浪增水值相當。考慮波浪—風暴潮耦合計算的方法，對推算沿岸設計潮位具有一定的參考意義。

表1 各測站潮位設計值

4 結(jié)語

1)本文在已建立風暴潮數(shù)值模型的基礎(chǔ)上，耦合波浪模塊，建立了珠江河口波浪—風暴潮耦合數(shù)值模型。采用實測波浪資料對“山竹”臺風波浪過程進行模擬驗證，計算有效波高平均絕對誤差為0.52 m，最大有效波高相對誤差僅5.3%，模型計算結(jié)果較好地反映了臺風期間的波浪過程，可用于珠江河口波浪—風暴潮耦合數(shù)值模擬研究。

2)基于本文建立的波浪—風暴潮耦合數(shù)值模型，以1822號臺風“山竹”為例，討論了波浪對風暴潮增水的影響，結(jié)果表明波浪對風暴潮增水的影響僅限于近岸海域，最大值為10～20 cm，占風暴潮增水的6.5%?？紤]波浪—風暴潮耦合作用后，計算結(jié)果更接近觀測值，可提高風暴潮數(shù)值預報的精度，同時對沿岸設計潮位的確定具有一定的參考意義。