珠江口沿岸極值增水的空間分布

王燕妮，呂曉鳳，郝嘉凌，吳亞男，邱婷，陳鈺祥，左軍成

（1.河海大學(xué)海洋學(xué)院，江蘇南京 210098；2.北京航天宏圖信息技術(shù)有限公司，北京 100195；3.河海大學(xué)港口海岸與近海工程學(xué)院，江蘇南京 210098；4.國家海洋局東海預(yù)報中心，上海 201200；5.中國水產(chǎn)科學(xué)研究院南海水產(chǎn)研究所，廣東廣州 510300）

珠江口沿岸極值增水的空間分布

王燕妮1，呂曉鳳2，郝嘉凌3，吳亞男1，邱婷4，陳鈺祥5，左軍成1

（1.河海大學(xué)海洋學(xué)院，江蘇南京 210098；2.北京航天宏圖信息技術(shù)有限公司，北京 100195；3.河海大學(xué)港口海岸與近海工程學(xué)院，江蘇南京 210098；4.國家海洋局東海預(yù)報中心，上海 201200；5.中國水產(chǎn)科學(xué)研究院南海水產(chǎn)研究所，廣東廣州 510300）

利用FVCOM模式建立了珠江口區(qū)域的天文潮和風(fēng)暴潮數(shù)值模型，采用調(diào)和分析方法得到該區(qū)域K1、O1、M2、S24個主要分潮的等振幅線和等遲角線的分布特征，并選取典型臺風(fēng)對影響珠江口附近海域的臺風(fēng)暴潮過程進(jìn)行了模擬計算，驗(yàn)證了模型的準(zhǔn)確性，并基于該模型模擬了1993—2013年影響珠江口附近海域的共48場臺風(fēng)過程；利用Gumbel分布，計算了百年一遇的風(fēng)暴潮極值增水，珠江口門外西側(cè)的極值增水大于東側(cè)，其中該區(qū)域百年一遇風(fēng)暴增水的最大值出現(xiàn)在赤溪鎮(zhèn)和三灶鎮(zhèn)中間的海域?yàn)?.19 m，可為珠江口沿岸地區(qū)的風(fēng)暴潮防護(hù)和海岸工程設(shè)計提供參考。

珠江口；百年一遇；風(fēng)暴增水；耿貝爾分布

1 引言

海洋環(huán)境動力要素的重現(xiàn)期值是海岸工程設(shè)計中需要考慮的重要參數(shù)，而由于通常可供利用的觀測資料時限較短，且受空間范圍的局限性，無法從實(shí)測數(shù)據(jù)出發(fā)利用極值分布理論和適應(yīng)法來解決[1-3]。由于多年一遇極值水位要求有19 a以上的觀測資料，而長期驗(yàn)潮站的分布多局限于沿岸，因此多年一遇極值水位的空間分布，必須借助于數(shù)值模擬[4]。珠江口地區(qū)臨近西北太平洋，受臺風(fēng)襲擊頻率高且地勢低平，三江匯集，八口入海這些地理特征都使得該區(qū)域成為我國海岸帶中風(fēng)暴潮災(zāi)害最集中的區(qū)域之一；與此同時，珠江口沿海港口資源豐富、人口密集、經(jīng)濟(jì)發(fā)達(dá)，臺風(fēng)暴潮的襲擊會對當(dāng)?shù)卦斐芍卮髶p失，開展該區(qū)域風(fēng)暴潮增水特征參數(shù)的研究對于當(dāng)?shù)匕l(fā)展意義重大[5]。

臺風(fēng)是影響臺風(fēng)暴潮強(qiáng)度及災(zāi)害等級的主要因素[6]。常用的臺風(fēng)氣壓場模型有高橋模型[7]、Bjerknes模型[8]、Jelesninanski模型[9]、Holland模型[10]以及一些半理論半經(jīng)驗(yàn)的橢圓模型[11]；臺風(fēng)風(fēng)場大多根據(jù)氣壓場計算臺風(fēng)作用范圍內(nèi)的中心對稱風(fēng)場，然后將其與臺風(fēng)移行風(fēng)場進(jìn)行矢量疊加。中心對稱風(fēng)場的計算可以利用藤田公式[12]或者Jelesnianski給出的公式，或者根據(jù)梯度風(fēng)原理得到；移行風(fēng)場的計算可以根據(jù)宮崎正衛(wèi)公式、上野武夫公式和Jelesnianski公式，通過無因次分析方法比較了這幾種臺風(fēng)模型[13]，本文將參照該方法選用高橋和藤田模型來嵌套模擬同一臺風(fēng)域中的氣壓場分布并采用梯度風(fēng)原理計算中心對稱風(fēng)場并采用宮崎正衛(wèi)公式計算移行臺風(fēng)場。

2 數(shù)值模式介紹

FVCOM（The Unstructured Grid Finite VolumeCommunityOceanMode）是一種近岸海洋數(shù)值模式[14]，F(xiàn)VCOM采用無結(jié)構(gòu)三角形網(wǎng)格，適用于復(fù)雜岸線區(qū)域；采用有限體積法求解方程符合能量、動量以及質(zhì)量的守恒條件；引入潮灘動邊界處理，使其在近岸淺灘模擬效果更好，并且可以借助大型服務(wù)器并行計算，極大的提升了計算效率，縮短了調(diào)試與計算時間，目前已經(jīng)被越來越多的學(xué)者用于近岸海洋數(shù)值模擬研究。

σ坐標(biāo)下的連續(xù)方程和水平動量方程如下：

式中：u、v、w是x、y、z方向的速度分量；f為科氏參量；P為壓強(qiáng)；ρ為密度；g為重力加速度；Fx、Fy代表水平動量擴(kuò)散項；Km、Am分別是垂向和水平方向的混合參數(shù)；海表面和海底的邊界條件如下：

在海表面處：

在海底處：

3 模型設(shè)置和驗(yàn)證

3.1 模型設(shè)置

珠江口附近海域的潮汐類型大部分為混合潮型，并且半日潮占優(yōu)，平均潮差在0.85～1.62 m之間，屬于弱潮型河口。本文模擬的區(qū)域范圍（105°～127°E，12°～30°N）南到南海中部，北至杭州灣，東向伸入西北太平洋，基本包括影響珠江口區(qū)域的臺風(fēng)所作用的海域。為了更準(zhǔn)確的貼合研究海域的地形和岸線，根據(jù)計算精度要求對研究區(qū)域進(jìn)行了逐層加密處理[15]，模型網(wǎng)格單元有68 782個，網(wǎng)格節(jié)點(diǎn)有35 041個，開邊界分辨率最小為70 km，網(wǎng)格加密的珠江口近岸海域最高分辨率約200 m，呈現(xiàn)出近岸密集，外海稀疏的網(wǎng)格分布特征（見圖1）。模型的垂向分7層并采用σ坐標(biāo)，時間步長外模為0.1 s，內(nèi)模為1 s，底部粗糙度最小值為0.002 5 m，溫鹽初始場采用常值，模型給定初始水位和初始零流速。模型開邊界給由TMD（Tide Model Drvier）提供M2、S2、K1、O1、N2、P1、Q1、K2這8個主分潮強(qiáng)迫，考慮天文潮和風(fēng)暴潮之間的耦合作用，先對無風(fēng)強(qiáng)迫作用下的天文潮模型進(jìn)行了模擬并利用鰂魚涌站（114.2133°E，22.2911°N）和深圳站（113.88°E，22.47°N）的實(shí)測潮位資料對天文大中小潮進(jìn)行了模擬驗(yàn)證，兩個站點(diǎn)的振幅和相位的擬合效果都比較好，說明模型的開邊界強(qiáng)迫的設(shè)置和網(wǎng)格設(shè)置都比較合理（見圖2）。

對4個主要分潮（K1、O1、M2、S2）進(jìn)行調(diào)和分析，計算出每個分潮的調(diào)和常數(shù)，并繪出同潮圖（見圖3）?？梢钥闯鲞@些調(diào)和常數(shù)的等振幅線及等遲角線在研究海域中表現(xiàn)出的空間分布特征均為近岸密、外海稀，并且兩個半日分潮分布特征相似，兩個全日分潮分布特征相似。對分布趨勢進(jìn)行分析，4個分潮均來源于太平洋潮波，半日分潮M2、S2通過呂宋海峽傳入；全日分潮K1、O1通過呂宋海峽和臺灣海峽傳入。這些分潮都是逆時針方向的旋轉(zhuǎn)潮波系統(tǒng)；在研究海域中，從外海到廣東沿岸，水深逐漸變淺，振幅逐漸增大，在廣東省沿岸，振幅最大的分潮是 K1分潮，O1、M2分潮次之，S2分潮最小。此外，在臺灣海峽中部臺灣淺灘附近，兩個半日分潮的南北兩支在此匯合，振幅隨著水深的迅速變淺而急劇增大，相對的，全日分潮在傳播時耗散的能量較少，因此在臺灣淺灘附近，半日分潮的等振幅線和等遲角線比全日分潮的密集。

圖1 模式網(wǎng)格圖

圖2 天文潮計算值與實(shí)測值比較圖

3.2 風(fēng)暴潮模型驗(yàn)證

本文統(tǒng)計了1993—2013年影響珠江口及其附近海域的共48場臺風(fēng)（見表1、圖4），從這些臺風(fēng)過程的路徑分布可以得出該區(qū)域受風(fēng)暴潮災(zāi)害的影響較大。同時考慮天文潮和風(fēng)暴潮之間的耦合作用，利用FVCOM模型分別模擬了這48場臺風(fēng)期間有臺風(fēng)場作用和無風(fēng)強(qiáng)迫兩種情況下的潮過程，并將兩種情況下得到的逐時風(fēng)暴潮位和天文潮位結(jié)果作差得到臺風(fēng)過程增水值。

圖3 南海北部4個主要分潮的同潮圖（單位：振幅/m，遲角/°）

表1 1993—2013年影響珠江口及其附近海域的臺風(fēng)統(tǒng)計表

圖4 1993—2013年影響珠江口附近海域的歷年臺風(fēng)路徑

為了使模型模擬的增水過程盡可能的接近真實(shí)的臺風(fēng)暴潮過程，本文選取了整個南海北部海域作為模型的研究區(qū)域。同時為了保證模型的準(zhǔn)確性，本文選取了歷史上引起較大增水并對珠江口區(qū)域造成較大影響的典型臺風(fēng)風(fēng)暴潮進(jìn)行了增水過程驗(yàn)證。選取的臺風(fēng)分別是2012年7月20日在菲律賓以東洋面生成并于當(dāng)月24日4時（北京時）登陸廣東赤溪鎮(zhèn)的1208號臺風(fēng)“VICENTE”和2013年9月16日同樣在菲律賓以東洋面生成并以中心氣壓935 hPa的強(qiáng)度登陸廣東汕尾的1319號臺風(fēng)“USAGI”。利用4個珠江口附近區(qū)域的4個驗(yàn)潮站鹽田港站、赤灣站、珠海站和大萬山站實(shí)測資料進(jìn)行了對比驗(yàn)證，并計算分析了風(fēng)暴潮過程中増水的模擬誤差，在實(shí)測增水最大值大于100 cm時采用相對誤差，實(shí)測增水最大值小于等于100 cm時采用絕對誤差，具體誤差分析見表2。

2012年的VICENTE臺風(fēng)過程期間4個站點(diǎn)都出現(xiàn)較大增水，鹽田港站增水的最大值出現(xiàn)在臺風(fēng)將要登陸時，實(shí)測值為1.32 m，模式值為1.08 m，相對誤差為0.18，峰值出現(xiàn)時刻的時間差為3 h；赤灣站最大增水實(shí)測值為1.53 m，模式值為1.33 m，相對誤差為0.13，峰值出現(xiàn)時刻的時間差為1 h；珠海站最大增水實(shí)測值為1.67 m，模式值為1.66 m，相對誤差為0.01，峰值出現(xiàn)時刻的時間差為1 h；大萬山最大增水實(shí)測值為1.0 m，模式值為0.88 m，絕對誤差為0.12 m，峰值出現(xiàn)時刻的時間差為1 h（見圖5）。

2013年的USAGI臺風(fēng)過程期間，鹽田港站最大增水實(shí)測值為0.84 m，模式值為0.75 m，絕對誤差為0.09 m；赤灣站最大實(shí)測增水值為0.51 m，模式值為0.43 m，絕對誤差為0.08 m，峰值出現(xiàn)時刻的時間差為1 h；珠海站最大增水實(shí)測值為0.36 m，模擬值為0.40 m，絕對誤差為0.04 m，峰值出現(xiàn)時刻的時間差為1 h；大萬山站最大增水實(shí)測值為0.37 m，模式值為0.34 m，絕對誤差為0.03 m，峰值出現(xiàn)時刻的時間差為6 h（見圖6）。由于地形資料的誤差和臺風(fēng)模型的不完美，這2個典型臺風(fēng)過程中4個站點(diǎn)模式計算的增水值與實(shí)測數(shù)據(jù)之間存在一定的誤差，但這些誤差都在模型的允許范圍內(nèi)，因此，本文建立的風(fēng)暴潮模型模擬得到的增水結(jié)果是可信的。

4 增水極值統(tǒng)計分析

Gumbel分布方法比較適用于水文氣象等要素多年一遇重現(xiàn)期值的計算，而風(fēng)暴潮增水的經(jīng)驗(yàn)累積頻率點(diǎn)分布近似為一條直線，所以風(fēng)暴潮增水的極值分布符合Gumbel分布[16]。

Gumbel分布函數(shù)為：

其密度函數(shù)為：

式中：α和β為參數(shù)，其值可以根據(jù)下面的方法計算得到：將增水極值x的一組樣本數(shù)據(jù)x1，x2，…，xN從大到小進(jìn)行排序，計算x的數(shù)學(xué)期望和標(biāo)準(zhǔn)差σx，然后可以按照式（10）求出參數(shù)α和β。

將計算得到的α和β代入到式（8）中，經(jīng)整理可得不同重現(xiàn)期極值增水：

表2 臺風(fēng)增水過程對比驗(yàn)證

圖5 1208號臺風(fēng)“VICENTE”期間的風(fēng)暴增水驗(yàn)證

圖6 1309號臺風(fēng)“USAGI”期間的風(fēng)暴增水驗(yàn)證

式中：p′是設(shè)計頻率，也是重現(xiàn)期T的倒數(shù)，φ為Gumbel離差系數(shù)，是 p′的函數(shù)。Cv是偏差系數(shù)[17-18]。

利用本文建立的風(fēng)暴潮模型計算模擬了1993—2013年21 a間影響珠江口附近海域的所有臺風(fēng)增水過程共48場，由于篇幅限制下文只給出了珠江口區(qū)域在0806號臺風(fēng)、0814號、1208號和1319號臺風(fēng)期間的增水極值分布（見圖8）。我們統(tǒng)計得到21 a的風(fēng)暴增水年極值，如此每一個空間點(diǎn)上對應(yīng)21個增水極值，對每一個空間點(diǎn)利用Gumbel分布模式統(tǒng)計計算得到每一個空間點(diǎn)的重現(xiàn)期為100 a的增水極值，從而就得到了該區(qū)域重現(xiàn)期為100 a的增水極值的空間分布（見圖9）。

圖7 4次典型臺風(fēng)風(fēng)暴增水極值分布

在圖7的4個臺風(fēng)過程中，珠江口區(qū)域風(fēng)暴增水極值的最大值分別為1.12 m、2.52 m、2.49 m、2.03 m。這幾場臺風(fēng)增水極值最大值出現(xiàn)的位置分別在深圳灣、赤溪鎮(zhèn)東側(cè)海域、赤溪鎮(zhèn)東側(cè)海域和三灶鎮(zhèn)西側(cè)海域、惠東縣沿海。臺風(fēng)增水最大值出現(xiàn)的位置主要受臺風(fēng)路徑和當(dāng)?shù)氐匦蔚墓餐绊?。該區(qū)域百年一遇風(fēng)暴增水的最大值出現(xiàn)在赤溪鎮(zhèn)和三灶鎮(zhèn)中間的海域?yàn)?.19 m，珠江口內(nèi)的最大值出現(xiàn)在南沙港附近為3.08 m，口門處的重現(xiàn)期值在1.52～2.02 m的范圍內(nèi)，珠江口東側(cè)的重現(xiàn)期值在1.71～2.68 m的范圍內(nèi)（見圖8），與前人利用實(shí)測潮位資料進(jìn)行統(tǒng)計分析得到的多年一遇的極值水位的結(jié)果相近[19]，珠江口內(nèi)及口門西側(cè)百年一遇的增水極值較大，沿岸受風(fēng)暴潮災(zāi)害的影響明顯大于東側(cè)沿岸區(qū)域，應(yīng)重點(diǎn)加強(qiáng)該區(qū)域的風(fēng)暴潮防護(hù)工作。

圖8 重現(xiàn)期為100 a的風(fēng)暴增水極值分布

5 結(jié)論

本研究基于FVCOM模式建立珠江口附近海域的風(fēng)暴潮模型，對珠江口沿岸極值增水的空間分布經(jīng)過計算分析，得出如下結(jié)論：

（1）通過鰂魚涌站和深圳站的實(shí)測潮位資料對天文大中小潮進(jìn)行了模擬驗(yàn)證，結(jié)果良好，表明模型的網(wǎng)格及開邊界設(shè)置合理，可以進(jìn)一步用來模擬風(fēng)暴潮；

（2）研究區(qū)域各分潮均成逆時針方向旋轉(zhuǎn)，廣東沿岸分潮振幅大于外海，其中全日分潮K1的振幅最大，O1、M2分潮次之，S2分潮最小，風(fēng)暴增水極易出現(xiàn)在近海沿岸分潮振幅較大區(qū)域；

（3）對引起較大增水并對珠江口區(qū)域造成較大影響的兩個典型臺風(fēng)“VICENTE”和“USAGI”進(jìn)行了風(fēng)暴增水過程驗(yàn)證，結(jié)果良好，表明風(fēng)暴潮模型具有較高的準(zhǔn)確度；

（4）基于風(fēng)暴潮模型計算出1993—2013年期間影響珠江口及其附近海域的共48場臺風(fēng)的臺風(fēng)過程增水值，并利用Gumbel分布得到研究海域內(nèi)重現(xiàn)期為100 a的增水極值的空間分布；

（5）重新期為100 a的風(fēng)暴增水極值在珠江口西側(cè)和珠江口內(nèi)的值較大，最大值出現(xiàn)在赤溪鎮(zhèn)和三灶鎮(zhèn)中間的海域?yàn)?.19 m，南沙港附近的增水極值為3.08 m。在制定風(fēng)暴潮防護(hù)規(guī)劃和工程設(shè)計時應(yīng)對這些區(qū)域予以重視。

[1]侍茂崇,高郭平,鮑獻(xiàn)文.海洋調(diào)查方法[M].青島:青島海洋大學(xué)出版社,2000:312-318.

[2]張月霞,王輝.臺風(fēng)風(fēng)暴潮災(zāi)害風(fēng)險評估研究綜述[J].海洋預(yù)報,2016,33(2):81-88.

[3]王喜年.風(fēng)暴潮預(yù)報知識講座——第七講風(fēng)暴潮風(fēng)險分析與計算[J].海洋預(yù)報,2002,19(4):73-76.

[4]于宜法,俞聿修.渤海天文-風(fēng)暴潮數(shù)值模擬和一種多年一遇極值水位的計算方法[J].海洋學(xué)報,2003,25(4):10-17.

[5]高義,王輝,劉桂梅.基于CVISS的風(fēng)暴潮災(zāi)害脆弱性評估研究[J].海洋預(yù)報,2013,30(6):1-6.

[6]梁必騏.熱帶氣象學(xué)[M].廣州:中山大學(xué)出版社,1990:213-219.

[7]馮士筰.風(fēng)暴潮導(dǎo)論[M].北京:科學(xué)出版社,1982.

[8]Bjerknes V.On the dynamics of the circular vortex:with applications to the atmosphere and atmospheric vortex and wave motions[M].Kristiana (Oslo),Norway:Gr?ndahl& S?ns Boktrykkeri,1921:704-710.

[9]Jelesnianski C P.A numerical calculation of storm tides induced by a tropical storm impinging on a continental shelf[J].Monthly Weather Review,1965,93(6):343-358.

[10]Holland G J.An analytic model of the wind and pressure profiles in Hurricanes[J].Monthly Weather Review,1980,108(8):1212-1218.

[11]陳紹偉.使用多普勒速度進(jìn)行臺風(fēng)擬合定位的試驗(yàn)[J].氣象研究與應(yīng)用,2007,28(S2):95-97.

[12]Fujita T.Pressure distribution within typhoon[J].Geophysical Magazine,1952,23(4):437-451.

[13]王喜年.第三講風(fēng)暴潮數(shù)值模式計算中氣壓場和風(fēng)場的處理[J].海洋預(yù)報,1986,3(4):56-64.

[14]Chen C S,Beardsley R C,Cowles G.An unstructured grid,finitevolume coastal ocean model(FVCOM)system[J].Oceanography,2006,19(1):78-89.

[15]朱軍政,曹穎.FVCOM模型在象山港三維潮流鹽度計算中的應(yīng)用[J].海洋環(huán)境科學(xué),2010,29(6):899-903.

[16]吳亞楠,王智峰,董勝,等.山東沿海臺風(fēng)暴潮數(shù)值模擬與統(tǒng)計分析[J].自然災(zāi)害學(xué)報,2015,24(3):169-176.

[17]李明杰,齊鵬,侯一筠.山東沿岸多年一遇最高水位計算[J].海洋科學(xué),2009,33(11):78-81.

[18]候瑞科.利用耿貝爾極值分布計算年最高水位[J].海洋通報,1993,12(3):126-129.

[19]李闊,李國勝.珠江三角洲地區(qū)風(fēng)暴潮重現(xiàn)期及增水與環(huán)境要素的關(guān)系[J].地理科學(xué)進(jìn)展,2010,29(4):433-438.

Spatial distribution of extreme water enhancement along the Pearl River Estuary

WANG Yan-ni1，LV Xiao-feng2，HAO Jia-ling3，WU Ya-nan1，QIU Ting4，CHEN Yu-xiang5，ZUO Jun-cheng1

（1.College of Oceanography,Hohai University,Nanjing 210098 China;2.Beijing Piesat Information Technology Co,.Ltd,Beijing 100195 China;3.Collage of Harbor Coastal and Offshore Engineering,Hohai University,Nanjing 210098 China;4.Forecast Center for East China Sea,SOA,Shanghai 201200 China;5.South China Sea Fisheries Research Institute,Chinese Academy of Fishery Sciences,Guangzhou 510300 China）

The numerical model of the astronomical tide and storm surge in the Pearl River estuary region is established by using FVCOM model.The distribution characteristics of K1,O1,S2and M2are obtained by using harmonic analysis method.The typical typhoon storm surge processed near the Pearl River Estuary is simulated to verify the accuracy of the model.Based on the model,the influence of 100 typhoons during the years 1993—2013 in the vicinity of the Pearl River estuary is simulated.By using the Gumbel distribution,the storm surge extreme water increment is calculated.The storm surge in the west side of the Pearl River estuary is more extreme than that in the east.The once-in-a-century storm surge maximum value 3.19m appears in the middle of the town San Zao and the town Chi Xi.This paper can provide reference for storm surge protection and coastal engineering design in the Pearl River Estuary.

Pearl River Estuary；one hundred years；storm surge；Gumbel distribution

P731.23

A

1003-0239（2017）05-0074-09

10.11737/j.issn.1003-0239.2017.05.008

2017-04-14；

2017-06-15。

國家自然科學(xué)基金（41576020）；國家社會科學(xué)基金重大項目（15ZDB170）；國際合作項目課題（GASI-IPOVAI-06）。

王燕妮（1991-），女，助理工程師，碩士，主要從事海岸海洋環(huán)境與動力災(zāi)害研究。E-mail：410366004@qq.com

呂曉鳳（1990-），女，助理工程師，碩士，主要從事海岸災(zāi)害防護(hù)研究。E-mail：764655240@qq.com