冯 鑫 王 刚 张甲波 刘会欣
(1.河北省海洋地质资源调查中心,秦皇岛 066000;2.河北省海洋岸线生态修复与智慧海洋监测工程研究中心,秦皇岛 066000)
随着人类认知水平的不断提高和科学技术快速发展,人类对自然资源的开发利用程度正逐渐加深[1]。近年来,我国海洋经济增长迅猛,成为我国国民经济新的增长极,并初步形成了北部、东部、南部三大海洋经济圈[2]。但是,海洋在给人类带来福祉的同时,也带来了风暴潮、海冰等灾害,其中风暴潮是我国水文气象灾害中最严重的海洋灾害。根据河北省自然资源厅2020年3月发布的《2019年河北省海洋灾害公报》,2019年河北省沿海共发生风暴潮过程2次,其中1次台风风暴潮和1次温带风暴潮,造成沿海地区直接经济损失3.34亿元。尽管风暴潮造成的人员伤亡和财产损失在总体上都具有不可避免性,但人类可通过自己的努力使灾害与灾害的损失得到减轻,或者有效地缓和灾害损失上升的趋势[3]。风暴潮数值模拟是一种依靠计算机手段以达到灾前有效防治的方法,可有效减轻海洋灾害损失。因此,开展台风风暴潮数值模拟研究工作是至关重要的。目前,国内外已有相关学者对台风风暴潮的数值模拟进行了研究。
国际上自20世纪70年代起开展了风暴潮灾害数值预报预警技术的研究[4],并于90年代形成较为完善的预警预报系统[5-6]。而我国在风暴潮预报和计算方面的研究较国际上起步较晚,国家海洋局于1974年5月在厦门召开了中国首次风暴潮经验与预报交流会,并首次开展风暴潮预报工作[7]。经过近些年的不懈努力,我国在风暴潮数值模拟研究方面得到了迅速地发展,其中,李健等[8]利用实测资料对比和FVCOM(an Unstructured Grid,Finite-Volume Coastal Ocean Model)数值模拟等方法,研究了201909 号台风“利奇马”风暴潮在渤海的演变规律。朱婧等[9]基于FVCOM 风暴潮模式,利用重建的台风风场资料,模拟了201614号台风“莫兰蒂”过程中厦门湾及其附近海域的风暴潮。王雪迎等[10]基于ADCIRC(an ADvanced CIRCulation model for oceanic,coastal and estuarine waters)和SWAN(Simulating WAves Nearshore)模型构建风暴潮和台风浪耦合模型,研究了长江口台风期间波致增水空间分布特征研究。
曹妃甸区地处环渤海中心地带、唐山南部,毗邻京津两大城市,是京津冀协同发展的战略核心区,同样又是渤海风暴潮严重地区之一。根据多年统计资料,渤海特大风暴潮发生概率大约10 年一次。近年来,201723 号台风“达维”、201810号台风“安比”、台风“利奇马”均对曹妃甸区产生较大影响。对此,本文以曹妃甸区为例,采用SMS(Surface Water Modeling System)软件开展前处理工作,实现黄渤海大尺度和曹妃甸区精细化非结构网格的绘制,选取近岸海洋数值模型ADCIRC 水动力模块,完成风暴潮漫滩和漫堤数值模拟,运用GIS(Geographic Information System)系统及数字可视化技术绘制了曹妃甸区淹没水深分布、危险性等级分布和风险等级分布图,并提出综合防潮对策建议以提高曹妃甸区风暴潮灾害防治能力,最大限度减轻风暴潮灾害带来的人员伤亡和财产损失。
ADCIRC模型是由美国北卡罗来纳大学海洋研究所联合多所大学共同研发的有限元海洋模型,具有高效、准确、守恒和自动捕捉间断流动等特点,不仅可用于模拟常规的河道和近海浅水流动,如天文潮、风暴潮、漫堤淹没等,还可以模拟涌潮、海啸和溃坝等特殊的水流现象。适用于海岸地区的平面二维水动力数值模拟,如海啸、天文潮和风暴潮等,可应用于水利工程设计及规划、复杂条件下的水流计算、治理规划和洪水风暴潮防灾减灾等。
控制方程包括1个连续性方程和2个动量方程:
式中:U、V分别为x、y方向上的垂向平均分量;H为总水深;f为科氏力系数;ps为表面大气压力;ρ0为海水密度;g为重力加速度;(η+γ)为牛顿引潮力势和固体潮作用项;τsx和τsy分别为表面风应力和波浪辐射应力梯度的x和y方向分量;τbx和τby、Dx和Dy、Bx和By分别代表底部切应力、扩散项以及斜压梯度的x和y方向分量。
本文按照《海洋灾害风险评估和区划技术导则第1部分:风暴潮》[11]的要求,选择历史上对曹妃甸区影响最严重、风暴增水最显著的“7203”台风作为典型台风过程,该台风于1972年7月26日(农历六月十六)15时登陆山东省荣成沿海地区,登陆时其最大风力达到10级,中心气压是971 hPa,穿过渤海海峡后,又于1972年7月27日上午7时到8时,在天津市塘沽地区登陆,此次登陆时最大风力为7级,中心气压980 hPa,在曹妃甸附近海域其中心气压可达到975 hPa,台风相关信息见表1。
表1“7203”台风路径及中心气压
在风暴潮数值模型建立之前,台风场和气压场的计算是一个重要的环节。台风模型主要包括气压模型、环形风速模型、移动风速模型以及其他要素,我国常用的台风模型主要有Fujita-Takahashi 模型、Holland 模型以及Jelesnianski模型。在目前开展的研究中,石洪源等[12]探究了不同台风模型在南海的适用性,结果表明南海地区采用Jelesnianski 模型得出的风场和实际风场吻合度最好。魏晓宇等[13]基于Jelesnianski 模型开展风暴潮预报方法在珠海及粤西海域的应用,结果表明预报准确率分别达95.0%和87.6%。对此,本文在台风场和气压场的模拟计算中选用Jelesnianski模型,计算公式如下:
式中:R为最大风速半径;r为计算点到台风中心的距离;V0为台风移动速度;WR为台风移动速度;(x,y)、(xc,yc)分别为计算点坐标和台风中心坐标;θ为流入角(计算中当r≤R时θ取10°,当r>1.2R时θ取25°,其余的θ在10°和25°之间线性内插);P0为台风中心气压,P∞为无穷远处的大气压(计算中取1 010 hPa);β为台风风速距离衰减系数。
WR为台风最大风速,使用ATKINSON-HOLLIDY(1977)[14]提出的风-压关系式来计算。
风暴潮数值模型计算域的选取主要考虑3个方面:首先,计算范围较大,囊括台风影响区域,边界远离研究区域,避免边界扰动对研究区域的水流特征产生影响;其次,预估风暴潮漫滩或漫堤后水流的最大可能淹没范围,确定有阻水作用建筑物的位置及高程,如海塘和高速公路等;最后,确定周边重点保护对象,如核电站、重要化工厂等。
本文重点关注流经曹妃甸区的主要河道及其邻近海域,其中,流经曹妃甸区周边的自然河流及人工开挖的排水主干渠有9条,包括8条自然河流和1条人工河流,人工河东西贯穿溯河和小青龙河,除人工河外,其他河流由北向南穿境入海,自西向东依次是沙河、黑沿子排干、双龙河、一排干、小青龙河、溯河及滦南第二泄洪道。以上河流均发源于上游各县,河流总面积1 596 km2,其中域内(含唐海部分农场及滦南柳赞镇部分,下同)流域面积169.2 km2,域内河道总长度198.1 km,最大泄洪流量657.9 m3/s,年径流量1.309亿m3;外海边界设在台风第二警戒线附近,囊括所有对曹妃甸区有影响的台风路径,包括渤海、黄海等,考虑双龙河口、小青龙河口、溯河口等沿海主要河口地区。陆地部分覆盖整个曹妃甸区域,包括唐海镇、柳赞镇等及邻近乡镇,考虑环城高速、唐曹高速、省道和县道等主要阻水建筑。曹妃甸主要入海河流位置分布见图1。
图1 曹妃甸主要入海河流位置分布图
整个计算域采用三角形网格剖分,总计453 860个单元和227 509个节点。网格尺寸由外海约30 000 m逐渐过渡至莱州-黄骅附近海域约800 m,曹妃甸区平均网格边长约300 m,重点加密所关注的主要河道及其邻近海域区域,近岸网格分辨率为60~300 m。外海边界条件由TPXO 7.2 全球潮波模型提供,所有数据均采用1985国家高程基准,投影坐标采用经纬度。计算网格及高程分布如图2所示。
图2 计算区域和网格图
根据资料收集情况,在黄渤海主要潮位站点中选取连云港、日照、烟台、黄骅港、塘沽、京唐港、秦皇岛等12个潮位站,验证模型在黄渤海间传播的可靠性。根据典型台风发生时间选取2011年7月黄渤海各潮位站的潮位数值,将各潮位站实测资料进行调和分析作为天文潮实测值,分析模型天文潮计算值与实测值间的误差以评价模型的优劣。
为了定量评价水动力模型模拟结果的优劣,需要寻求一个评价标准来进行衡量。对水动力模型计算结果的评价,传统的做法是通过计算模拟值和实测值的相关系数和标准差来进行。但相关系数只能应用于线性模型,有一定的局限性。本文主要采用目前常用的Wilmott 统计学方法对模型的模拟结果进行评价,Wilmott 统计学方法考虑了实测值与实测平均值的偏差、模型计算值和实测平均值的偏差这两者的相关程度,其计算方法为:
计算所得到的skill代表了实测值D与实测平均值-D的偏差、模型计算值M和实测平均值-D的偏差这两者的相关程度,计算所得的skill的范围在0~1。skill=1 时代表模型计算值和实测值之间完全相符,skill>0.65 时结果为极好,0.5<skill<0.65 为非常好,0.2<skill<0.5 为好,小于skill<0.2为差,0代表模型计算值和实测值之间完全不符。水动力模型效率系数见表2。
由表2可以看出,除东风港外文中选取的黄渤海验潮站评价结果均为极好,skill均高于0.86,证明模型模拟值与实测值吻合较好。东风港潮位站模拟值与实测值的skill系数偏低,是由于东风港潮位站地处黄河入海口附近,受黄河入海流量流速和地形冲刷淤积变化等原因的影响,使得东风港潮位站的实测值和模拟计算值的吻合度较差,skill系数相对较低。但文中重点研究台风过经曹妃甸区邻近海域,台风到达天津塘沽地区后二次登陆未对东风港邻近海域产生影响,对此仅东风港1站的模拟结果对整体模型的影响较小,文中所建立的数值模型可以准确地对曹妃甸区邻近海域的水动力条件进行模拟计算。
表2 水动力模型效率系数
图3为渤海各典型潮位站天文潮验证结果,黄骅港、塘沽、京唐港和秦皇岛4站为曹妃甸区邻近海域具有重要研究意义的潮位站,图中天文潮表示各潮位站逐小时记录的实测资料进行调和分析后的天文潮实测值,计算值表示由文中所建立的数值模型得出的潮位计算结果。潮位的验证一般将振幅和相位两项作为重点考核指标,其中振幅表示潮位所能达到的最高和最低潮位,相位表示潮波在传播过程中与传播时间的符合程度,从图3中可以看出,计算值与预报值的振幅和相位均吻合良好,表明模型能较好地模拟潮波在黄渤海内的传播。
图3 渤海各潮位站天文潮验证结果图
由于曹妃甸及其邻近海域是强潮海区,大潮、中潮和小潮高潮位之间的差异较大,在台风路径和中心气压已知条件下,风暴潮潮位高低主要取决于台风增水与天文潮潮位遭遇时刻。对此,本文通过设计外海天文潮过程,使曹妃甸站天文高潮位为6-10月高潮位的多年平均值,并调整风场时间,确保风暴潮最大增水与天文潮最高潮位在曹妃甸站遭遇。
台风在曹妃甸邻近海域运动时,低压台风场驱动潮流高速运动,同时中心气压场附近的水位发生明显增加,曹妃甸潮位站有2.754 m 的最大增水,已超过当地红色警戒潮位0.214 m,影响范围和强度相比台风“利奇马”剧烈,曹妃甸潮位站最大增水时刻水位及流场如图4 所示。
图4 曹妃甸潮位站最大增水时刻水位和风生流场图
当曹妃甸潮位站预警到最大增水时,未来2~8 h内台风引起的增水会迅速淹没浅滩和堤防相对较低的地带,随后在地势低洼的区域扩散开来,导致城市积水和内涝现象发生,模拟台风过程主要影响范围集中在曹妃甸东部区域柳赞镇附近,曹妃甸漫滩漫堤淹没过程如图5所示。
图5 台风影响曹妃甸区漫滩漫堤淹没过程图
本文选取历史上对曹妃甸区影响最严重、风暴增水最显著的“7203”台风作为典型台风过程。基于曹妃甸区现存的防潮设施,并按照以上台风参数和外海边界条件,计算台风风暴潮影响下曹妃甸区的淹没水深分布、危险性等级、脆弱性等级和风险等级评估与区划,其中在确定曹妃甸区各级淹没水深分布的面积时,本文根据曹妃甸区受台风实际淹没情况共设置4 级统计区间分类,分别为(0.00 m,0.15 m)、[0.15 m,0.50 m)、[0.50 m,1.20 m)和[1.20 m,3.00 m);在确定危险性等级分布的面积时,本文采用淹没水深作为风暴潮危险性等级划分依据,危险性等级与淹没水深对应关系为:当淹没水深高于3.00 m 时,危险性等级为Ⅰ级,危险性最高;当淹没水深介于1.20 m 和3.00 m区间时,危险性等级为Ⅱ级;当淹没水深介于0.50 m和1.20 m 区间时,危险性等级为Ⅲ级;当淹没水深介于0.15 m和0.50 m区间时,危险性等级为Ⅳ级,危险性最低。在确定风险等级分布的面积时,本文依据曹妃甸区风暴潮灾害危险性和脆弱性分析结果,进行风暴潮灾害风险评估。最后以镇或街道为单元统计曹妃甸区各淹没水深、风险等级分布的面积,结果如图6所示。
图6 曹妃甸区风险评估与区划过程图
风暴潮灾害风险评估计算基于区域灾害系统理论,综合危险性和脆弱性评估结果,评价风暴潮灾害风险。脆弱性评估采用定性评估,风险等级评估公式为:
式中:R(Risk)代表风险;H(Hazard)代表危险性等级分布,V(Vulnerability)代表脆弱性等级分布,主要依据不同二级土地利用类型斑块所占面积比例确定社区(村)脆弱性等级。“×”为风险等级识别矩阵,依据危险性等级与脆弱性等级按表3 综合形成风险等级。为了方便使用,本文依据危险性和脆弱性等级划定标准,采用赋分的方式依次对不同评价等级进行表示,将等级分别赋值为1、2、3、4,其中1代表评价等级Ⅳ级,4代表评价等级Ⅰ级,最后根据风险评估公式可将风险等级R划分为低(Ⅳ级,1≤R≤3)、中(Ⅲ级,3<R≤6)、高(Ⅱ级,6<R≤9)、极高(Ⅰ级,9<R≤16),其中Ⅰ级为高风险,Ⅳ级为低风险。
表3 风险等级R划分标准
本文通过对台风风暴潮影响下研究区域内淹没水深分布、脆弱性等级、危险性等级、风险等级评估与区划结果的分析,可得出曹妃甸区现存的大部分防潮设施可抵御该次风暴潮的影响,仅有纳潮河南岸等少部分区域被淹没。曹妃甸北部和东部区域整体风险较小,主要风险为Ⅲ级和Ⅳ级。Ⅰ级和Ⅱ级高风险区域主要集中在入海河道两侧。
本文根据研究区域风暴潮数值模拟、淹没水深及风险评估与区划的结果,基于研究区域目前风暴潮防灾减灾现状,以减轻灾害风险为目的,提出如下对策建议。
(1)根据曹妃甸区风险评估与区划的结果可知,绝大部分的Ⅰ级和Ⅱ级高风险区域集中在入海河口及河道两侧,建议重新核查曹妃甸区入海河口及河道两侧现有防潮防台风减灾工程和非工程体系,加强海防工程建设,主要内容包括以下4 方面:①核查现有入海河口闸门信息,针对主要入海河口,补充建设闸门防护系统;②核查现有沿海防护堤堤顶高程,核实沿岸海堤是否发生沉降,对已发生沉降的沿岸海堤进行加固;③核查主要入海河道两侧高程,对高程偏低的区域进行防护堤建设;④为安全起见,现有各入海河道闸门应做好防潮防台风应急预案。
(2)根据曹妃甸区相关管理部门调研可知,曹妃甸区现存的应急避灾点较少,建议根据本应急疏散路径的结果对避灾体系进行梳理规划,优化避灾体系,加强对海洋灾害的安全避灾点选址与建设,选址离易灾区域近、地势高、结构坚固、耐水浸泡的场所作为海洋灾害避灾点。
(3)根据风暴潮数值模型和淹没水深分布的结果,建议核验现有警戒潮位标志设施的警戒潮位,构建风暴潮灾害警示告知体系,标明各台风等级下该地可能出现的最高潮位值和最大淹没水深;根据风险评估与区划成果,在曹妃甸区南部沿海区域设置警示杆,并核验现有警戒潮位标志设施的警戒潮位;在避灾点显眼区域挂避灾点铭牌,在主干路口和拐点处设立疏散指示牌。
本文以曹妃甸区为例,采用SMS 软件开展前处理工作,实现了黄渤海大尺度和曹妃甸区精细化非结构网格的绘制,选取近岸海洋数值模型ADCIRC 水动力模型,完成了风暴潮漫滩和漫堤数值模拟,并运用GIS系统及数字可视化技术绘制了曹妃甸区淹没水深分布、危险性等级和风险等级分布图,得到如下结论:
(1)在台风风暴潮的影响下,曹妃甸区绝大部分被淹没区域的形成是因入海河道两侧高程偏低而发生漫滩导致,仅有纳潮河南岸等少部分被淹没区域的形成是由高程偏低而发生漫堤导致的。
(2)曹妃甸北部区域主要以水田和水浇地为主,东部区域主要以沿海滩涂为主,尚未布置工业企业,该区域整体风险较小,主要风险为Ⅲ级和Ⅳ级,Ⅰ级和Ⅱ级高风险区域主要集中在入海河道两侧。