基于WRF-波浪风暴潮耦合模型的深圳市极端风暴潮安全分析

林泽辉

(惠州市惠阳区永良堤围管理所，广东 惠州 516007)

我国东南沿海地区气象水文条件复杂，水文方面，海水水位受天体引潮力作用而周期性涨落。气象方面，西太平洋热带地区气压受到扰动作用成涡，并不断发展移动，这些气旋登陆后会对人类社会造成巨大的破坏。另一方面，台风所造成的大气压强场变化也会导致海面异常起伏，产生风暴潮，当风暴潮恰好叠加上天文大潮，海平面会巨幅上升，这对海岸防护、河口挡潮等水利工程带来了巨大挑战。因此研究台风过程前后的波浪以及增水具有重要意义。文章利用…模型对深圳地区水域在台风作用下的水动力特征进行了数值模拟研究。

1 研究方法

1.1 区域概况

深圳市海岸线长约250km，东起大亚湾，西至珠江口，增水>50cm的风暴潮过程年平均出现2.1次。沿岸分布有机场、重点港口、大型住宅区、核电站等大型企业，人口密度极高。

深圳市年平均风速2.6m/s，冬季各月风速较大为3.0m/s左右，夏季各月较小约2.0m/s。大风日数(风速≥17m/s，即风力相当于8级或8级以上)年平均为7.3d。全年各月均出现大风，其中以7-9月居多，共占全年大风日数的44%，其次为6月和10月，由于受台风或热带低压影响，常出现大风，因此，大风出现多的月份，也是台风季节。

1.2 数据资料

影响风暴增水量级的因素有很多，包括台风气压场分布、台风风场的方向与大小以及近岸地形变化因素。在从我国东南沿海登陆的台风里，影响到阳江至汕尾海域的台风均可引起深圳海域的风暴潮增水。目前在深圳海域，国家海洋局共设有两个长期验潮站，分别是赤湾和盐田站，紧邻深圳大亚湾岸线的有位于惠州市的惠州海洋站，各站位置见图1。收集上述站点自建站以来台风期间的潮位资料统计历史风暴潮增水情况。广东省水文局1964年在深圳市赤湾港内建设有长期水文观测站，1986年国家海洋局南海分局在相距不到1km处建设了赤湾海洋站，并于1986年开始投入业务化观测。盐田海洋站位于大鹏湾盐田港海域，由国家海洋局南海分局建设和管理，2003年投入使用。惠州海洋站位于大亚湾湾顶西侧荃湾港区，由国家海洋局南海分局建设和管理，2006年建成投入使用，他们的分布位置如图1所示。

图1 资料站位位置示意图

1.3 研究方法

早期人们主要通过经验总结归纳的方法研究风暴潮，观察台风天时海塘外潮位变化研究风暴潮的成因。20世纪中期，随着计算机技术的发展，人们开始利用数值模拟研究风暴潮。1956年W·Hansen首次利用计算机对欧洲北海的最大风暴潮增水进行了研究。随着数值模拟技术的进步，Jelesnianski建立了研究风暴潮的SPLASH模式以及基于二维水动力的SLOSH模式，在不断地摸索研究后，该模式的预报精度不断提高。近年来，不断发展的模型计算实现了对复杂地形影响的考虑，提高了计算的准确性[1]。其次，有关风暴增水与天文潮的非线性拟合机制也变得更加清晰，这也提高了风暴潮数学模型的准确性。我国自20世纪60年代也开展了关于风暴潮的相关研究。冯士筰所著的《风暴潮导论》系统阐述了风暴潮的概念，形成机制，并介绍了早期的风暴潮数值预报方法。赵长进等利用ADCIRC潮流与SWAN波浪耦合模式对长江口附近海域的风暴潮过程进行了模拟。黄本胜等用SELFE模型构建了南海—珠江河口双重嵌套模型，研究了“山竹”台风过程中珠江河口的风暴增水与波浪增水分布。

在波浪模拟方面，目前工程界使用较为广泛的数学模型有两类，一类是波浪流体动力学模型，另一类是动谱密度方程。缓坡方程模型与Boussinesq方程模型是主要的波浪流体力学模型。经典的Boussinesq仅仅适用于浅水，并不适用于水深变化大的水域和不规则波的情况。因此在工程中适用范围不广。动谱密度方程的基础是能量平衡理论，在模拟过程中可提高海岸波浪场的动力作用过程的计算效率。该方程以动谱密度来描述波浪变形，考虑了对波-波非线性相互作用、波流相互作用、波浪破碎和底部摩擦引起的能量耗散等过程，通过将源项添加到方程中，可以全面考虑动力作用的复杂过程，使其能完成大尺度区域的波浪产生、成长、传播和变形的计算过程。其中对物理源项的处理和研究是波作用量守恒的关键，波作用量守恒方程发展到现在已经到了第三代，逐渐开发出WAVEWATCH、SWAN、MIKE21SW、TOMAWAC等模型。

1.4 模型设置

1.4.1 模型介绍

文章采用的是DELFT3D-FLOW与SWAN的耦合模型，DELFT3D是一款由荷兰DELTRA开发的开源水动力模拟软件，在风暴潮模拟方面已经广泛应用。

Delft3D-Flow是基于浅水方程和Boussinesq假设，采用有限差分法求解斜压Navier-Stokes方程的模型，模型控制方程组分别为：

质量守恒方程：

(1)

式中：ζ为水位(总水深)；d为到参考平面的当前水深(净水深)；U、V-x、y为方向上的速度分量；Q为单位面积上质量强度。

动量守恒方程：

(2)

(3)

式中：f为科氏力参数；g为重力加速度；νh为动态水平涡流黏度；τsx,τsy为x、y方向上的海平面风压分量；τbx,τby为x、y方向上的底部的剪应力分量；ρ0为参考密度；ρ'为不规则密度。

SWAN模型的控制方程为：

(4)

式中：Cx，Cy，Cσ，Cθ为二维动谱密度在x、y、σ、θ为空间上的传播速度；σ为相对频率；θ为波向；S为能量源项。

1.4.2 模型设置

1)模型范围：

文章所采用的模型囊括整个珠江三角洲区域，见图2。

图2 模型区域示意图

2)地形水深：

模拟范围内的地形水深采用的是中华人民共和国海事局2016年绘制的海图插值数据。将地形水深数据转化到统一基面(平均海平面)后，用IDW插值方法插值到网格点上，得到模型的水深地形。

3)开边界设置：

开边界采用天文潮水位，主要考虑M2、N2、S2、K2、K1、O1、P1、Q1等八个天文分潮以及MS4、M4、M6等浅水分潮共11个分潮的作用。分潮的调和分析常数由2016-2019年的潮位实测数据经MATLAB中的T_TIDE工具包调和分析所得。模型上部采用流量边界，输入数据为梧州、石角、博罗三个水文站的2016-2019洪季平均流量。

4)底摩擦设置：

模型根据水深范围设置不同大小的糙率，水深越大糙率越小。

1.4.3 风场输入

文章的风场采用美国国家大气研究中心(NCAR)等机构共同开发的中尺度大气预报系统。WRF模式适用于台风模拟，可以通过嵌套网格运算来达到高精度模拟的目的。本次采用了WSM6微物理方案、KainFritsch卷积云方案、Monin-Obukhov表面层方案、Noah陆面层方案、MYJ大气边界层方案、RRTMG辐射方案模拟台风。

2 模型验证

2.1 风场验证

首先利用实测资料对所建立的风暴潮-波浪耦合模型进行率定验证。选取1713号台风(2017年第13号台风—天鸽，其移行路径如图3所示)作为模型边界条件，为了证明WRF模式的准确性其在风暴潮模拟方面的优势，文章也用其他论文常用的经验风场模型做了验证对比。经验风场采用的是宫崎正卫模型作为移动风场，HOLLAND模型作为梯度风场，他们的公式分别为：

图3 “天鸽”台风移行路径

(5)

(5)

式中：B为HOLLAND模型气压参数；f为科氏力参数，与模拟区域纬度有关；R为最大风速半径，与中心气压有关。

图4为WRF模型模拟风场值以及经验风场模拟值与实测风场的对比，可以看出，WRF模型模拟的风场在峰值以及过程上都比经验风场模型更贴近实际，因此利用WRF模型生成的台风风场作为风暴潮-波浪耦合模型的边界条件更为合理。

图4 WRF风场及模型风场与实测风场对比

2.2 风暴潮及波浪率定验证

以1713号台风(2017年第13号台风)为边界条件分别驱动珠江口风暴潮模型以及珠江口波浪-风暴潮耦合模型，得到如图5所示的赤湾站与盐田站水位验证曲线，可以看出，考虑波浪后的水位过程曲线更符合实际，最大误差绝对值为13cm，平均相对误差为11.3%，标准差为2.94，符合模拟精度要求。

图5(a) 赤湾站不同模型下水位曲线

图5(b) 盐田站不同模型下水位曲线

3 深圳市沿岸漫堤风险评估

3.1 可能最大风暴潮模型计算

为了寻求可能影响深圳市的最大风暴潮对深圳市沿岸的影响，采用概率论法进行计算可能最大台风中的近中心最大风速Vmax，台风中心气压P0。

利用中国气象局(CMA)1949-2015年热带气旋最佳路径集(Best-track)数据，选取出研究区域400km直径范围内对本研究区域有重要影响的台风，以所选取台风的最小台风中心气压值作样本。采用极值I型分布计算1000年一遇的台风中心气压值作为可能最大台风的中心气压，计算结果如图6(a)所示，即深圳区域1000a遇的台风中心气压P0=880hPa。

图6(a) 影响深圳的台风中心最低气压极值I型分布

图6(b) 影响深圳的台风中心最大风速极值I型分布

利用中国气象局(CMA)1949-2015年热带气旋最佳路径集(Best-track)数据，取深圳市海湾区域400km直径范围内历年路经本海区的台风最大风速值作样本，对于没有台风经过研究区域的年份，该年Vmax取为台风Vmax系列中的平均值。采用极值I型分布计算1000a一遇的最大风速值作为可能最大台风的最大风速值，计算结果如图6(b)所示，深圳区域1000a一遇的台风中心最大风速Vmax=85m/s。另外，台风路径选取为东南向西北由大鹏湾登陆的路径。

输出堤前代表点(具体见图7)的最高潮位与堤顶高程进行比对，了解沿岸淹没的情况，由于模式计算结果是相对平均海平面的，为便于比对，需要根据当地平均海平面与85高程的关系，将计算结果转成以85黄海基面为起算面，其中西段提防采用赤湾站高程关系，即计算结果加上0.54m；东段提防采用盐田站高程关系进行转换，即计算结果加上0.66m。表1列出堤前统计的最高潮位及堤顶高程。从代表点的最高潮位和海堤高程比较来看，在可能最大热带气旋影响下，深圳沿岸均可能发生漫堤。

图7 堤前代表点位置示意图

表1 可能最大风暴潮漫堤分析 (单位：m，基面：85高程)

续表1 可能最大风暴潮漫堤分析 (单位：m，基面：85高程)

3.2 海堤防御标准

对于大部分沿海地区来说，1000a一遇的标准是过高的。但是深圳作为特大型城市，其风暴潮安全具有重要意义，因此考虑1000a一遇情形下的风暴潮安全是有必要的，在此情形下得出的淹没情况有利于应急管理部门采取相应的防范措施。

4 结 论

1)利用WRF中长尺度气象预报模式作为波浪-风暴潮耦合模型的风场边界条件输入，提高了风暴潮模型的精确度

2)在珠江口波浪-风暴潮模型验证实测资料准确的基础上，利用极值I型分布确定1000a一遇台风参数对影响深圳市海岸堤防未来发生的最大可能风暴潮进行计算分析。文章研究成果可为深圳市风暴潮安全防范提供重要技术支撑，同时也为海洋水利防灾应急部门提供决策参考。