王燕妮，吕晓凤，郝嘉凌，吴亚男，邱婷，陈钰祥，左军成

（1.河海大学海洋学院，江苏南京 210098；2.北京航天宏图信息技术有限公司，北京 100195；3.河海大学港口海岸与近海工程学院，江苏南京 210098；4.国家海洋局东海预报中心，上海 201200；5.中国水产科学研究院南海水产研究所，广东广州 510300）

珠江口沿岸极值增水的空间分布

王燕妮1，吕晓凤2，郝嘉凌3，吴亚男1，邱婷4，陈钰祥5，左军成1

（1.河海大学海洋学院，江苏南京 210098；2.北京航天宏图信息技术有限公司，北京 100195；3.河海大学港口海岸与近海工程学院，江苏南京 210098；4.国家海洋局东海预报中心，上海 201200；5.中国水产科学研究院南海水产研究所，广东广州 510300）

利用FVCOM模式建立了珠江口区域的天文潮和风暴潮数值模型，采用调和分析方法得到该区域K1、O1、M2、S24个主要分潮的等振幅线和等迟角线的分布特征，并选取典型台风对影响珠江口附近海域的台风暴潮过程进行了模拟计算，验证了模型的准确性，并基于该模型模拟了1993—2013年影响珠江口附近海域的共48场台风过程；利用Gumbel分布，计算了百年一遇的风暴潮极值增水，珠江口门外西侧的极值增水大于东侧，其中该区域百年一遇风暴增水的最大值出现在赤溪镇和三灶镇中间的海域为3.19 m，可为珠江口沿岸地区的风暴潮防护和海岸工程设计提供参考。

珠江口；百年一遇；风暴增水；耿贝尔分布

1 引言

海洋环境动力要素的重现期值是海岸工程设计中需要考虑的重要参数，而由于通常可供利用的观测资料时限较短，且受空间范围的局限性，无法从实测数据出发利用极值分布理论和适应法来解决[1-3]。由于多年一遇极值水位要求有19 a以上的观测资料，而长期验潮站的分布多局限于沿岸，因此多年一遇极值水位的空间分布，必须借助于数值模拟[4]。珠江口地区临近西北太平洋，受台风袭击频率高且地势低平，三江汇集，八口入海这些地理特征都使得该区域成为我国海岸带中风暴潮灾害最集中的区域之一；与此同时，珠江口沿海港口资源丰富、人口密集、经济发达，台风暴潮的袭击会对当地造成重大损失，开展该区域风暴潮增水特征参数的研究对于当地发展意义重大[5]。

台风是影响台风暴潮强度及灾害等级的主要因素[6]。常用的台风气压场模型有高桥模型[7]、Bjerknes模型[8]、Jelesninanski模型[9]、Holland模型[10]以及一些半理论半经验的椭圆模型[11]；台风风场大多根据气压场计算台风作用范围内的中心对称风场，然后将其与台风移行风场进行矢量叠加。中心对称风场的计算可以利用藤田公式[12]或者Jelesnianski给出的公式，或者根据梯度风原理得到；移行风场的计算可以根据宫崎正卫公式、上野武夫公式和Jelesnianski公式，通过无因次分析方法比较了这几种台风模型[13]，本文将参照该方法选用高桥和藤田模型来嵌套模拟同一台风域中的气压场分布并采用梯度风原理计算中心对称风场并采用宫崎正卫公式计算移行台风场。

2 数值模式介绍

FVCOM（The Unstructured Grid Finite VolumeCommunityOceanMode）是一种近岸海洋数值模式[14]，FVCOM采用无结构三角形网格，适用于复杂岸线区域；采用有限体积法求解方程符合能量、动量以及质量的守恒条件；引入潮滩动边界处理，使其在近岸浅滩模拟效果更好，并且可以借助大型服务器并行计算，极大的提升了计算效率，缩短了调试与计算时间，目前已经被越来越多的学者用于近岸海洋数值模拟研究。

σ坐标下的连续方程和水平动量方程如下：

式中：u、v、w是x、y、z方向的速度分量；f为科氏参量；P为压强；ρ为密度；g为重力加速度；Fx、Fy代表水平动量扩散项；Km、Am分别是垂向和水平方向的混合参数；海表面和海底的边界条件如下：

在海表面处：

在海底处：

3 模型设置和验证

3.1 模型设置

珠江口附近海域的潮汐类型大部分为混合潮型，并且半日潮占优，平均潮差在0.85～1.62 m之间，属于弱潮型河口。本文模拟的区域范围（105°～127°E，12°～30°N）南到南海中部，北至杭州湾，东向伸入西北太平洋，基本包括影响珠江口区域的台风所作用的海域。为了更准确的贴合研究海域的地形和岸线，根据计算精度要求对研究区域进行了逐层加密处理[15]，模型网格单元有68 782个，网格节点有35 041个，开边界分辨率最小为70 km，网格加密的珠江口近岸海域最高分辨率约200 m，呈现出近岸密集，外海稀疏的网格分布特征（见图1）。模型的垂向分7层并采用σ坐标，时间步长外模为0.1 s，内模为1 s，底部粗糙度最小值为0.002 5 m，温盐初始场采用常值，模型给定初始水位和初始零流速。模型开边界给由TMD（Tide Model Drvier）提供M2、S2、K1、O1、N2、P1、Q1、K2这8个主分潮强迫，考虑天文潮和风暴潮之间的耦合作用，先对无风强迫作用下的天文潮模型进行了模拟并利用鰂鱼涌站（114.2133°E，22.2911°N）和深圳站（113.88°E，22.47°N）的实测潮位资料对天文大中小潮进行了模拟验证，两个站点的振幅和相位的拟合效果都比较好，说明模型的开边界强迫的设置和网格设置都比较合理（见图2）。

对4个主要分潮（K1、O1、M2、S2）进行调和分析，计算出每个分潮的调和常数，并绘出同潮图（见图3）。可以看出这些调和常数的等振幅线及等迟角线在研究海域中表现出的空间分布特征均为近岸密、外海稀，并且两个半日分潮分布特征相似，两个全日分潮分布特征相似。对分布趋势进行分析，4个分潮均来源于太平洋潮波，半日分潮M2、S2通过吕宋海峡传入；全日分潮K1、O1通过吕宋海峡和台湾海峡传入。这些分潮都是逆时针方向的旋转潮波系统；在研究海域中，从外海到广东沿岸，水深逐渐变浅，振幅逐渐增大，在广东省沿岸，振幅最大的分潮是 K1分潮，O1、M2分潮次之，S2分潮最小。此外，在台湾海峡中部台湾浅滩附近，两个半日分潮的南北两支在此汇合，振幅随着水深的迅速变浅而急剧增大，相对的，全日分潮在传播时耗散的能量较少，因此在台湾浅滩附近，半日分潮的等振幅线和等迟角线比全日分潮的密集。

图1 模式网格图

图2 天文潮计算值与实测值比较图

3.2 风暴潮模型验证

本文统计了1993—2013年影响珠江口及其附近海域的共48场台风（见表1、图4），从这些台风过程的路径分布可以得出该区域受风暴潮灾害的影响较大。同时考虑天文潮和风暴潮之间的耦合作用，利用FVCOM模型分别模拟了这48场台风期间有台风场作用和无风强迫两种情况下的潮过程，并将两种情况下得到的逐时风暴潮位和天文潮位结果作差得到台风过程增水值。

图3 南海北部4个主要分潮的同潮图（单位：振幅/m，迟角/°）

表1 1993—2013年影响珠江口及其附近海域的台风统计表

图4 1993—2013年影响珠江口附近海域的历年台风路径

为了使模型模拟的增水过程尽可能的接近真实的台风暴潮过程，本文选取了整个南海北部海域作为模型的研究区域。同时为了保证模型的准确性，本文选取了历史上引起较大增水并对珠江口区域造成较大影响的典型台风风暴潮进行了增水过程验证。选取的台风分别是2012年7月20日在菲律宾以东洋面生成并于当月24日4时（北京时）登陆广东赤溪镇的1208号台风“VICENTE”和2013年9月16日同样在菲律宾以东洋面生成并以中心气压935 hPa的强度登陆广东汕尾的1319号台风“USAGI”。利用4个珠江口附近区域的4个验潮站盐田港站、赤湾站、珠海站和大万山站实测资料进行了对比验证，并计算分析了风暴潮过程中増水的模拟误差，在实测增水最大值大于100 cm时采用相对误差，实测增水最大值小于等于100 cm时采用绝对误差，具体误差分析见表2。

2012年的VICENTE台风过程期间4个站点都出现较大增水，盐田港站增水的最大值出现在台风将要登陆时，实测值为1.32 m，模式值为1.08 m，相对误差为0.18，峰值出现时刻的时间差为3 h；赤湾站最大增水实测值为1.53 m，模式值为1.33 m，相对误差为0.13，峰值出现时刻的时间差为1 h；珠海站最大增水实测值为1.67 m，模式值为1.66 m，相对误差为0.01，峰值出现时刻的时间差为1 h；大万山最大增水实测值为1.0 m，模式值为0.88 m，绝对误差为0.12 m，峰值出现时刻的时间差为1 h（见图5）。

2013年的USAGI台风过程期间，盐田港站最大增水实测值为0.84 m，模式值为0.75 m，绝对误差为0.09 m；赤湾站最大实测增水值为0.51 m，模式值为0.43 m，绝对误差为0.08 m，峰值出现时刻的时间差为1 h；珠海站最大增水实测值为0.36 m，模拟值为0.40 m，绝对误差为0.04 m，峰值出现时刻的时间差为1 h；大万山站最大增水实测值为0.37 m，模式值为0.34 m，绝对误差为0.03 m，峰值出现时刻的时间差为6 h（见图6）。由于地形资料的误差和台风模型的不完美，这2个典型台风过程中4个站点模式计算的增水值与实测数据之间存在一定的误差，但这些误差都在模型的允许范围内，因此，本文建立的风暴潮模型模拟得到的增水结果是可信的。

4 增水极值统计分析

Gumbel分布方法比较适用于水文气象等要素多年一遇重现期值的计算，而风暴潮增水的经验累积频率点分布近似为一条直线，所以风暴潮增水的极值分布符合Gumbel分布[16]。

Gumbel分布函数为：

其密度函数为：

式中：α和β为参数，其值可以根据下面的方法计算得到：将增水极值x的一组样本数据x1，x2，…，xN从大到小进行排序，计算x的数学期望和标准差σx，然后可以按照式（10）求出参数α和β。

将计算得到的α和β代入到式（8）中，经整理可得不同重现期极值增水：

表2 台风增水过程对比验证

图5 1208号台风“VICENTE”期间的风暴增水验证

图6 1309号台风“USAGI”期间的风暴增水验证

式中：p′是设计频率，也是重现期T的倒数，φ为Gumbel离差系数，是 p′的函数。Cv是偏差系数[17-18]。

利用本文建立的风暴潮模型计算模拟了1993—2013年21 a间影响珠江口附近海域的所有台风增水过程共48场，由于篇幅限制下文只给出了珠江口区域在0806号台风、0814号、1208号和1319号台风期间的增水极值分布（见图8）。我们统计得到21 a的风暴增水年极值，如此每一个空间点上对应21个增水极值，对每一个空间点利用Gumbel分布模式统计计算得到每一个空间点的重现期为100 a的增水极值，从而就得到了该区域重现期为100 a的增水极值的空间分布（见图9）。

图7 4次典型台风风暴增水极值分布

在图7的4个台风过程中，珠江口区域风暴增水极值的最大值分别为1.12 m、2.52 m、2.49 m、2.03 m。这几场台风增水极值最大值出现的位置分别在深圳湾、赤溪镇东侧海域、赤溪镇东侧海域和三灶镇西侧海域、惠东县沿海。台风增水最大值出现的位置主要受台风路径和当地地形的共同影响。该区域百年一遇风暴增水的最大值出现在赤溪镇和三灶镇中间的海域为3.19 m，珠江口内的最大值出现在南沙港附近为3.08 m，口门处的重现期值在1.52～2.02 m的范围内，珠江口东侧的重现期值在1.71～2.68 m的范围内（见图8），与前人利用实测潮位资料进行统计分析得到的多年一遇的极值水位的结果相近[19]，珠江口内及口门西侧百年一遇的增水极值较大，沿岸受风暴潮灾害的影响明显大于东侧沿岸区域，应重点加强该区域的风暴潮防护工作。

图8 重现期为100 a的风暴增水极值分布

5 结论

本研究基于FVCOM模式建立珠江口附近海域的风暴潮模型，对珠江口沿岸极值增水的空间分布经过计算分析，得出如下结论：

（1）通过鰂鱼涌站和深圳站的实测潮位资料对天文大中小潮进行了模拟验证，结果良好，表明模型的网格及开边界设置合理，可以进一步用来模拟风暴潮；

（2）研究区域各分潮均成逆时针方向旋转，广东沿岸分潮振幅大于外海，其中全日分潮K1的振幅最大，O1、M2分潮次之，S2分潮最小，风暴增水极易出现在近海沿岸分潮振幅较大区域；

（3）对引起较大增水并对珠江口区域造成较大影响的两个典型台风“VICENTE”和“USAGI”进行了风暴增水过程验证，结果良好，表明风暴潮模型具有较高的准确度；

（4）基于风暴潮模型计算出1993—2013年期间影响珠江口及其附近海域的共48场台风的台风过程增水值，并利用Gumbel分布得到研究海域内重现期为100 a的增水极值的空间分布；

（5）重新期为100 a的风暴增水极值在珠江口西侧和珠江口内的值较大，最大值出现在赤溪镇和三灶镇中间的海域为3.19 m，南沙港附近的增水极值为3.08 m。在制定风暴潮防护规划和工程设计时应对这些区域予以重视。

[1]侍茂崇,高郭平,鲍献文.海洋调查方法[M].青岛:青岛海洋大学出版社,2000:312-318.

[2]张月霞,王辉.台风风暴潮灾害风险评估研究综述[J].海洋预报,2016,33(2):81-88.

[3]王喜年.风暴潮预报知识讲座——第七讲风暴潮风险分析与计算[J].海洋预报,2002,19(4):73-76.

[4]于宜法,俞聿修.渤海天文-风暴潮数值模拟和一种多年一遇极值水位的计算方法[J].海洋学报,2003,25(4):10-17.

[5]高义,王辉,刘桂梅.基于CVISS的风暴潮灾害脆弱性评估研究[J].海洋预报,2013,30(6):1-6.

[6]梁必骐.热带气象学[M].广州:中山大学出版社,1990:213-219.

[7]冯士筰.风暴潮导论[M].北京:科学出版社,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].气象研究与应用,2007,28(S2):95-97.

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

[13]王喜年.第三讲风暴潮数值模式计算中气压场和风场的处理[J].海洋预报,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模型在象山港三维潮流盐度计算中的应用[J].海洋环境科学,2010,29(6):899-903.

[16]吴亚楠,王智峰,董胜,等.山东沿海台风暴潮数值模拟与统计分析[J].自然灾害学报,2015,24(3):169-176.

[17]李明杰,齐鹏,侯一筠.山东沿岸多年一遇最高水位计算[J].海洋科学,2009,33(11):78-81.

[18]候瑞科.利用耿贝尔极值分布计算年最高水位[J].海洋通报,1993,12(3):126-129.

[19]李阔,李国胜.珠江三角洲地区风暴潮重现期及增水与环境要素的关系[J].地理科学进展,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。

国家自然科学基金（41576020）；国家社会科学基金重大项目（15ZDB170）；国际合作项目课题（GASI-IPOVAI-06）。

王燕妮（1991-），女，助理工程师，硕士，主要从事海岸海洋环境与动力灾害研究。E-mail：410366004@qq.com

吕晓凤（1990-），女，助理工程师，硕士，主要从事海岸灾害防护研究。E-mail：764655240@qq.com