0 引言

茅洲河流域位于珠江三角洲东南部，伶仃洋东岸，茅洲河干流全长41.61 km，其中下游11.40 km河段为深圳市与东莞市的界河，对深圳和东莞两市影响重大^[1]。

茅洲河两侧的深圳市和东莞市作为粤港澳大湾区的重要组成城市，具有重要的政治、经济地位。然而茅洲河两岸未进行过综合整治^{[2, 3]}，当发生影响较大风暴潮时，茅洲河两岸经常被淹，对当地社会经济及人民财产安全造成重大损失。随着《粤港澳大湾区规划纲要》^[4]的公布，茅洲河作为深圳、东莞两市的界河，得到政府高度重视。在此背景下，开展茅洲河风暴潮风险评估分析是进一步深化流域综合管理^[5]，切实落实《粤港澳大湾区规划纲要》^[4]的当务之急。

国内外学者已对风暴潮风险评估开展了大量的研究，部分研究者采用风暴潮数学模型^{[6, 7]}研究风暴潮的控制因子和发生机理，预测某一场风暴潮的强度和范围，评估风暴潮造成的危险，但不能对风暴潮危险性的长期趋势做出预测。目前大部分地区历史台风的信息记录时间序列仅为数十年，对于百年一遇或者千年一遇历史台风的记录很少或者信息缺失严重，这就导致在应用于实际的长期灾害风险评估时会出现样本不足的问题，而重现期风暴潮的研究可以解决研究资料时间序列不足的问题^[8]。受全球气候变暖和海平面上升的影响，不同重现期的风暴潮极值水位将发生很大变化，随着近几年超强台风的频繁发生，台风引起的近海风暴潮最大可能增水是防汛部门和海洋工程设计单位研究的热点，也是近海风暴潮风险评估关注的重点^[9]。开展基于频率分析的典型重现期风暴潮评估及最大可能风暴潮预测，可兼顾短期危险和长期危险，有效识别重点区域未来可能遭遇风暴潮的频率和强度^[10-12]。

本研究基于茅洲河风暴潮实测资料的分析，将影响该水域的特大重现期风暴潮及未来可能发生的最不利风暴潮两种尺度相结合，评估茅洲河遭遇风暴潮的可能性风险，这不仅可以提供风暴潮袭击时的范围和强度，而且可以为当地的城市布局规划和堤防建设提供科学依据。

1 茅洲河口风暴潮特征

通过统计1949-2017年共69 a的国家气象局台风资料，总计有303个热带风暴级以上的热带气旋对茅洲河造成影响。主要影响是造成茅洲河口潮位急剧升高，对茅洲河水位造成顶托。69 a间，平均每年4.39个台风，其中超强台风41个，强台风67个，台风73个。

经统计，影响茅洲河的热带气旋有珠江口正面登陆型热带风暴33次，粤西登陆型热带风暴102次，珠江口外以东登陆型热带风暴111次，海南、广西登陆型热带风暴32次。此外，还有26次未在我国华南海岸登陆。由此可见，影响茅洲河口的台风类型以粤东及粤西登陆型占比较大。而以历史统计资料来看，尤以在珠江口和粤西沿岸登陆的台风成灾最为严重^[13-15]。

2 风暴潮数值模型

风暴潮是一种来自海上的巨大自然灾害，指由于强大的大气扰动，强风和气压骤变所导致的海平面异常升高现象^{[16, 17]}。风暴潮是海洋与大气扰动的大尺度运动而成，在构建数学模型时需要考虑足够大的海气范围才能拟合出较为准确的风暴潮大气扰动及海洋增水过程^[18-22]。珠江河口槽滩相间，河网复杂，岸线曲折多变，水动力条件十分复杂，需要较精细的网格才能高精度地模拟珠江三角洲的地形和岸线^[23]。因此，本研究采用ADCIRC数学模型^[24]，建立大尺度的南海潮汐风暴潮耦合数学模型，在考虑风暴潮作用的同时，考虑洪水、天文潮协同动力机制。

2.1 模型构建

南海潮汐风暴潮模型覆盖南海大部分海域，计算域为12.7°-29.4°N，105.6°-124.5°E，深海处网格尺寸为2′×2′，水深起算基面为平均海平面，模型外边界采用62个天文分潮调和常数控制。计算水深采用美国国家地球物理资料中心(NGDC)数据。

珠江河口上游边界至珠江河口八大口门，茅洲河至上游15 km处，网格最小尺寸为5 m，近岸海域地形采用的是2013-2016年实测地形(图 1)。

2.2 叠加风场

台风风场是产生风暴潮的决定因素，现有的国内外风暴潮模式一般都是采用台风经验模型来模拟台风风场^[25-27]，受台风路径、风速、时间、海域等因素影响，经验风场对台风中心模拟较为准确，但对台风外围的背景风场模拟精度较差^[28-31]。本文采用经验风场与National Centers for Environmental Prediction (NCEP)再分析风场相叠加的方式构造一种叠加风场，在距离台风中心较近处采用经验风场，准确描述台风中心位置、中心风速、最大风速半径处风场特征；在远离台风中心处采用NCEP再分析风场，准确描述风场的背景风速特征，中间采用平缓的过渡，保证两种风场衔接的连续性。从模拟的风场计算结果来看，叠加风场更接近于真实台风风场的变化趋势，风速极值误差为±15%^[29]，该风场已多次应用于南海风暴潮及台风浪预报，且均取得较好的计算成果。

2.3 模型率定与验证

天文潮率定主要采用2008年9月和2017年8月实测潮汐资料。风暴潮验证选取对工程位置造成明显影响的10场台风进行验证，包括6903号超强台风、7908号强热带风暴、9107号台风、9316号台风、9504号强热带风暴、0104号台风“尤特”、0313号超强台风“杜鹃”、0814号强台风“黑格比”、1319号超强台风“天兔”及1713号强台风“天鸽”等。由于文章篇幅限制，仅提供部分工程附近水文站点验证图(图 2)。

在各台风条件下，各主要潮位站模型与原型的潮位过程线吻合良好，模拟高潮位与实测值较为接近，相对误差在13%以内，该风暴潮模型可用于茅洲河风暴潮风险分析(图 3)。

3 茅洲河漫堤风险分析 3.1 重现期风暴潮漫堤风险分析

重现期风暴潮研究主要是预测一个区域在未来10-1 000 a发生风暴潮灾害的可能性，它是一种长时间尺度的趋势预测研究，对于沿海工程的设计及灾害评估意义重大^[32-34]。由于茅洲河附近无长时间序列的水位监测站，本研究采用南海潮汐风暴潮耦合模型对连续20 a对茅洲河影响较大的106场风暴潮进行数值模拟，得到每年茅洲河口风暴潮最高潮位。通过极值Ⅰ型统计学方法对连续20 a风暴潮最高潮位进行频率统计(图 4)，得到茅洲河口百年一遇风暴潮高潮位为2.98 m，五十年一遇风暴潮高潮位为2.78 m(本文所采用的基面均为珠基)。

茅洲河左岸有多处地势低洼处，最低处高程为2.17 m。从表 1可以看出，当发生超十年一遇风暴潮时，即有漫堤风险，发生千年一遇风暴潮时，左岸多处漫堤，漫堤总长度约为1 524 m，随着风暴潮设计频率的降低，漫堤长度逐渐减少。

茅洲河右岸地势较左岸高，高程为3.20-6.58 m。当发生超百年一遇风暴潮时即有漫堤风险，当遭遇千年一遇风暴潮时，右岸漫堤总长度为668.14 m，漫堤比例为5.08%。

3.2 最大可能风暴潮漫堤风险分析

在沿海重大工程项目建设时，应充分考虑风暴潮袭击的频繁程度和最大可能台风增水。深圳和东莞作为粤港澳大湾区的重要组成城市，为保障地区经济安全高效发展，防潮工程中适当考虑可能最大台风增水是必要的^[35]。本研究主要方法为根据历史上在特定海岸地区发生最大持续风速所选择的气象参数值的组合，以这些气象参数值构造可能最大台风，其主要参数包括台风中心气压、最大风速半径、台风移动速度以及台风登陆时与海岸的交角。这些参数和台风可能移动的路径不是虚构的，而是依据统计方法计算后得到的，是可能发生的各种不利情况的组合。

3.2.1 统计法构造可能最大台风

(1) 台风中心气压

ΔP为台风外围不受台风影响处的气压，或正常气压和台风中心气压两者之差，取千年一遇的数值，P_∞取我国沿海夏季海平面气压的平均值(1 008 hPa)。资料系列取1960-2017年共计58 a。对ΔP系列进行PⅢ型频率分析(图 5)，其中千年一遇中心气压为911 hPa，百年一遇中心气压为927 hPa。

(2) 台风的移动速度

台风中心的移动速度V_F是重要的台风参数之一。台风风暴数值模式的研究成果表明，对于开阔水域，登陆台风或紧靠岸边又平行于海岸移动的台风，移动速度越快，风暴增水越大。对于茅洲河附近海域，上承径流影响，外受海水顶托，台风增水情况要复杂得多，一般采用统计结果确定该参数。根据历史资料统计，台风最大移动速度为37 km/h，最小移动速度为10 km/h，平均移动速度为25 km/h。在进行可能最大台风计算时，台风中心移动速度取其平均值25 km/h。

(3) 最大风速半径

台风最大风速半径R是指台风中心到出现最大风速处的半径距离，是衡量台风尺度的重要参数之一，也是风暴潮模式计算时所采用的台风气压场和风场的重要参数之一。

最大风速半径采用经验公式：

$ \begin{array}{*{20}{c}} {R = 28.52{\mathop{\rm th}\nolimits} [0.0873(\phi - 28)] + }\\ {{\rm{12}}{\rm{. }}22\exp \left[ {\frac{{{P_0} - 1013.2}}{{33.86}}} \right] + 0.2{V_F} + 37.22} \end{array} $

式中，V_F为台风中心移动速度，取25 km/h；地理纬度ϕ取22.73°N；P₀为台风中心气压，取911 hPa。计算得最大风速半径R为29 km。

(4) 最大可能台风移动路径

设计台风路径时，需确定台风登陆地点和台风移动方向两个要素。为构造假想台风，搜寻近70 a影响南海的台风，从台风强度和致灾严重程度来看，1713号台风、1208号台风及0814号台风造成舢板洲站超过180 cm的增水。在设计登陆点时，考虑1713号台风的珠海登陆点、1208号台风的台山登陆点及0814号台风的电白登陆点，以及在1960-2017年间，共有13场台风登陆阳江，例如7610号台风、8404号台风、0218号台风及1311号台风等，表明阳江遭受台风侵袭的次数频繁，因此，取1311号台风在阳江的登陆位置作为设计台风登陆点。珠江口西侧的澳门是经常受台风侵扰的区域，其中登陆澳门的6411号台风“艾达”强度大，最低气压925 hPa，登陆时中心气压972 hPa，中心附近最大风速33 m/s，取澳门作为登陆点之一。

研究发现，导致茅洲河口出现大幅增水的台风均为珠江口以西登陆。在台风路径设计时，选取270°-350°、间隔10°共计8个台风行进方向、5个登陆点进行风暴潮增水数值实验，以获取最大可能台风的路径。数值试验方案见图 6，图中红色五角星为5个设计台风登陆点，不同计算方案下，茅洲河口台风增水数值试验结果如表 2所示。

(5) 最大可能台风增水数值模拟

根据以上研究，确定登陆地点为澳门，移动方向为340°是最大可能台风移动路径。设置台风的参数：P₀=911 hPa，R=29 km，V_F=25 km/h，登陆地点为澳门，移动方向为340°，进行风暴潮数值模拟，绘制最大可能台风期间茅洲河口台风增水过程线，结果如图 7所示，最大可能台风增水发生在台风登陆后2-3 h，茅洲河口最大可能增水达到3.64 m。图 7为最大可能台风登陆3 h珠江口内台风增水分布图。在台风推动下海水向珠江口内涌入，台风增水由伶仃洋1.0 m向珠江口内递增，口内台风增水达到3.5 m以上。

3.2.2 最大可能风暴潮情况下茅洲河风险分析

根据最高潮位与洪水遭遇的分析结果，茅洲河口风暴潮与大洪水遭遇的可能性不大^[3]。因此，本次茅洲河沿程潮位计算组合采用最大可能风暴潮与茅洲河上游2 a一遇设计洪水(550 m³/s)遭遇，基本能涵盖茅洲河口最高潮位与上游洪水遭遇情况，是一种合理的遭遇组合。

通过计算，茅洲河口的风暴潮极值潮位为4.78 m，沿程水位较高，现状堤防设计标准较低，茅洲河两岸会发生大面积漫堤(图 8)。茅洲河左岸漫堤总长度为9 830 m，漫堤岸线比为74.68%；右岸漫堤总长度为7 310 m，漫堤岸线比为55.50%。茅洲河左岸最大漫水高度为2.79 m，最小漫水高度为0.04 m，其中利源石油—瑞云科技园段、海滨路南段、润安汽车检测站附近及洋涌路东段漫水高度超过1.00 m。茅洲河右岸最大漫水高度为1.97 m，最小漫水高度为0.02 m，其中丽城科技工业园—澳维力工业园段、基达利工业园段、松源楼工业园—碧头平泰昌科技园段及大洋洲高新产业园—燕罗湿地公园段漫水高度超过1.00 m。

4 结论

通过实测资料分析，在南海潮汐风暴潮数学模型验证准确基础上，对影响茅洲河重现期风暴潮及未来发生的最大可能风暴潮进行灾害计算分析。茅洲河口百年一遇风暴潮高潮位为2.98 m，五十年一遇风暴潮高潮位为2.78 m。茅洲河左岸有多处地势低洼处，当发生超十年一遇风暴潮时，即有漫堤风险，发生千年一遇风暴潮时，茅洲河多处漫堤，漫堤总长度约为1 524 m，随着风暴潮设计频率的降低，漫堤长度逐渐减少。当台风参数为P₀=911 hPa，R=29 km，V_F=25 km/h，登陆地点为澳门，移动方向为340°时，将发生最大可能风暴潮，茅洲河口的风暴潮极值潮位为4.78 m，茅洲河两岸多处漫堤，左岸漫堤总长度为9 830 m，漫堤岸线比为74.68%。

从计算结果可得，重现期潮位的计算方法可以较好地解决无测潮站地区的设计水位推算问题，以及最大可能风暴潮的设计，能最大可能估算特定区域的风暴潮风险情况，并在此基础上制定风暴潮风险防范策略及减灾策略。本文研究成果可为茅洲河综合治理提供重要技术支撑，同时也为海洋及水利防灾减灾部门应急减灾提供决策参考。