舟山复杂海岛条件下某交通码头溢油环境风险评估

李献丽，陆凡，王臣，沈良朵

(1.浙江国际海运职业技术学院海洋装备工程学院,浙江舟山 316021;2.浙江海洋大学船舶与海运学院,浙江舟山 316022;3.舟山市港航工程规划设计院有限公司,浙江舟山 316000;4.浙江海洋大学海洋工程装备学院，浙江舟山 316022)

舟山群岛位于浙江省东北部，岛屿众多，航道纵横，水上交通较为发达。多年来，为满足城市发展需要，港口码头不断建设发展。岛内航线密集而繁忙，船舶航行碰撞发生溢油的因素加大，可能产生的海洋环境危害也愈发得到重视[1]。

港口码头海域环境危害主要源于海上溢油事故的发生，受风向和潮流的影响其危害范围较广。如何采取有效手段预测溢油风险，做到未雨绸缪，就显得尤为重要。海上石油发生泄漏后，油膜受海水密度、温度、盐度等因素的影响发生多种物理、化学和生物变化过程，如扩展、分散、蒸发、乳化、沉降及生物降解[2-3]。基于此，国内外学者开展了一系列溢油模拟分析。20 世纪60 年代，FAY[4]提出了经典的三阶理论对油膜扩散进行解释，将溢油分为3 个阶段；JOHANSEN,et al[5]在油粒子模型中加入Lagrange 追踪法，使得在风浪等作用下油膜破碎分离现象更清晰；CHO,et al[6]同时考虑潮流和潮汐，建立了溢油分散模型，对河北精神号超大型油轮在黄海发生的漏油事故进行计算扩散分布，为防治灾害提供直观的数据材料；ELIZARYEV,et al[7]使用GNOME 模型和HAZAT 轨迹模型2 种不同的模型对Exxon Valdez 溢油轨迹进行了分析，发现影响油在水体中分散的物理化学特性变化的主要因素是风速风向、水温和波高；游立新等[8]通过建立二维水动力模型并耦合了油粒子模型预测了枝城区铁水联运码头前沿发生溢油事故后油膜的漂移扩散过程；吴琼等[9]采用油粒子模型以湛江港航道疏浚工程为例，对其施工期和运营期可能对海域水环境的影响和存在风险进行预测及分析，为大型航道工程建设提供参考。目前，溢油事故预测研究最受青睐的就是“油粒子”模型[5]，其认为油膜是以大量独立的油粒子组成，每个油粒子含有一定油量。互不干扰的特性使得“油粒子”模型可以确切地表现出较厚的油向油膜边缘扩展的过程以及油膜形状在风向上明显拉长的现象，同时相较于传统模式，“油粒子”模型在油膜断裂和迎风压缩等方面也更具合理性。

因此，本文建立了工程区域二维水动力模型，在模型适用性得到验证的基础上，耦合油粒子模型对舟山定海某交通码头工程进行溢油的漂移扩散、扫海面积以及到达敏感区域的时间进行预测分析，为溢油事故应急方案制定提供理论指导[9-12]。

1 二维水动力模型建立

1.1 控制方程

式中：t为时间，s；d为静水深，h；η 为水面高程，m；u、v分别为x、y方向上的分量，m·s-1；pa为当地大气压；ρ为水的密度，kg·m-3；ρ0为参考水密度，kg·m-3；f=2Ωsinφ 指柯氏力参数，s-1(φ 为所处地理纬度，°；Ω 为自转角速率，取0.729×10-4rad·s-1)；fv¯和fu¯为地球自转引起的加速度，m·s-2；sxx、sxy、syx、syy为波辐射应力分量，N·m-2；Txx、Txy、Tyx、Tyy为水平粘滞应力项；S为源汇项；(ps,qs)源汇项水流流速，m·s-1。

1.2 计算区域及网格

本文数值模型计算区域设置范围较大，为28.5°N——32°N，120°E——126°E。其中长江水域上游边界设在江阴处，杭州湾上游边界设在澉浦潮位站。以非结构三角网格划分，在远海区域设置最大网格尺度5 000 m，工程区则进行网格局部加密，最小网格距20 m，两者之间采用模型平滑迭代过度，便于更好地描绘出岸线和水下地形，减少模型计算时间。模型共计41 041 个节点，78 473 个单元，水深地形及网格分别如图1、2 所示。

图1 计算域网格和水深示意图Fig.1 Grid and water depth diagram

图2 局部网格示意图Fig.2 Local grid diagram

1.3 模型验证

模型验证选择中交第三航务工程勘察设计院有限公司水文测验期间的实测资料，对定海和镇海2 个长期潮位站进行潮位验证，对定点测站N1 和N14 分大小潮进行流速流向验证(大潮：2015.05.04——2015.05.05；小潮：2015.05.12——2015.05.13)。测站位置如表1 所示。

表1 验证点位置Tab.1 Verification point location

1.3.1 潮位验证

图3 潮位验证结果显示，整体数值结果趋势与实测值拟合较好，误差仅在6 cm 以内，相对误差控制在5%左右，仅在部分最低潮处误差存在15 cm 左右，但总体吻合度较高。

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

1.3.2 潮流验证

图4 潮流验证结果显示，大潮阶段计算结果与实测值基本保持同步，相比之下，小潮阶段2 个测站个别最高潮处有所跳动，但整体误差基本在10%左右。流向结果拟合度较高，相对误差在7%以内，并能看出潮流运动呈现往复流特征。

图4 潮流验证Fig.4 Power flow verification

综上所述，模型精度能够基本上满足《海岸与河口潮流泥沙模拟技术规程》[13]的要求，能够反映出涨落潮时刻变化以及潮流运动特征，可以为溢油模型提供较好的流场驱动力。

1.4 潮流场分析

图5 是工程区附近海域大潮涨落急时刻的流场图。由图可知，涨急时刻潮流主要是来自东海的前进波，经舟山东南外海岛屿间水道而分流抵达工程区域。工程区东西侧海域均较为开阔，地形起伏变化不大，故整体上涨落时刻的流路呈现出往复流的形式，且基本是平行于岸线方向。而流速在岛屿狭道处受到束水作用，过水断面减小，流速最大可达1.2 m·s-1；对于落急时刻，仅有方向上反向的差异。

2 溢油模型预测

2.1 计算方程

(1)在风和流的共同作用下，油粒子群的每一个油粒子的运动可用下式表示：

式中：X、Y为经过后油粒子的位置；X0、Y0为某质点的初始坐标；U、V分别为X、Y方向的流速分量，m·s-1，包括潮流和风海流2 部分，流场由前述潮流模式计算得到；W10为海面上的风速，m·s-1；A为风向，°；α 为风拖曳系数；r为随机走动距离(扩散项)，是由水流的随机性脉动所导致每个油粒子的空间位移，r=RE，R为0～1 之间的随机数，E为扩散系数；B为随机扩散方向，B=2πR。

(2)风拖曳系数选取

风海流采用如下计算公式：

式中：α 为风拖曳系数；f(θ)为科氏力引起的偏转角的函数；θ 为偏转角，°。在北半球，风海流向右偏转于风向，此处偏转角取15°，α 取0.03。

风拖曳系数是海洋大气物理学中的重要参数，本文中采用JIN Wu[14]的公式，即：

式中：Ca=1.255×10-3；Cb=2.425×10-3；Wa=7 m·s-1；Wb=25 m·s-1。

2.2 预测方案

本次研究对象的溢油风险事故主要来源于码头装卸作业产生的溢油和来港船只航道上受海损事故产生的溢油。根据调查[15]发现，操作性事故造成溢油量相对较小，而海损性事故虽然发生概率较小，但往往溢油量较大，对环境的危害较为严重。此码头不停泊油船，故以操作性事故作为码头风险事故源项，一般性海损事故作为航道溢油事故源项。

2.2.1 研究区域

舟山群岛内航道多而繁忙，岛屿狭道效应和波流作用明显，涨落潮时刻流速变化较大，对船舶航行造成一定的影响。故针本文确定的2 处事故源项，选择具有代表性的码头前沿以及航道处作为溢油点。同时，选择3 个较为重要且典型的敏感目标作为对象，分析溢油影响及开始影响时间。选取的敏感目标分别为双礁与黄牛礁特殊利用区、西蟹峙特殊利用区及沈家门农渔业区(溢油点位置和敏感区域如图6 所示)。

图6 溢油点位置和敏感区域Fig.6 Location of oil spill point and peripheral sensitive targets

2.2.2 预测条件

根据油膜在潮流作用下输移的物理过程，选择高平潮和低平潮2 个时刻分别进行溢油释放计算[16]。考虑对敏感区域的不利风向(SE)，风速均按建议值取13.8 m·s-1。基于以上条件，根据《船舶污染海洋环境风险评价技术规范(试行)》[17]，码头前沿选择2 t 溢油量进行模拟，航道处由于缺少工程区主力船型数据，取建议值10 t 进行模拟。

2.2.3 溢油工况设计

综合考虑潮流和不利风向因素，对各溢油点按风向和潮流类型进行组合，溢油计算工况如表2、3 所示。

表2 码头前沿溢油预测方案Tab.2 Oil spill prediction scheme for wharf front

表3 航道处溢油预测方案Tab.3 Oil spill prediction scheme at waterway

3 预测结果分析

3.1 码头前沿溢油数值分析

油膜初期受潮流作用明显，在涨潮流作用下迅速向杭州湾西侧方向漂移。随后油膜受往复流作用向东南方向震荡漂移，7 h 后开始对西蟹峙特殊利用区产生环境影响，24 h 内大量油粒子集中扩散到西蟹峙特殊利用区东南侧，随后在SE 风向作用下沿杭州湾方向运动趋势明显，其油膜面积逐渐增大，并在30 h 内开始对双礁与黄牛礁特殊利用区环境有不利影响。最终在不利风向SE 的主导下，72 h 内油膜影响范围达605.2 km2，如图7、8 所示。

图7 溢油发生24 h 后油粒子轨迹(涨潮时刻)Fig.7 Oil particle trajectory after 24 h of oil spill(rising tidal)

图8 溢油发生72 h 后油粒子轨迹(涨潮时刻)Fig.8 Oil particle trajectory after 72 h of oil spill(rising tidal)

落潮时刻溢油在不利风(SE)条件下的预测如图9、10 所示：油膜初期同样在潮流作用下迅速向东南侧方向漂移。大部分油粒子12 h 内往西北飘移，故到达双礁与黄牛礁特殊利用区比涨潮时刻早7 h，对西蟹峙特殊利用区比涨潮时刻晚1 h。油膜在SE 风向以及往复流作用下向西北——东南方向扩散面积逐渐增大，72 h 内溢油扫海面积内影响范围为926.3 km2。

图9 溢油发生24 h 后油粒子轨迹(落潮)Fig.9 Oil particle trajectory after 24 h of oil spill(falling tide)

图10 溢油发生72 h 后油粒子轨迹(落潮)Fig.10 Oil particle trajectory after 72 h of oil spill(falling tide)

3.2 航道处溢油数值分析

航道处溢油初期情况与码头处有所不同，油膜在涨潮流作用下会迅速沿岸线向西北侧漂移，14 h 后油膜到达西蟹峙特殊利用区，前期24 h 内油粒子受岸线吸附影响明显，呈现出长带状油膜，并在27 h后逐渐对双礁与黄牛礁特殊利用区产生环境影响。随后在不利风向和往复流双重作用下，油膜向西北-东南方向震荡漂移，扩散面积不断增大，最终72 h 内油膜影响的海域面积为1 442.3 km2，如图11、12 所示。

图11 溢油发生24 h 后油粒子轨迹(涨潮)Fig.11 Oil particle trajectory after 24 h of oil spill(rising tidal)

图12 溢油发生72 h 后油粒子轨迹(涨潮)Fig.12 Oil particle trajectory after 24 h of oil spill(rising tidal)

航道处不利风(SE)条件下落潮时刻溢油预测如图13、14 所示，油膜初期在落潮流作用下迅速向东南侧方向漂移，与涨潮时不同的是岸线吸附作用减弱，油粒子多是聚集状态，8 h 后开始影响西蟹峙特殊利用区。同时12 h 内大量油粒子聚集在双礁与黄牛礁特殊利用区北侧，故影响时间早于涨潮时刻3 h。油膜随后在SE 风向以及往复流作用下向西北-东南方向震荡漂移，其油膜面积逐渐增大，最终72 h 内溢油扫海面积为1 010.1 km2。

图13 溢油发生24 h 后油粒子轨迹(落潮)Fig.13 Oil particle trajectory after 24 h of oil spill (falling tide)

图14 溢油发生72 h 后油粒子轨迹(落潮)Fig.14 Oil particle trajectory after 72 h of oil spill(falling tide)

综上所述，由于流况不同，风速风向也有一定变化，所以在不同地点、不同时刻发生溢油后所追踪到的油膜轨迹就不尽相同[18]。对于码头前沿处溢油情况，相较于涨潮时刻风向的主导作用，复杂海岛地形也会影响油膜扩散轨迹。其中，金塘岛会对涨潮时刻油膜向杭州湾方向运动起到一定的阻碍作用，扩散动力减弱导致72 h 内油膜影响范围小于落潮时刻；航道处溢油风险预测表明，落潮时刻24 h 内油粒子受潮流与风向作用，多在舟山本岛与金塘水道内漂移，向西南外海区域漂移较弱导致油膜影响面积小于涨潮时刻。同时受水深的影响，航道处水动力条件强于码头前沿，故整体上航道处发生溢油后的影响区域面积大于码头前沿处溢油面积。

4 结论

本文在二维水动力模型基础上耦合油粒子模型对舟山定海某码头前沿及航道处进行溢油预测分析，结果表明：

(1)初期阶段2 处溢油点主要受到潮流作用，沿涨落潮方向进行漂移。随后在不利风向SE 主导作用下，进行震荡漂移扩散。72 h 内涨落时刻码头前沿处溢油最大油膜扫海面积达到926.3 km2，航道处溢油最大油膜扫海面积达到1 442.3 km2。

(2)码头前沿处涨潮时刻，复杂的海岛地形对溢油扩散轨迹起到阻碍作用，扩散动力减弱导致72 h 内油膜影响范围小于落潮时刻；航道处溢油发生后，落潮时刻油粒子受众多岛屿的影响，在舟山本岛与金塘水道内漂移，对外海影响面积减小。

(3)该工程区域发生溢油事故影响敏感目标主要为西蟹峙特殊利用区和双礁与黄牛礁特殊利用区，对沈家门农渔业区影响极小。同时在涨落时刻，溢油到达敏感区域时间存在明显差异，码头前沿最大时间差出现在双礁与黄牛礁特殊利用区为7 h，航道处则在西蟹峙特殊利用区为6 h。