街区尺度风暴潮漫滩数值模拟

张露，傅赐福，董剑希,2，于福江,2

( 1. 国家海洋环境预报中心，北京 100081；2. 国家海洋局海洋灾害预报技术研究重点实验室，北京 100081)

1 引言

我国是世界上少数几个海洋灾害十分严重的国家之一，尤其以风暴潮灾害为主，约占总灾害损失的90%以上，风暴潮不仅会对码头、港口和堤岸造成破坏，还会在洪水冲毁堤坝后，漫滩淹没和损坏沿海水产养殖、农田和房屋[1]。21世纪以来，每年因风暴潮所带来的经济损失，都高达人民币数百亿元，我国沿海风暴潮灾害影响出现日益严重的趋势[2–3]。

目前国内对风暴潮漫滩进行数值模拟模型的网格分辨率，最高可达到50 m左右。郑国诞等[4]利用MIKE21FM模块建立了适用于台州温岭市的高分辨率风暴潮漫滩数值模式，台州近海及陆上网格边长约50 m。傅赐福等[5]利用ADCIRC模型建立了一套非结构三角形网格，适用于滨海新区的高分辨率风暴潮漫滩数值模式，陆地区域分辨率达到50～80 m，对两次典型的温带风暴潮进行模拟。李勇等[6]对渤海湾西部风暴潮漫滩进行数值模拟，计算区域覆盖了整个东海海域，对普陀海域进行重点加密，其分辨率为30～200 m。Yang等[7]利用台风“天鸽”（2017年）和台风“山竹”（2018年）台风对澳门的淹没图，结合实测气象和波浪资料验证了一套完整的模拟台风期间潮波场耦合的数值模型，网格最高分辨率为20～100 m。

国内尚未有在地形复杂、街区密集情况下，基于街区尺度5～10 m分辨率的风暴潮漫滩数值模拟。但是随着城市建设步伐的加快，建筑建设鳞次栉比，街巷阡陌纵横紧密，在近岸和复杂的海岸地形区域中使用高分辨率网格，对于精细化的模拟和预报风暴潮漫滩都十分必要。这将是未来风暴潮预报精细化的趋势，对于制定城市的防灾减灾战略和完善风暴潮预警机制有着重要意义。

本文将利用ADCIRC模型，依据宁波市北仑区街区尺度的高分辨率高程数据，建立街区尺度风暴潮漫滩模型，北仑区陆地分辨率达到5～10 m。利用1211号台风“海葵”和5612号台风“温黛”对街区尺度风暴潮模型进行风暴潮模拟验证。以5612号台风“温黛”为例，在北仑区开展基于街区尺度5～10 m分辨率的风暴潮漫滩数值模拟，并考虑下垫面底摩擦的变化对风暴潮漫滩的影响。

2 街区尺度风暴潮漫滩模型的建立

2.1 ADCIRC模型介绍和参数设置

本文风暴潮模拟采用的是ADCIRC二维模型。其在球坐标系下通过基于垂直平均的原始连续方程和海水动量方程来求解自由表面起伏、二维流速等3个变量，即 ζ ,u,v。在球坐标系下海水的连续方程和海水原始动量方程为

式中，λ和φ为经度和纬度；ζ为从平均海平面起算的自由表面高度；H为海水水柱的总水深；U和V为深度平均的海水水平流速；R为地球的半径；η为牛顿引潮势；g为重力加速度；f为科氏参数；Ps为海水自由表面的大气压强；ρ0为海水密度；Dλ，Dφ为动量方程的水平扩散项；τbλ,τbφ为海底摩擦力分量；τsλ,τsφ海表面应力分量。

初始条件为： ζ =u=v=0；海岸边界条件为：Vn=0这里Vn为 岸边界的法向深度平均流速。在求解所需物理变量的过程中，空间采用有限元法离散，时间采用有限差分法，将连续方程和运动方程通过引入了一个时空加权参数进行结合求解[8–9]，提高了计算结果的稳定性。参数设置采用冷启动、球坐标、二维模式、混合底摩擦形式、考虑有限振幅项并采用干湿法、科氏参数基于β平面近似、考虑潮汐、时间步长为 0.02 s，满足计算稳定和收敛要求。由 K1、O1、P1Q1、N2、M2、S2、K2等 8 个分潮驱动计算。

2.1.1 气压场和风场模型

本文选用的台风风场模型为Holland（1980）[10]台风气压场分布公式为

式中，P(r,θ)为径向距离r和方位角 θ的函数，是距台风中心r的海表面气压值；Rmax是台风最大风速半径；台风中心气压为Pc；台风以外不受干扰的背景气压为Pn，设为 1 012 hPa。

风应力在ADCIRC中的表达式为

式中，ρa为空气密度；Cd为拖曳系数采用Garratt[11]公式

若Cd超过了0.003，则Cd=0.003。式（5）中第1个是代表海面上10 m处风速，第2个是海面上10 m处风速的大小。

2.1.2 底摩擦设置

底摩擦的混合形式与二次形式表达式相同，均为

Cf作为底摩擦系数是变化的。混合形式为

式中，Cfmin和Hbreak都是常数；λ和θ也是常数，λ =1/3，θ=10，因此Cf只随总水深H的变化而变化。二次形式为

式中，g为重力加速度；n为曼宁系数；H0为临界水深，即H0=Hbreak；λ=1/3。

2.1.3 干湿网格处理

风暴潮漫滩采用干湿网格法进行计算。此算法在连续方程之后和动量方程之前进行计算，位于一次时间步长的中间部分。

首先，在模式中设定一个为较小正数的最小水深Hmin，计算范围中总水深大于Hmin的某个格点，模式将该点判断为湿点。第二，当干点变为湿点后，底摩擦力和水位梯度相平衡产生定态速度U，模式要求设定的参数Umin（湿点最小速度）需小于它，为了约束邻近格点的水位，设定Umin的值。

2.1.4 风暴潮漫滩溢流的计算

水位高度超过堤坝的顶部时，水流将流过堤坝。单位时间内的流量表达式为

流量系数Cd的表达式为

式中，v1为前进流速；g为重力加速度；h1为高于堰顶的水位高度。

2.2 街区尺度风暴潮漫滩模型

宁波市北仑区地处浙江最东岸，受风暴潮和海浪等海洋灾害影响较大。在历史上受到台风风暴潮灾害影响次数较多，台风路径也有详细历史资料，并且在北仑区有高分辨率的高程数据支撑，因此我们选取北仑区作为实验区域。

2.2.1 岸线、水深以及高程数据

采用的岸线数据为1997–2001年经国家测绘局海岸线资料订正后，比例尺为1∶250 000的海岸线数据。宁波市北仑区以外的水深数据是通过范围5°S～52°N，99°～157°E 之间，分辨率为 2′×2′的水深数据，插值得到的。北仑区近海以及甬江口的水深数据采用海图及实测水深。本文使用的北仑地区数字高程模型（DEM）数据空间分辨率为5 m，比例尺为1∶10 000，这对于高分辨率的风暴潮漫滩模拟提供了重要的基础数据，选取范围为29.891°～ 29.962°N，121.811°～121.905°E，如图1所示。

图1 所选区域Fig. 1 Selected area

在选取的区域中，低层建筑多为二层楼房，高程高于5 m的不做处理，低于5 m的高程统一赋值为5 m。高层建筑物所在地高程高于10 m的不做处理，低于10 m的高程统一赋值为10 m。对高层和低层建筑进行划分如图2左图所示，右图为处理后的效果。

图2 高低层建筑的划分示意图Fig. 2 Division diagram of high and low buildings

2.2.2 网格的建立

本文所选区域选用SMS v8.1生成街区尺度非结构三角网格。考虑到为了较好地模拟大洋潮波、风暴潮的传播过程，减小开边界对模拟结果的影响，选取的计算范围包括东海、黄海和南海部分海域，经纬度为 20°～34°N 和 115°～131°E。所选的北仑区街区分布密集，城市建筑错综复杂，为达到街区尺度的分辨率，将北仑区网格分辨率设置为5～10 m，沿岸分辨率逐渐向外海增大，能够满足对重点地物的刻画。整套街区尺度网格包括978 783个三角形网格，共计502 063个节点，北仑区陆地占据了668 162个三角形网格。考虑到所选台风以及漫滩模拟的要求，不考虑堤坝。

图3为所选网格区域、北仑区附近海域以及北仑区局部放大的陆地网格示意图。由于陆地网格分辨率过高，为了满足CFL条件，时间步长定为0.02 s，但是这样使计算时间过长，因此将部分陆地高低层建筑的边界设置为固体边界，即图3中北仑区陆地网格的空白区域，减少网格数量以节省运算时间。

图3 大洋所选网格区域（a）、北仑区附近海域（b）以及北仑区局部放大的陆地网格（c）Fig. 3 The selected grid of ocean (a),the sea area near the Beilun District (b),and the land grid magnified locally in the Beilun District (c)

图4展示了所选区域、北仑区附近海域水深分布。图5为北仑区高程图，采用分辨率为5 m的北仑区高程数据进行插值，能进一步细致刻画北仑区的城市建筑和街区分布，清晰的分辨高低层建筑。

图4 大洋（a）和北仑区附近海域（b）水深分布Fig. 4 The distribution of water depth in the ocean (a) and nearby seas of Beilun District (b)

图5 陆地高程插值图Fig. 5 Interpolation map of land elevation

3 街区尺度风暴潮模拟验证

3.1 1211号台风“海葵”风暴潮过程验证

1211号强台风“海葵”于2012年8月3日在西北太平洋上生成，5日17时加强为强热带风暴，6日加强为台风，8日凌晨03时20分在浙江省象山县鹤浦镇登陆，登陆时中心气压为965 hPa，近中心风力为14级。8日16时在浙江省杭州市境内减弱为强热带风暴。分别选取镇海、吴淞、乍浦和定海4个验潮站。台风风场资料来自台风年鉴，台风路径和验潮站分布如图6所示。

图6 1211号台风“海葵”路径（a）和验潮站分布（b）Fig. 6 Track of No. 1211 Typhoon Haikui (a) and distribution of tide gauge stations (b)

采用台风与天文潮的非线性耦合，模拟1211号台风“海葵”风暴潮过程。选取时间8月6日0时至8月8日20时，将4个验潮站风暴增水的实测值（蓝色散点）与模拟值（红色曲线）对比（图7）。表1为4个验潮站最大增水模拟值与实测值相对误差分析。

图7 镇海、吴淞、乍浦和定海站风暴潮实测值和模拟值对比Fig. 7 Comparison of measured and simulated storm surge values at Zhenhai,Wusong,Zhapu and Dinghai stations

表1 4个验潮站最大增水实测值和模拟值相对误差分析统计（单位：m）Table 1 Statistical table of relative error analysis of the measured and simulated values of the maximum water increasing at 4 tide gauge stations (unit: m)

由图7可以看出，4个测站的增水趋势大体一致，对比最大增水的发生时刻，镇海站模拟与实测值的最大增水发生时刻基本一致，定海站模拟的最大增水推后1 h，乍浦站模拟值的最大增水时刻提前6 h，吴淞站提前7 h。结合表1来看，镇海、吴淞、乍浦和定海4个验潮站的相对误差都较小，说明模型较好地重现了此次台风风暴潮的发生过程。

3.2 5612号台风“温黛”风暴潮过程验证

1956年5612号台风“温黛”生成于菲律宾以东洋面，7月28日02时达到台风强度，30日14时台风强度达到最强（台风中心气压905 hPa，近中心最大风速90 m/s）。8月1日由石垣岛和冲绳之间进入东海，8月2日00时登陆浙江象山，台风中心最大风力达到16级，登陆时中心气压达到923 hPa，登陆后继续向西北方向移动，8月5日消失于陕西境内。选取镇海、吴淞、乍浦和高桥4个验潮站，验潮站都位于台风路径右侧。台风路径和验潮站分布如图8所示。

图8 5612号台风“温黛”路径（a）和验潮站位置分布（b）Fig. 8 Track of No. 5612 Typhoon Wanda (a) and location distribution of tide gauge stations (b)

采用台风与天文潮的非线性耦合，模拟5612号台风“温黛”风暴潮过程（图9）。将4个验潮站的最大增水实测值和模拟值进行相对误差分析，如表2所示。从图9和表2可以看出，各验潮站增水变化趋势基本一致，增水峰值时间基本同步。最大增水时间对比，镇海站提前2.5 h，乍浦站提前0.67 h，吴淞站推后1.75 h，高桥站推后1.17 h。结合台风路径分析，镇海、吴淞、乍浦和高桥站位于台风路径右侧，离台风登陆点位置较近，模拟值偏大一点，相对误差分别为9.13%、12.96%、5.6%和5.44%，平均误差为8.28%，模拟效果较好。总体来说，街区尺度风暴潮漫滩模型能较准确模拟5612号台风“温黛”风暴潮情况。

图9 镇海、吴淞、乍浦和高桥站风暴潮实测值与模拟值对比Fig. 9 Comparison of measured and simulated storm surge values at Zhenhai,Wusong,Zhapu and Gaoqiao stations

表2 4个验潮站最大增水实测值与模拟值相对误差分析统计表（单位：m）Table 2 Statistical table of relative error analysis between the measured and simulated values of the maximum water increasing at 4 tide gauge stations (unit: m)

4 街区尺度风暴潮漫滩模拟

5612号台风“温黛”登陆强度大，深入内陆深，浙江遭受台风损失最为严重。全省共有75个县（市）和不同程度受灾作物735万亩。此次风暴潮过程破坏了85万间房屋，超过4 900人被砸死、淹没或者触电身亡，超过15 000人受伤[12]。产生的严重风暴潮灾害影响是近几十年来罕见的，因此选取5612号台风“温黛”进行漫滩数值模拟。

4.1 径流的设置

甬江位于浙江东部，4 518 km2的流域面积，横贯宁波市区，在宁波市区东北部由奉化江、姚江两条江汇合后入海，中下游主要为平原河流。由于姚江下游建成姚江大闸，如图10所示，本文选取了澄浪堰和姚江大闸所处位置为姚江和奉化江集中入流的径流边界。

图10 甬江澄浪堰和姚江大闸径流边界位置分布Fig. 10 Runoff boundary location distribution of Chenglang Weir on Yongjiang River and Yaojiang Gate

根据李文杰和邵学强[13]关于甬江径流流量的资料整理，采用1983−2003年21年间年均尺度径流流量作为姚江大闸和澄浪堰边界的径流输入数据。澄浪堰和姚江大闸的年均流量分别为16.475×108m3和12.321×108m3，本文选用这两个数据作为径流边界流量，并把径流流量当作为常数。

4.2 总潮位数值模拟

利用已建立的街区尺度风暴潮漫滩模型模拟5612号台风“温黛”影响下，北仑区附近验潮站总潮位的情况，并以此来验证5612号台风“温黛”过程中，所建立的街区尺度风暴潮漫滩模型对北仑区风暴潮漫滩模拟的准确性。

将镇海、吴淞、乍浦和高桥4个验潮站总潮位实测值与模拟值进行对比（图11）。由图11可以看出，对于各站总潮位的模拟，趋势与实测数据大体一致，最高潮位时间稍有偏差。将镇海、吴淞、乍浦、高桥4个验潮站模拟的最高潮位和实测最高潮位进行对比，发现各站绝对误差为0.03 m、0.34 m、0.79 m、0.36 m，各站相对误差为1.6%、14.5%、17.97%、15.28%。从模拟结果以及误差分析看，镇海站模拟的误差较小，吴淞、乍浦和高桥的误差稍大，是由于吴淞、乍浦和高桥3个站位于台风的右侧，模拟的最高潮位相对于实测值偏大一点。因此从总潮位的趋势看，建立的街区尺度模型可以较好地模拟北仑区的最高潮位。

图11 镇海、吴淞、乍浦和高桥站总潮位实测值与模拟值对比Fig. 11 Comparison of measured and simulated total tide level at Zhenhai,Wusong,Zhapu and Gaoqiao stations

4.3 漫滩数值模拟

在5612号台风“温黛”影响期间，浙江沿海产生特大海潮，象山县的最高潮位为4.7 m。据不完全统计，浙江、上海、安徽和河南省共影响作物6 946万亩，摧毁220万户房屋。但是当时没有对漫滩淹没范围和水深进行系统的统计，缺乏当年实测数据与模拟值对比。所以在风暴增水和总潮位模拟较好的效果下，认为此漫滩模拟结果有一定合理性。

在不考虑海堤的情况下，利用已建立的街区尺度风暴潮漫滩模型对5612号台风“温黛”作用下北仑区的漫滩情况进行模拟，其结果如图12所示。

图12 5612号台风“温黛”漫滩模拟的最大淹没范围Fig. 12 Maximum inundation area simulated by No. 5612 Typhoon Wanda

从图12中可以看出，在5612号台风“温黛”作用下，北仑区选定区域漫滩面积约为628.9 km2，淹没了近岸大部分区域，淹没水深大部分在0.1～0.4 m，最大水深在2.0 m左右。水流绕过高低层建筑，淹没了建筑周围的街道，更加清晰地看到淹没区域的具体情形。在边界处淹没水深较大，是因为此处为陆地边界，从而导致水流在这附近区域堆积。在街区尺度网格下，可以直观清晰地看到水势漫延趋势，以及在高低层建筑物的阻挡下水势的走向，体现出街区尺度网格的优势。

为了更直观地体现5612号台风“温黛”漫滩过程对北仑区街区和建筑的淹没情况，利用GIS三维可视化工具，在数字地面模型基础上叠加高分辨率遥感影像制作了此次漫滩过程的三维立体图。选取时间是8月1日14时至8月2日10时，时间间隔为4 h。从图13可以看出此次漫滩过程水深逐渐加深和降低的细致变化，对北仑区沿海街道和建筑的淹没比以往的漫滩模拟更加精细。可以看出水流在街区和建筑错综复杂分布时的流动情况，以及不同建筑物附近水深的变化情况，体现了街区尺度漫滩模拟的显著特点。

图13 5612号台风“温黛”漫滩过程三维图Fig. 13 Three-dimensional diagram of the inundation process of No. 5612 Typhoon Wanda

5 下垫面底摩擦变化对风暴潮漫滩的影响

不同下垫面具有不同的底摩擦系数，在海水的流动过程中影响着海水的流速和流向，这都直接影响漫滩模拟的结果[14–16]。

将下垫面的不同地物分类，修改控制方程中底摩擦系数的表达形式。此次采用二次形式的底摩擦，将下垫面底摩擦的变化考虑在ADCIRC模型中，底摩擦项为

式中，g为重力加速度；n为曼宁系数；H0为临界水深，即H0=Hbreak。

利用北仑区遥感图像对北仑区地物进行分类。根据美国NLCD数据的分类标准，我们选取21、22、23、31、43、51、83等 7类地物，不同地物的曼宁系数

如表3所示。即高、低密度居民区、商业区、裸露的岩石和沙地、四季常青和随季节变化的两种森林的混合、灌木丛和农作物稀少地，图14为所选区域曼宁系数n的分布示意图。

图14 北仑区曼宁系数n示意图Fig. 14 Schematic diagram of Manning-n in the Beilun District

表3 7类地物的曼宁系数（n）设置Table 3 Manning-n settings for seven types of features

在考虑下垫面底摩擦变化的条件下，模拟北仑区在5612号台风“温黛”作用下最大水深淹没范围，如图15所示。淹没面积达到494.37 km2，对比未考虑下垫面底摩擦情况的模拟实验结果，淹没面积减少了21.4%。沿岸大部分区域，最大淹没水深在0.2 m左右，部分区域达到0.5 m，沿岸淹没较深的局部区域可以达到1 m以上。

图15 考虑下垫面底摩擦变化后，北仑区最大淹没范围Fig. 15 The maximum submerged area of Beilun District after considering the change of the underlying surface friction

图16为对比了考虑下垫面底摩擦变化前后，北仑区淹没水深的变化，即考虑下垫面变化后的最大淹没水深减去未考虑下垫面变化时模拟的最大淹没水深的差值。可以看出在考虑下垫面底摩擦变化之后，大部分区域水深降低0.1～0.2 m。其中大部分水深降低区域集中在沿岸的低密度居民区，降低0.1 m左右。在农作物、商业区局部区域，水深降低了0.5 m左右。其中水深降低最大的区域是在北仑区内的两条河流两侧。因此可以看出，在考虑下垫面的底摩擦变化之后，其对最大淹没水深有降低作用，在不同的区域水深降低程度不同。

图16 考虑下垫面底摩擦变化前后，北仑区最大淹没水深差值Fig. 16 Maximum submerged depth difference in the Beilun District before and after considering the underlying surface friction

图17为对比考虑底摩擦前后，漫滩面积的变化。蓝色为面积减少区域，约减少138.76 km2，该区域多为低密度居民区。红色为面积增加区域，约增加15.61 km2，该区域为岩石区和低密度居民区。总体而言，淹没面积减少。极小的面积增加区域可能是由于在考虑下垫面底摩擦变化之后，街区细致刻画使得水流流速变化，致使部分区域漫滩面积增加。因此在考虑底摩擦情况后，淹没面积总体减少，下垫面的底摩擦变化对淹没面积有减少作用。

图17 考虑下垫面底摩擦变化前后，淹没区域变化Fig. 17 The change of the submerged area before and after considering the change of the underlying surface friction

6 结论

本文采用5～10 m的街区尺度网格，依据宁波市北仑区高分辨率的街区尺度高程数据，细致刻画街区道路和建筑分布，建立街区尺度风暴潮漫滩模型，开展基于街区尺度的风暴潮漫滩数值模拟，并考虑下垫面底摩擦的变化对漫滩的影响。结论如下：

（1）选取对北仑区影响严重的5612号台风“温黛”，在考虑甬江径流和不考虑堤坝的条件下，街区尺度风暴潮漫滩模型可以较好地模拟风暴潮和总潮位。依托街区尺度网格，可以看到水流沿着街道向城市内部漫延的情况以及在高低层建筑物的阻挡下水势的走向，体现出高分辨率街区尺度网格的优势。

（2）考虑下垫面底摩擦变化对漫滩的影响，漫滩淹没面积达到494.37 km2，对比未考虑底摩擦情况的模拟结果，淹没面积减少了21.4%。在街区尺度网格下，细致展示了淹没水深和淹没面积在不同地物分布情况下的变化。大部分水深降低区域集中在沿岸的低密度居民区，降低0.1 m左右。在农作物、商业区局部区域，水深降低了0.5 m左右。其中水深降低最大的区域集中在北仑区内的两条河流两侧。淹没面积减少的区域多为低密度居民区。总体而言，考虑下垫面底摩擦变化后，淹没水深降低，淹没面积减少。