潮汐是海洋中最普遍和常见的运动过程, 对海洋内波、海洋环流、海洋生态系统等有不可忽视的影响^[1]。而东海作为边缘海, 拥有广阔的大陆架^[2], 且沿岸地区大多经济发达、人口密度大, 因此深入了解东海潮汐信息对科学、经济、军事等具有重要意义^[3]。20世纪30年代Ogura^[4]利用多个验潮站的实测数据绘制了渤、黄、东海的潮汐分布图; 20世纪50年代末, 我国为124°E西侧海域绘制了更加精细的潮汐分布图; 70年代初, 我国对124°E以东海域的潮汐分析计算^[5]。多年前我国就开始了对渤、黄、东海潮汐的研究探索。20世纪80年代以来, 伴随计算机技术的进步, 数值模拟成为了研究潮汐的重要手段^[6]。Fang^[7]利用精细的数值模拟数据和实测数据, 制作了精确度和完整度都较高的潮汐潮流图。在二维模拟的基础上, 万振文等人^[8]使用三维POM模式首次得到了渤、黄、东海的潮流水平和垂直涡动系数; 张衡等人^[9]利用球面坐标系下ECOM数值模式, 对渤、黄、东海的8个主要分潮进行数值模拟; 朱学明等^[10-11]利用FVCOM海洋数值模式模拟了渤、黄、东海以及日本海、鄂霍次克海、南海北部和部分西北太平洋海域的潮汐和潮流; 罗丹^[12]用POM三维海流数值模式对渤、黄、东海潮波进行模拟, 得到M₂和K₁分潮在该海域的特征; 宋泽坤等^[13]通过中国沿海的14个潮位站的调和常数验证了在渤、黄、东海及琉球群岛构建数值模型, 得到了渤、黄、东海的8个主要分潮M₂、S₂、K₁、O₁、N₂、K₂、P₁、Q₁的分布和传播特征; 袁帅等^[14]利用FVCOM模型和DELFT模型对渤海进行数值模拟, 并证实了FVCOM模式进行局部加密的可行性; 袁帅等^[15]利用FVCOM海洋模式对渤、黄海的八大分潮进行数值模拟, 得到了该海域的潮汐和潮流情况; 董宇航等^[16]分析了舟山群岛附近的潮汐和潮流情况, 并且具体估量在利用潮流能方面的应用; 涂成东^[17]利用SCHISM海洋模式对全球大洋潮波进行数值模拟, 得到全球大洋的潮汐和潮流情况。

基于前人对潮汐的探索, 考虑到海洋数值模型对复杂岸线的拟合情况, 本文使用两个验潮站在2016年的实测海面高度数据对一个FVCOM模式的结果进行验证, 并利用该结果得到我国东海区域的同潮图, 展现了该海域的潮汐特征, 为该海域海洋环境保护和相关工程提供依据。

1 资料及方法 1.1 模式介绍

FVCOM模式是由Chen等^[18]开发的基于自由海面、非结构化三角形网格, 采用有限体积法求解三维原始方程的海洋数值模式。该模式在水平方向采用非结构化正交网格, 垂向采用σ坐标系, 对复杂岸线和海底地形有很好的拟合效果, 被广泛应用于近岸海域。

1.2 模式设置

模式计算区域为26°~37°N, 119°~127°E, 水深数据来源于海图和美国地球物理中心发布的ETOPO2数据, 垂向上分为40层, 模式计算区域共有70 761个节点(图 1), 网格分辨率在河口和近岸区域较高, 在近岸小区域的分辨率约500 m, 在距离岸线较远的大区域分辨率在1~10 km, 内模和外模时间步长为10 s和1 s。模式考虑风场、径流的影响, 并且在开边界考虑温盐驱动(由HYCOM模式得到的海洋温盐再分析数据)和潮位驱动(由全球潮汐数值模式Naotide得到)。模式采用冷启动, 设置水位和流速初始值为0, 计算了2016年4月19日至2017年12月31日期间研究区域的潮位情况, 计算过程中每小时输出一次结果, 选取2017年全年的结果进行验证。

1.3 使用资料

本文使用东海大桥南(30.65875°N, 121.9879°E)和洋山(30.6304°N, 122.0194°E)两个验潮站(图 2)2016年1月1日0时至12月31日23时的逐时潮位数据。

模式结果为FVCOM海洋数值模式在计算区域(26°~37°N, 119°~126°E)得到的海面高度数据, 数据时间分辨率为1 h。选取模式计算的北京时间2017年1月1日08时至12月31日08时的结果进行分析。为达到在保证数值模式计算网格节点分布较均匀的同时降低计算量的目的, 本文的验证中在计算区域构建以0.1°为分辨率的网格, 考虑到研究区域内河口可能导致周围陆地出现数值, 本文仅保留与构建网格节点最邻近且距离小于10 km的数值模式计算网格节点, 处理后的网格节点分布如图 3所示。

1.4 评价模式潮汐模拟准确度的方法

Munk等^[19]利用潮汐的振幅和位相(格林威治时间, 下同)表示与潮汐相关的海表面高度, Wang等^[20]在此基础上, 忽略节点校正和潮汐频率的Doodson数, 得到

$ {\text{ζ}}\left( {j,t} \right) = \mathop \sum \limits_{k = 1}^K {C_k}\left( j \right){\text{cos}}\left[ {{\omega _k}t + {\theta _k}\left( j \right)} \right] , $

(1)

式中K为j位置处总的分潮个数, $ {\text{ζ }}\left( {j,t} \right) $表示在j位置处t时刻的海表面高度, $ {C_k}\left( j \right) $和$ {\theta _k}\left( j \right) $分别表示j位置处分潮k的振幅和位相。

Wang等^[20]在Munk等^[19]以及Fok等^[21]评价潮汐模式和准确性方法的基础上进一步简化, 得到在相对较长的时间T内, 海表面高度的均方误差与K个潮汐分量的振幅和位相有关, 并定义实测数据与模式结果海表高度的均方误差D

$ D = {\left\{ {\frac{1}{T}\mathop \smallint \nolimits_0^T {{\left[ \begin{gathered} \mathop \sum \limits_{k = 1}^K C_k^{\text{o}}\cos \left( {{\omega _k}t + \theta _k^{\text{o}}} \right) - \hfill \\ \mathop \sum \limits_{k = 1}^K C_k^{\text{m}}\cos \left( {{\omega _k}t + \theta _k^{\text{m}}} \right) \hfill \\ \end{gathered} \right]}^2}{\text{d}}t} \right\}^{1/2}}, $

(2)

式中ω表示潮汐频率, $C_k^o$和$\theta _k^o$表示观测得到的潮汐振幅和位相, $C_k^m$和$\theta _k^m$表示利用模式结果计算得到的潮汐振幅和位相。

如果忽略不同频率的潮汐分量之间的交叉项, 当时间T远大于潮汐的周期时有

$ D = {\left\{ {\frac{1}{2}\mathop \sum \limits_{k = 1}^K \left[ \begin{gathered} {\left( {C_k^{\text{o}}} \right)^2} + {\left( {C_k^{\text{m}}} \right)^2} - \hfill \\ 2C_k^{\text{o}}C_k^{\text{m}}{\text{cos}}\left( {\theta _k^{\text{o}} - \theta _k^{\text{m}}} \right) \hfill \\ \end{gathered} \right]} \right\}^{1/2}} , $

(3)

当在相对较长的时间T内, K个潮汐分量引起的海表高度的变化幅度V可以表示为

$ V = {\left\{ {\frac{1}{T}\mathop \smallint \nolimits_0^T {{\left[ {\mathop \sum \limits_{k = 1}^K C_k^{\text{o}}\cos \left( {{\omega _k}t + \theta _k^{\text{o}}} \right)} \right]}^2}{\text{d}}t} \right\}^{1/2}}, $

(4)

由(1)(2)得

$ V = {\left( {\frac{1}{2}\mathop \sum \limits_{k = 1}^K {{\left( {C_K^{\text{O}}} \right)}^2}} \right)^{1/2}}. $

(5)

实测数据与模式结果海表高度的均方误差与潮汐产生的海表高度变化的比值D/V可以用来表示潮汐的相对误差。

2 结果 2.1 潮汐调和分析

潮汐调和分析是指将潮位资料分解成多个分潮, 求出每个分潮的平均振幅和位相的值。Pawlowicz等^[22]给出了一个包括节点校正、推断和分析指定分潮的T_TIDE程序包, 它可用于对1 a或者更短的时间段进行经典调和分析, 采用节点的修正解释一些未定潮汐成分, 以及计算所分析成分的置信区间。

潮汐是海水受天体引潮力的作用而产生的周期性涨落现象^[23], 东海海平面上升十分缓慢, 平均速率约为3.02 mm/a^[24], 因此可以假设相同位置的潮汐在每年的相同时间段内的振幅和位相相同, 在本文中利用T_TIDE程序包对数值模式结果和实测数据进行分析处理, 并对这一过程进行验证。

2.1.1 模式结果进行调和分析过程的验证

用距离东海大桥南最近的模式点数据作为该站点的模式结果, 对该模式结果进行调和分析, 得到东海大桥南在2017年海面高度的平均值为0.59 m, 潮差约为4 m。图 4为东海大桥南2017年模式结果与调和分析结果对比图, 图中蓝线表示的是模式结果中该点处每小时的海面高度, 红线表示调和分析回报的潮位与平均海平面的和(以下简称回报海面高度)。从图中可以看出, 模式结果经过调和分析后的回报海面高度和模式结果中海面高度的变化趋势基本一致; 从数值看, 在这一时段内的高潮时刻, 该点处模式结果的海面高度比回报海面高度高0.5 m左右, 在这一时段内的低潮时刻, 该点处模式结果的海面高度比回报海面高度低0.5 m左右, 而在除了高潮和低潮之外的其他绝大多数时刻模式结果的海面高度和回报海面高度基本一致。

同样对洋山站处的模式数据点进行调和分析, 得到该站点在2017年海面高度的平均值为0.575 m, 潮差约为4 m。图 5为洋山站2017年模式结果与调和分析结果对比图, 图中蓝线表示的是模式结果中该点处每小时的海面高度, 红线表示调和分析回报的潮位与平均海平面的和。由图 5可知, 该位置模式结果的海面高度与回报海面高度的关系和东海大桥南处的关系基本相同, 模式结果的海面高度和回报海面高度变化趋势相同, 绝大多数时刻的模式结果和回报海面高度相同, 与东海大桥南相比在部分高潮和低潮位置处的海面高度差有很小的变化。因此, 在东海大桥南和洋山站处使用T_TIDE程序包对模式结果进行调和分析是可行的。

2.1.2 实测数据进行调和分析过程的验证

通过对实测资料进行调和分析, 得到东海大桥南在2016年海面高度的平均值为2.96 m, 潮差约为4 m。图 6为东海大桥南2016年实测资料与调和分析结果对比图, 图中蓝线表示的是实测数据中该点处每小时的海面高度, 红线表示调和分析的回报海面高度。从图中可以看出, 实测海面高度的变化曲线与调和分析得到的回报海面高度变化曲线高度重合, 除极少数的低潮时刻实测海面高度低于回报的海面高度, 剩下的绝大多数时刻回报和实测资料的海面高度都基本相同。

对洋山站的海面高度进行调和分析, 得到洋山站在2016年海面高度的平均值为2.62 m, 潮差约为4 m。图 7为洋山站2016年实测资料与调和分析结果对比图, 图中蓝线表示的是实测数据中该点处每小时的海面高度, 红线表示调和分析的回报海面高度。由图可得, 在绝大多数时刻实测海面高度的变化曲线与调和分析后得到的回报的海面高度变化曲线重合, 个别日期实测与回报海面高度的极小值有较大差距, 这可能是特殊天气或环境因素导致实测出现的特殊情况引起的, 但总体看, 实测与回报海面高度差值的均值仅为0.003 3 m, 方差仅为0.12, 在东海大桥南和洋山站处使用T_TIDE程序包对实测数据进行调和分析是可信的。

2.2 模式与实测对比

本文利用东海大桥南和洋山两个站点的模式结果和实测数据, 在进行潮汐调和分析时主要考虑渤、黄、东海的4个主要分潮M₂、S₂、K₁、O₁和4个次要分潮N₂、K₂、P₁、Q₁, 计算这8个分潮的潮汐调和常数^[14]。计算结果如表 1和表 2所示, 表中信噪比为振幅与振幅误差比值的平方, 一般认为信噪比大于2的分潮是显著的。

从表 1可以发现, 在东海大桥南处, 将实测资料分解为8个主要的分潮后, 最显著的分潮为M₂、S₂、K₁; 将模式结果分解为8个主要的分潮后, 最显著的分潮为K₁、M₂。从分潮振幅来看, 实测数据和模式结果分解计算后各分潮振幅差值最小的为Q₁, 仅相差0.11 cm, 相差最大的为M₂, 相差23.72 cm, M₂、S₂分潮的振幅差较大, 其他分潮的振幅差较小。从位相来看, 实测数据和模式结果分解计算后各位相差最小的为M₂, 仅相差8.79°, 相差最大的为Q₁, 相差39.88°, 其他分潮中位相差在15°~30°。

从表 2可以发现, 将洋山站的实测资料分解为8个主要的分潮后, 最显著的分潮为M₂、S₂、K₁; 将模式结果分解为8个主要的分潮后, 最显著的分潮为K₁、M₂、O₁。从分潮振幅来看, 实测数据和模式结果分解计算后各分潮振幅差值最大的为M₂, 为31.17 cm, 最小的为Q1, 仅相差0.63 cm。M₂、S₂分潮振幅差为20~30 cm, 其他分潮中分潮振幅差为0~10 cm。从位相来看, 实测数据和模式结果分解计算后各分潮的位相差最小的为M₂, 仅相差2.47°, 相差最大的为Q₁, 相差26.64°, 除N₂、O₁、K₂分潮的位相差约20°~30°, 其他分潮的位相差为2°~12°。

利用1.3中的式(3)和式(5), 对东海大桥南和洋山两个站点, 计算模式和观测海表高度的均方误差D和K个潮汐分量引起的海表高度的变化幅度V。结果如表 3所示。

在两站实测数据与模式结果计算所得的各分潮位相差最大的三个分潮为K₁、Q₂、N₂。与洋山站相比, 东海大桥南实测数据与模式结果计算各分潮位相的差值略大, 但两站的绝大多数实测数据与模式结果计算分潮位相相差都较小。由表 3可知, 两站模式和观测海表高度的均方误差与潮汐产生的海表高度变化的比值较大, 潮汐引起的相对误差分别为25.99%和25.29%, 猜想是由三部分原因造成的: 一是T_TIDE程序包所使用的经典潮汐调和分析方法在非平衡潮上表现较差^[23], 两站所处位置可能会受长江水输入影响, 其中东海大桥南的位置距离长江入海水体更近, 使得相应的误差也稍大; 二是数值模式计算结果受到网格加密、地形平滑等因素影响与真实值可能存在差距; 三是使用距离站点最邻近的数值模式计算网格点数据作为站点的模式数据, 虽然二者邻, 但仍会产生一定误差。

2.3 半日分潮和全日分潮的潮汐特征

通过东海大桥南和洋山两个站实测和模式计算结果的对比, 发现在这两个站点的占主导的分潮为M₂、S₂、K₁, 从赵保仁等^[5]的研究中可知在本文选择的27°N~36°N, 120°E~126°E的海上区域最主要的半日分潮为M₂型分潮, 最主要的全日分潮为M₁型分潮(K₁和O₁的平均值)。由于实测站点仅有两个, 本文将利用模式结果对所研究区域的半日分潮M₂和全日分潮K₁进行讨论。

图 8给出模式结果计算所得M₂分潮同潮图, 与Fang^[7]文中计算所得M₂分潮同潮图对比, 在本文所研究的区域存在一个无潮点, 且该无潮点位于连云港东侧海域, 与Fang计算所得结果相符合。从振幅来看, 在无潮点附近, 图 8和Fang的等振幅线分布和对应位置等振幅线的值基本一致, 等振幅线的值自无潮点向四周逐渐增加; 自33°N以南部分, 两图等振幅线的趋势一致, 等振幅线的值随距陆地距离减小逐渐增加, 在杭州湾附近存在最大值, 但在对应位置图 8的等振幅线的值比Fang的结果小20~30 cm。从同潮时线来看, 图 8中同潮时线的分布、变化趋势和对应位置的值与Fang的结果基本相同。总体来看, 本文绘制的M₂分潮同潮图与Fang的结果在趋势上基本相同, 但在33°N以南区域M₂分潮的振幅比Fang的结果偏小。

图 9为模式结果计算所得K₁分潮同潮图, 与Fang^[7]文中计算所得K₁分潮同潮图对比, 二者都在(123°E, 34°N)附近存在一个无潮点, 等振幅线对应的振幅自无潮点向四周逐渐增加, 除32°N以南对应位置等振幅线的值比Fang的结果小5 cm左右外, 其他位置的等振幅线的分布、变化趋势和对应位置的振幅一致。从同潮时线来看, 图 9与Fang的结果在同潮时线分布上基本一致, 对应位置的同潮时线对应值相同。总体来看, 该模式结果绘制的K₁分潮同潮图与Fang所绘制的基本一致, 在K₁分潮的模拟上较好。

3 结论

本文利用东海大桥南和洋山站两个验潮站的实测海面高度数据对FVCOM数值模式在东海区域的海面高度结果进行验证。通过计算研究区域的8个主要分潮: M₂、S₂、K₁、O₁、N₂、K₂、P₁、Q₁的潮汐调和常数, 得到两站点数值模式计算潮汐的相对误差分别为25.99%和25.29%。使用模式计算结果所绘同潮图与已有研究结果在无潮点的位置、同潮时线、等振幅线的分布和趋势上十分接近。不同之处在于, 与前人研究对比, 本文所绘M₂分潮同潮图中对应位置等振幅线的值低20~30 cm, K₁分潮同潮图中对应位置等振幅线的值低5 cm左右。

与前人直接使用潮汐表或直接使用验潮站对应时间各分潮的调和常数资料不同, 在本文的研究中利用与数值模式计算结果有着时段长度相同但时间范围不同的实测数据进行验证, 这一思路削弱了对于实测数据的限制, 且最终结果能够对于模式精度进行评价。但受经典潮汐调和分析在非平衡潮上的限制、数值模式设置引起的模式计算结果与真实值的差异以及使用最临近数值模式计算网络节点数据作为站点的模式数据产生的误差等因素的影响, 本文得到的结果可能存在不确定性。后续考虑使用更多的实测站点数据, 比较相同位置相邻年份或不同年份潮汐调和常数的差异, 以及使用Pan等^[25]开发的能够适用于所有潮汐的S_TIDE进行进一步的比较分析。

致谢: 论文写作得到王晓春教授的建议和指导, 在此致谢。