舟山绿色石化三维潮流特性数值研究

2022-06-18陆凡沈良朵高郁王晋宝

海洋学报 2022年5期

陆凡，沈良朵*，高郁，王晋宝

(1.浙江海洋大学船舶与海运学院，浙江舟山 316021)

1 引言

在经济全球化快速发展的大背景下，世界各国对海洋资源的开发利用越来越重视。特别地，我国大部分沿海城市在港口码头、物流集散中心和绿色石化基地的建设投入逐步加大[1]。然而，部分以围填海方案来获取工业用地的建设方式由于缺少科学的指导依据，对周边海域的生态环境、水动力条件和泥沙冲淤造成了不良影响[2]。基于此，目前通过采用数值模拟分析法以水工建筑物对周边海洋环境的影响范围为目标进行科学预测，能够合理有效地解决如何绿色开发海洋资源的问题。

围填海工程前能够科学地模拟出建设海域相关的潮流运动形态差异和泥沙冲淤的显著地区，对工程实施有着积极的推动作用。因此，一些国内外专家学者借助数值模拟进行了相关研究。任一晗等[3]基于FVCOM 三维水动力数值模型，选取1984 年、2010 年、2019 年3 个代表年份，探讨围垦工程影响下舟山群岛海域潮流结构与潮能分布的时空变化状况；Safavi等[4]采用MIKE3 平台构建了乌鲁米耶湖的三维数值模型，研究了受湖中堤道分隔影响湖中南北部水中盐度的变化情况；窦明等[5]利用FVCOM 三维海洋数值模型，模拟嵊泗大洋山围垦工程建设填海前后的流场及泥沙运动，对该工程附近海域潮流及各层泥沙分布变化进行了探讨；刘嘉星和刘长根[6]利用验证后的三维潮流数学模型ROMS 探讨了天津港附近海域的潮流特征、潮流的分层特点和围海造地引起的潮流场变化情况；刘金鹏等[7]以龙口市人工岛为例，利用验证合理的潮流场和波浪场数学模型，模拟分析了在极限风况条件下人工岛群周围海域水环境的变化特征及人工岛岸线的波高特征值的变化规律。Zamani 和Koch[8]以伊朗西南部Maroon 水库为研究对象，通过与实测数据进行比较的基础上，对这两种三维水动力模型的适用性进行了对比研究；Soltani等[9]通过Mike21，2011 Version 模型，模拟了海岸带的波浪和水流，以确定奈班德海湾地区的泥沙输移速率和方向；Elshemy 等[10]使用MIKE21 建模系统开发了曼扎拉湖的水动力模型和水质模型，研究了未来气候变化对湖泊水动力和水质特征的影响。

虽然前人对围垦工程的数值模拟研究已经取得不错的结果，但大多数都是以FVCOM 三维模型或者MIKE 二维水动力模型对相关海域的水动力和温盐变化进行了细致研究，而采用MIKE3 和MIKE21 这两种水动力模型对国内近海岸岛屿围垦工程的对比性研究较少。基于此，本文以舟山绿色石化基地三期围填海工程为研究对象，基于MIKE3 平台构建了工程前后三维水动力数学模型，并考虑到MIKE21 中的二维模型也可以分析水平方向上的潮流运动情况，故对两种模型进行了潮位验证结果精度的比较分析。但结合到舟山群岛海域水道纵横、地形结构复杂的特性，水平方向上的二维研究无法呈现出三维模型展现出的流场空间结构变化，因此在三维模型下分析了采用“裁弯取直”和水道封堵方案下对潮流场特征的影响，并从表层、中间层和底层分别对工程前后的流速流向影响进行三维垂向剖析研究。

2 工程概况

本工程位于舟山群岛和岱山以西的大、小鱼山岛海域，北濒岱衢洋、南接灰鳖洋，周边岛屿岸线曲折，港湾众多。水下地形南北深、东西浅，海底高程基本在-65～-1 m，水下地形起伏明显[11]。三期围填海工程地处于121°55′11″～121°59′37″E 之间，其用海面积约为18.41 km2，成陆面积为17.09 km2。工程涉及海域潮流运动以往复流为主，而且基本为平行岸线方向的往复流。涨潮流主要为来自东海的前进波经由舟山群岛之间的狭道进入杭州湾，其中岱山海域的涨潮流主流较为平顺，整体上呈现由东南流向西北方向。落潮流则为杭州湾的落潮水流经舟山群岛之间的狭道流入外海，落潮流基本与涨潮流一致，岱山周边落潮流整体较为平顺。海区内泥沙主要来源于长江口以及杭州湾外泄泥沙随浙东沿岸流的输移扩散，含沙量总体较大。

3 三维水动力模型

3.1 控制方程

三维模型使用的是分层网格，即在水平域中使用非结构网格，而在垂直域中使用结构化网格。垂直网格是分为sigma 坐标或组合的sigma/z 级坐标。本文采用的是考虑了自由曲面的sigma 坐标。使用sigma坐标最重要的优势是它们能够准确地表示水深测量，并在海床附近提供一定的分辨率[12]。采用单元中心的有限体积方法进行了原始方程的空间离散化，方程式如下：

式中，t为时间；x、y、z为笛卡尔坐标；η为表面高度；d为静止水深；h=η+d为总水深；u、v、w为x、y、z方向的速度分量；f=2Ωsinφ，sinφ为科里力参数（Ω为角转速，φ为地理纬度）；g为重力加速度；ρ为水密度；sxx、sxy、syx、syy为辐射应力张量的分量；νt为垂直湍流（或涡流）黏度；pa为大气压力；ρ0为水的参考密度；S为由于点源引起的排放量大小；us和νs为水被排放到环境水中的速度；Fu和Fv分别为x和y方向上的水平应力分项。

3.2 模型建立

3.2.1 计算区域及网格

考虑到工程区域附近涨落潮的流路经杭州湾和部分东海海域，因此所选取区域范围远大于三期工程计算区域，北界至32°N 左右，南界至28.5°N，东到125°E 附近，潮流界顶点江阴作为长江上游边界，澉浦潮位站附近作为杭州湾上游边界，计算区域位于120°～126°E 之间。三维模型采用sigma 网格，垂向均匀分10 层，工程附近海域网格进行加密，网格最小尺度为18 m，工程外海区域，网格距为8 000 m，不同尺度网格之间通过设置实现平滑过渡，共计52 007 个节点，99 195 个单元。对于二维模型，共计41 041 个节点，78 473 个单元。模型计算区域及水深如图1 所示。

图1 模型区域和水深Fig.1 Computational domain and water depth

3.2.2 模型参数设置

模型计算时间步长设置根据克朗数CFL 条件进行动态调整，在确保模型计算稳定进行下，平均时间步长为2 s，最短时间步长为0.01 s。干湿水深判别设置为干水深0.005 m，淹没水深0.05 m，湿水深0.1 m。这有利于保护模型稳定运行。水平涡黏系数经率定取0.28，垂直涡黏系数使用对数公式确定。通过粗糙高度来表征底床糙率，率定后取0.003 m。边界采用水位控制，即用潮位预报的方法得到开边界条件。外海开边界潮位根据由16 个主要分潮（M2、S2、N2、K2、K1、O1、P1、Q1、MU2、NU2、T2、L2、2N2、J1、M1和OO1）的调和常数，以开边界潮位预报公式计算得到：

式中，A0为平均海面；Fi和 (ν0+u)i为天文要素；Hi和 gi为调和常数。在闭边界取流速的法向导数为0，在潮滩区采用漫滩边界处理。

4 模型验证

由于本工程周边缺少实测数据，故本次模型验证采用了距鱼山岛较近的小洋山海域，由中交第三航务工程勘察设计院有限公司在水文测验期间的潮流实测数据。其中，分大中小潮对临时测站T1 和T6 进行潮位验证（大潮测流时间为2015 年5 月4 日至5 日；中潮测流时间为2015 年5 月7 日至8 日；小潮测流时间为2015 年5 月12 日至13 日）；同样地，对定点测站N1、N14 以及走航测站WP05 流速流向数据进行潮流分层验证，各测站位置见图2。

图2 验证点位置Fig.2 Verification point location

采用MIKE3 的水动力模型计算数据与实测数据进行验证，潮位验证结果如图3 所示。图中显示，从整体上看T1 测站和T6 测站潮位拟合较好，其中小潮阶段拟合效果最好，相对误差结果在±10%左右，个别高潮位处相对误差稍大。总体来看，三维模型模拟的潮位验证是符合要求的。

图3 潮位验证Fig.3 Tide level verification

对潮流测站N1、N14 和WP05 分为表层（第10层）、中层（第5 层）和底层（第1 层）进行流速流向验证。在分层验证下，可以清楚看到表层小潮阶段的流速拟合结果是相对较大的，其他各层拟合结果都是极好的。表层模拟结果如图4 所示，各测站仅在高潮位由于流量加大的影响导致流速拟合结果有0.5 m/s左右的误差，其他时间段拟合相对误差控制在8%左右，表层流向模拟结果令人满意。中间层模拟结果如图5 所示，各测站除高潮位处流速存在0.15～0.3 m/s 的误差，其余各时间点相对误差在10%以内，流向拟合结果较满意。底层模拟结果如图6 所示，大潮和中潮阶段的模拟流速仅在N1 测站的高潮位时有0.35 m/s 的误差，其余拟合较好，小潮阶段除在个别时间点有0.45 m/s 的误差外，其余各点存在0.15～0.25 m/s 的误差范围；底层流向模拟结果较好，流向相对误差在7%以内，极个别点处由于测量过程存在失误而出现的误差，不予考虑。

图4 表层流速流向验证Fig.4 Verification of flow velocity and direction in the surface layer

图5 中间层流速流向验证Fig.5 Verification of flow velocity and direction in the middle layer

图6 底层流速流向验证Fig.6 Verification of flow velocity and direction in the bottom layer

为对模型验证结果有更加科学性的认识，此处通过Willmott[13]所提出的skill 方程进行评估，表达式如下：

式中，Pi是模拟值；Oi是实测值；是实测平均值。通过d值判断模型效率，d＞0.65 时，模型效率极好；0.5＜d≤ 0.65，表示模型效率非常好；0.2＜d≤ 0.5，表示模型效率好；d≤ 0.2，表示模型效率差。

结果表明，N1 和WP05 测站各层的大中小潮期间流速d值在0.80～0.97 之间，流向d值区间为0.73～0.99，其模型效率都是极好的。限于篇幅，仅具体列出了N14 测点各层的模型效率值，具体见表1。

表1 模型效率系数Table 1 The index of agreement

5 三维结果分析

5.1 MIKE3 与MIKE21 对比分析

对MIKE3 与MIKE21 两种拟合效果进行对比分析，分别将这两种水动力模型的计算数据与实测数据进行验证，潮位对比验证结果如图7 所示，由图可知：两种模型的潮位计算值与实测值的吻合度较高，特别是在小潮期间，其最大误差仅有0.15 m 左右。大潮和中潮期间，三维数值结果与实测值在个别高潮处有0.5 m 左右的误差，总体上看，在将二维模型中的初始条件同样地设置于三维模型的情况下，二维数值结果更接近于实测值。

图7 潮位对比Fig.7 Tide level comparison

考虑到理论方面上MIKE3 中间层的潮流模拟结果与MIKE21 的平均流速流向的结果相似，故在此处进行了流速流向的对比验证。中间层模拟结果对比如图8 所示，图中各测站流向验证结果良好。本文选择的水深0.6H的实测数据与垂向平均流速基本保持一致，在大潮和中潮阶段两种模型拟合结果较好，最大相对误差在20%左右，相较而言，小潮阶段的高潮位处误差稍大，相对误差在25%左右，其余各时间点两者模拟结果均良好。通过这4 组数据的对比发现：二维模型的大部分计算值更加接近于实测值，尤其是在涨急时刻，原因可能是二维模型是取值于垂向水深平均的流速进行模拟，而本文三维模型未考虑风场应力的影响，表层流速计算误差会对三维流速计算结果有影响。总体上看，三维模型中间层的模拟结果与二维计算结果是最接近的，且中间层拟合效果最为出色，模拟精度基本满足《海岸与河口潮流泥沙数值模拟技规程》的要求，模型可用于工程后的预测等各项工作[14-15]。

5.2 工程前后潮流场分析

从计算的潮流场分析得出，涨急时刻外海潮波沿东南向西北方向前进。在大鱼山的东南侧受岛屿阻挡，潮流分成两股，分别沿着大鱼山岛东南侧向南北流入杭州湾。大鱼山岛区域涨潮时刻受峡道效应，南岸流速大于北岸，落潮时刻则北岸较大，工程南北两端的矶头岸线对涨落潮流起到了挑流的作用，而且两侧都为深槽，因此涨落潮流速相对较大，但整体还较为平整。其中，涨潮时平均流向为230°～300°，落潮时平均流向为70°～110°，涨落潮平均流向约为180°，基本符合往复流性质。

考虑到围垦工程所处的大小鱼山岛周边范围内岛屿错落，水道纵横，潮流运动情况复杂，故对工程前后计算域内大潮涨落时刻的表层流场变化进行具体的对比分析，同时给出了工程后流速的等值线图以更直观看出流速大小变化，如图9 所示。整体上看，此围填海工程前后对工程周边海域的潮流运动没有较大影响，是明显的往复流形式。在大潮涨急阶段，由于受围垦工程“截弯取直”方法的影响，处于工程东北侧的岸线外伸，东北侧岸线存在一个向内的45°转角，受岸线角度和洋流反射作用，较于工程前平行于岸线的潮流运动方向，此处产生一个直径约为700 m 的逆时针旋涡；鱼山岛的西北侧和东南侧受围填影响，潮流运动遇到外扩海岸线，此处潮流运动方向变化呈现出30°转角；大范围流场前后变化不大。大潮落急阶段，工程西北角处的一处小岛屿围填时与主岛进行了连接处理，将此处原先经狭道效应的高速流与经岸线折射过来的洋流形成的旋转流影响削弱，故旋转流转变成了与岸线弯口一致的运动形式；鱼山岛东北侧、南侧以及西侧落潮流受外延海岸线的挑流作用，工程后流场运动均出现30°～45°的转角变化，这表明处于迎流面的岸线变化对流场的改变起主导作用。

图9 表层潮流场对比（上：涨急下：落急）Fig.9 The comparison of tidal current field on the surface layer (up:rising;down:falling)

5.3 垂直流场分析

为更立体的表现出工程区域的流速流向情况，本节同时选取了N1、N14 和WP05 这3 个典型的测站数据进行垂直结构的分析。由图10 分析得：N1 测站的流速在垂直方向上的实测值和涨落急阶段的计算值各自都较接近，并在中间层附近的水深0.6H处仅有4.6%左右的相对误差，计算与实测的最大误差出现在表层附近，相对误差在-25%左右，涨落急各分层流速较为一致，整体上都呈现出流速随水深增加而减小，且幅度较大，可能是受到底部海床的摩阻力较大的影响；对于各分层流向趋于一致，且与实测数据拟合较好。N14 测站处涨落急时刻的流速计算值与实测值最大误差出现在表层，其相对误差为-28.01%，同时此测站的涨急流速计算值大于落急时刻计算值，这可能是由于此处网格布置较疏的影响，加上缺少风场实测资料，未考虑风场的原因所导致。WP05 测站可能处在洋山港与壳子山之间的狭窄水道里，由宽阔海域进入该水道，过水断面突然减小，故同样是在表层处的流速与实测数据有-30%的相对误差，而0.6H以下水深的计算流速与实测值的误差范围为13～22 cm/s，且涨急时刻流速小于落急流速，同样地各分层的流向是近乎一致的。通过对垂直流场的分析，可以清楚地得到流速与水深呈负相关，流向受水深影响近乎不计。而针对表层拟合较大问题，初步考虑是受到了风场和波浪对潮流的驱动作用。王世澎等[16]在二维情况下波浪对潮流场的数值分析中指出，对于流速过程，波浪在没有破碎的情况下对流速过程的影响比对潮位过程的影响大。考虑到围填海工程后，周围海域受局部挑流作用，加之波浪辐射应力影响，波浪会在近岸破碎产生沿岸流叠加到原有潮流场上，对表层流速影响较大。针对此问题，将会在以后的研究中通过波浪耦合进行率定工作。