裂流危险性的数值预报方法及其在三亚大东海浴场的应用

杨怀玮，原野，高义，邢闯，高志一

（1.国家海洋环境预报中心，北京 100081；2.厦门大学，福建 厦门 361102；3.自然资源部海洋灾害预报技术研究重点实验室，北京 100081）

1 引言

裂流是传播至岸边的波浪在地形影响下发生变形和破碎，并在波浪辐射应力和近岸壅水的共同作用下，形成的一股向海的快速的条带状海流[1]。完整的裂流系统一般由裂流头、裂流颈、补充流和向岸流组成，是破波带波生环流的重要组成部分。裂流速度通常在0.3～1 m/s，最快可达1～2 m/s，离岸距离可达数百米[2-3]，具有流速快、流辐窄以及方向与海岸近乎垂直等特点[4]。裂流的存在对近岸物质交换、悬浮物及泥沙输运[5]等方面均有显著意义，裂流可以在很短的时间内将游泳者带至离岸很远的区域，导致游泳者在惊慌失措中溺水[6-7]。

在美国，裂流引起的事故在近海事故中占60%，每年有超过100 人因裂流溺亡[8]；在澳大利亚，近90%的海滩救援和溺水事故都与裂流有关[9-10]，每年都有40～50 起裂流相关的溺水事故[11]；在英国，裂流灾害事件达到海滩所有事件的2/3；2012年8月，韩国海云台浴场突发强烈裂流，142名游客被卷入水中[12]；除此之外，裂流也对印度[13]、巴西[14]和新西兰[15]等海滩游客的生命财产安全造成极大威胁。在我国，根据自然资源部海洋减灾中心初步开展的全国裂流风险调查和评估成果，福建、广东和海南等省份的部分海水浴场存在不同程度的裂流风险，有必要推动裂流观测预报和减灾研究工作。

裂流研究一直是国内外学者共同关注的热点问题。本文建立一种裂流危险性数值预报方法，并对我国三亚大东海浴场的裂流危险性进行讨论。

2 国内外研究现状

2.1 数值模拟

当前模拟裂流的数值手段主要分为两类：一类是波浪时均模型，另一类是时域波动模型（或波分辨模型）[16]。前者通过波浪谱模型（如海浪数值模式SWAN（Simulating Waves Nearshore）等）计算得到波高分布场和波浪辐射应力，进而对二维浅水方程（或三维水动力模型）添加波浪动量通量（辐射应力项）来模拟波致平均流场。通过在模型中考虑次重力波和波群效应，也可以对裂流的时空变化特征进行模拟。但此类方法在模拟波浪非线性效应和波波相互作用等方面有所欠缺，对小尺度的波生流现象模拟效果不佳。时域波动模型通过构建不同的造波器，采用非线性高阶频散方程对该造波器强迫下的二维波浪场的传播、波波相互作用和破碎等过程进行直接模拟。通过对波浪场和流速场进行1个至数个波周期的平均，可以得到破波带壅水和波生流的水平分布特征。最典型的时域波动模型包 括FUNWAVE-TVD（Total Variation Diminishing version of the FUlly Nonlinear Boussinesq wave model）和SWASH（Simulating WAves till SHore）等[17]。

2.2 裂流预报

Lushine[18]将美国佛罗里达州东南海滩的裂流溺亡和救援数据与风向、风速、浪高和潮时进行统计分析，建立了LURCS（LUshine Rip Current Scale）方法。Lascody[19]和Engle 等[20]在此基础上做了改进，降低甚至去除风的比重，增加对涌浪、波向和潮位的考虑，大大提高了预报的客观性。2013 年，Kumar 等[13]针对印度沿岸开发了一种新的经验预测方法CAP-LURCS（Coastal Andhra Pradesh LURCS），该方法能够进行每小时裂流自动预报，且达到了更高的预报精度，但无法对裂流的发生位置进行预测，可靠性依然有待提高。

需要注意的是，上述经验预测指标体系均不具有普适性，在预报应用中应结合海滩的实际地貌、水深特征、水位气象特征和裂流救援数据进行修订，因此开展裂流数值预报是今后发展的重要方向。目前，国际上部分国家已把裂流预报作为日常预报项目之一，如美国气象局把全美重要海滩的裂流风险预报作为冲浪区预报的一部分；韩国气象局也建立了裂流经验和数值预报一体的预报系统。此外，澳大利亚和印度等国家也开展了裂流预报业务化研究和应用。自然资源部近年来开展了裂流风险调查评估工作，但尚未针对重点区域建立有效的裂流预报系统。

2.3 裂流风险等级评价

裂流的生成与近岸破波带的地形变化特征息息相关，因此通过对海滩地形的变化趋势进行地貌分类成为近年来学者的研究重点。Dean 常数Ω可以作为裂流发生的一个重要评估指标[21]：

式中，H是破碎波高；W是泥沙沉降速度；T是入射波周期。李志强等[22-24]在基于地形动力学的裂流安全性评价方面做了很多工作。

基于水动力数值分析判定裂流风险等级也是一种重要方法。基于滨海旅游地段的实际水深地形，进行裂流及波流相互作用等水动力情况的精细化模拟推衍，入射波浪采用标准JONSWAP波浪谱。在满足模拟的精度、稳定及时长等要求下，结合裂流最大流速、长度及宽度的模拟结果进行等级判定。还可针对遥感影像直接进行解译分析，分辨率需优于1 m，重点选取弧形海滩、周期性滩角发育岸段和海岬等裂流高发区域，对裂流及地形特征进行风险指数判定。此外，也可利用无人机勘测裂流或水底沙槽地形，探测精度小于0.5 m，生成数字高程模型（Digital Elevation Model，DEM）；彭石等[25]和张尧等[26]对环保染料追踪裂流进行了研究。

自然资源部在2020年发布的《滨海旅游区裂流灾害风险排查技术规程》（HY/T 0298—2020）中提出，在判定目标岸段裂流灾害风险分析时，至少需要采用海滩地形动力学评价、水动力数值分析和遥感影像解析中两种及以上的方法进行初步判断。

3 裂流数值预报系统

建立裂流数值预报系统是防范裂流灾害风险的主要发展方向。与其他海洋灾害数值预报技术相比，由于裂流具有空间范围小、发生发展瞬变性、预测业务链条较长以及对计算资源需求较大等因素，开展裂流危险性数值预报具有较大难度。本文提出的裂流数值预报系统及其在三亚大东海浴场的应用技术路线如图1所示。在数值预报业务链条上，采用欧洲中期天气预报中心（European Centre for Medium-range Weather Forecasts，ECMWF）西北太平洋风场数值预报作为初始场，耦合中国近海WRF（Weather Research and Forecasting）中尺度风场，对研究区域所在海域开展72 h 海面风预报；采用中尺度风场作为驱动场，基于具有高性能和高效率的海浪数值预报模型，建立南海中北部海浪预报系统的72 h 海浪预报；实时提取预报点外海二维波浪谱，采用FUNWAVE-GPU 建立水平分辨率1 m 的波动及波生流预报模型；结合研究区域破波带地形特征以及风浪流要素，建立针对性的裂流危险性预报指数；最终建立区域大气强迫-区域波浪谱-近岸裂流为业务链条的嵌套数值预报模型。

图1 裂流预报系统技术路线Fig.1 The technical process of the rip current forecasting system

3.1 试点预报区域

本文选取位于三亚南侧的大东海浴场作为试点预报区域。大东海海滩长约2 km，位于兔尾岭与鹿回头半岛间，面向正南，为公共开放海滩。据统计，该海滩每年至少需要营救溺水人员上百人[23]。在Google Earth 上截取了2018 年8 月7 日的卫星图（见图2），海滩沿岸的沙坝-沟槽地型清晰可见。由于海底沙坝的存在，入射波浪在破波带沿岸方向存在显著的破碎不均一现象，这点可通过沙坝上方波浪破碎产生白色泡沫观察得到。海底沙坝之间的深色区域为裂流槽（Rip Channel），有间歇性裂流发生。海湾两侧和岬角水深较浅，波浪作用显著，底质主要为沙砾和礁石。西侧水面下有珊瑚礁分布[27]，东侧主要为基岩地形。通过查阅现场踏勘资料和历史卫片，大潮期海湾两侧会在退潮时出露海面。

本文基于截取自Google Earth 的卫星图片，根据水陆光谱差异提取各成像时刻的水边线，基于潮滩地形反演模型[28]得到高分辨率潮间带地形。依据海底沙坝-沟槽地形在遥感影像中的反射率（表现为水色不同），将其平均划分为5 个等级，根据划分等级提取相应等值线，从陆向海分别以0.2 m 为间隔赋高程值，水下区域地形通过在0 m等深线和5 m及10 m 等深线间等间距加密，然后通过空间插值得到水深地形数据。本文研究的大东海海域具体地形如图2b。

图2 大东海地形及水深Fig.2 Topography and bathymetry of Dadonghai

根据自然资源部三亚海洋观测站2018年水位、水温和波浪连续观测数据，结合国家海洋环境预报中心30 a海浪再分析数据集（见图3和图4），可以看出：大东海海域四季水温适宜，是游客冬季旅游的最佳去处之一；该区域为不正规半日潮和全日潮的混合潮型，以不规则半日潮为主，潮差一般小于2 m，为弱潮能区；受盛行季风和湾口朝向影响，有效波高和波周期呈现较为显著的季节变化，夏季盛行西南季风，入射的优势浪向为南向和东南向，有效波高和波周期较大；在正常的天气条件下，受南海开阔海域入射涌浪影响，日平均有效波高偶尔会达到1.0 m 以上；夏季受过境西行或西北行台风影响，有效波高呈现峰值，日均有效波高可达1.5 m 以上，且涌浪可持续一周以上。

图3 大东海浴场相关水文再分析数据统计Fig.3 Statistics of relevant hydrological reanalysis data in Dadonghai beach

图4 大东海浴场外海波玫瑰图Fig.4 Rose map of the wave of Dadonghai beach

3.2 裂流数值模拟

本文主要探讨裂流的数值预报，关于中国近海海面风数值预报模型和南海中北部大区域海浪数值预报模型不在此赘述，详情可参考中国海洋预报网（网址：http://www.oceanguide.org.cn）。南海中北部大区域波浪谱数值模型可对未来3 d 大东海浴场外海的有效波高、波向和波周期进行预测，并提供二维波浪谱参数，为嵌套大东海浴场裂流数值模型提供波浪边界条件。

裂流数值模拟主要基于时域（波分辨）波动模型FUNWAVE 开展。首先，通过历史高分辨率卫星图片基于水色进行预报区域破波带地形反演，反演地形包括沙坝-沟槽型、横向沙坝型和缓坡型等，并与现场观测水深进行校验；其次，采用经GPU 加速的FUNWAVE-GPU 模型建立水平分辨率1 m 的波动及波生流数值模型。通过施加不同潮位、不同入射方向和有效波高的不规则波，重点考察不同水动力条件下裂流的产生和时空分布特征。

FUNWAVE 模型最初基于Kirby 等[29]建立的完全非线性Boussinesq 方程，Shi 等[30]采用基于总变差减少（Total Variation Diminishing，TVD）的高精度数值格式对FUNWAVE 进行重构，提出了FUNWAVE-TVD 模型，该模型在数值模拟稳定性方面有很大的改善。其控制方程如下：

式中，∇是水平梯度算子；M是水平体积分量，表示为：

式中，总水深H=h+η，h代表静止水深，η代表波高；uα是位于参考高度zα处的流速；代表水平速度场的平均深度贡献，是一个二阶修正量。该模型可以在包含波的非线性和色散性的条件下应用，在模拟波传播过程中各因素的相互作用方面非常出色。此类模型在开展破波带海浪和裂流模拟中应用广泛，但也存在计算量大的弱点。Yuan 等[31]开发了FUNWAVE-GPU，可以通过GPU 并行加速对二维波浪场进行快速模拟，一定程度上克服了Boussinesq方程计算效率低的问题。

大东海浴场裂流数值模型的水平网格为1 678×659，地形逆时针旋转90°，边界条件采用JONSWAP二维波浪谱构建不规则入射波场，底摩擦系数设置为0.002 5，模型模拟时长为3 000 s。FUNWAVETVD 对波生流的数值预报所需时间较长，为保证裂流在规定时间内的可预报性，采用FUNWAVE-GPU进行建模。与36 核的计算节点相比，单GPU 卡的并行加速比可达到1:5，模拟用时小于30 min。

3.3 不同水动力条件产生裂流的数值试验

选取不同潮位、不同入射波浪有效波高及周期条件的裂流，对其产生和发育情况进行模拟，并选取点A（见图2b）的流速时间序列进行统计分析。不同数值试验的条件设置见表1。表中潮位以平均海平面为基准，入射角度以垂直向岸方向逆时针为正。

表1 数值模拟的不同水动力条件设置Tab.1 Setting of different hydrodynamic conditions for numerical simulations

3.4 裂流危险性指数

在3.3 数值模拟结果的基础上，选取模型稳定后（600 s）的模拟结果进行处理。以5 m为水平间隔进行网格划分，并对各网格中离岸流（定义为流向在45°～135°之间）的最大速度及此速度的持续时间进行统计和赋分，其中，持续时间定义为离岸流出现的时间（t）与总模拟时长（tm）的比值。对照表2，将两指数相乘的结果划分为4个危险性等级进行判定。结果在[1-4）为四级，用蓝色表示；结果在[4-8）为三级，用黄色表示；结果在[8-12）为二级，用橙色表示；结果在[12-16]为一级，用红色表示。应当注意，对研究区域的危险性指数判定应根据研究区域的不同情况进行调整，本文提出的判定方法仅针对3.3中的大东海浴场的数值模拟结果。

表2 裂流危险性等级判定Tab.2 Rip hazard level

4 结果与分析

4.1 不同水动力条件下的裂流强度变化特征

序号1—8（见表1）的数值试验结果如图5。由图可见由于波浪破碎效应，海岸两侧会产生较强的沿岸流。对于半封闭或者口袋型海湾，入射波浪在湾口两侧率先破碎，形成波生沿岸流，泥沙向湾顶净输运，这也导致海湾两侧主要为基岩底质分布。从图5a—d 可以看到y=950 m、800 m 和750 m 处均产生了明显的裂流。随着入射波高从0.7 m 增大至1.8 m，裂流的空间尺度和流速显著增加。如图5b所示，位于y=950 m 和800 m 处的离岸流发生汇聚。红框中选取典型裂流槽并放大置于图5e—l，可见裂流呈现复杂的空间分布特征，离岸流两侧有小尺度涡旋对分布，伴随着涡旋对的发展离岸流的流向左右摆动。通过对比图5f、5i和5j图中的裂流特征，我们发现不同入射波周期导致同一裂流槽中的裂流传播方向产生明显差异，入射波周期增大使得裂流强度增强。如图5k—l 所示，对于斜向入射波浪，波生流的沿岸流动分量显著增大。

图5 不同水动力条件下的模拟结果Fig.5 Simulation results under different hydrodynamic conditions

由图6 波生流时间序列可知，裂流具有不稳定特征，且呈现显著的周期性波动。这与破波带地形、入射不规则波波群效应以及波生涡旋的产生具有密切联系。一般情况下，普通游泳者的速度在0.1～0.3 m/s，当裂流流速超过0.3 m/s时就会有将游泳者带入海中或使其陷入流涡的危险。裂流流速瞬变性特征明显，可以在短时间内达到或超过0.3 m/s，这一特点很容易使溺水的游泳者产生心理恐慌。

图6 点A处波生流变化时间序列Fig.6 Time series graphs of wave-induced current changes at point A

根据图7可以看出同一入射波浪周期下裂流传播方向大致相同，与沙坝之间的裂流槽开口朝向一致。入射波浪的持续增大使得裂流强度不断增强，方向更集中。结合图7b、7g 和7h 观察，潮位升高会使裂流强度变弱，离岸方向分散。由于水深变深，波浪即便在沙坝上方也不易破碎，产生的裂流较弱。这也证实了低潮位会使裂流的发生概率大大增加。

图7 点A处波生流玫瑰图Fig.7 Wave-induced current rose diagrams at point A

4.2 裂流危险性等级预报

基于本文提出的裂流危险性等级划分标准，对2018 年8 月7 日的裂流危险性进行后报，绘制出我国首幅裂流危险性等级后报图，如图8 所示。总体上看，大东海浴场的沙坝-沟槽地形有利于产生较强的裂流，对比左侧的地形图可知，在较为发育的裂流槽向海一侧的40～200 m 距离内，裂流危险性达到三级以上。由于海湾西侧（图中下方）较强的沿岸流在流速剪切效应下向外海一侧偏转，与裂流槽产生的裂流汇聚，产生较大范围的离岸流，更易产生裂流危害。图中B 点所处位置均匀分布着3 个裂流槽，相比A 点处分布单个沟槽的地形来说，裂流危险区向海一侧延伸更远，需引起注意。C 点处裂流危险等级为二级，向海一侧延伸的距离约为50 m。上述特征说明在裂流灾害防御中需要更加重视沿岸分布密集的裂流槽。

图8 裂流危险性示意图Fig.8 Hazard level of rip current

5 结论

我国的裂流研究工作已经引起部分学者的关注，但相关工作还局限于实验室物理模型研究、海滩裂流安全性指数评价以及理想条件下裂流的数值模拟等阶段。本文首次提出了裂流危险性数值预报的技术路线并在三亚大东海浴场进行了后报试验。由于裂流数值预报的业务链条较长，本文采用基于GPU 并行加速的FUNWAVE-GPU 模型进行高分辨率裂流数值模拟，大大削减了数值预报耗时。模拟结果表明裂流对于不同水动力条件及地形等因素均具有较强的敏感性，入射波高及入射波周期的增大会使裂流空间尺度和流速明显增加；不同入射波周期会使裂流传播方向明显不同；入射角度的增加会使受地形影响产生的沿岸流特征明显；低潮位时更易产生裂流。

未来将对南海中北部大区域海浪谱预报模型进行GPU 加速，可进一步缩短裂流数值预报耗时。同时，需要注意破波带裂流与地形存在相互作用，沙坝和沟槽的演变对裂流数值预报的准确性有直接影响，因此须采用高精度卫星图片定期对破波带地形进行实时监测和反演。下一步将在更多的裂流高风险区开展裂流综合性现场观测，并建立相应的裂流危险性数值预报系统，系统地开展裂流危险性指数预报工作。