随着沿海石油贸易、港口建设、航运规模的发展，海上船舶的航行密度增大，船舶间或船舶与其他物体间存在碰撞并发生溢油事故的可能性加大。据统计，1990—2010年我国50 t以上船舶溢油事故共发生71起，溢油总量为22035 t，船舶事故主要发生在开阔海域和港口/码头，污染物类型主要为燃料油和柴油^[1]。溢油事故不仅会对海洋生态系统造成威胁，还会对渔业、海岸环境以及人类的可持续发展造成不可估量的严重后果^[2]。世界各国相继开展了溢油模拟的研究，发展了许多溢油模式，模拟和预测海上溢油的时空分布，为溢油事故发生后的应急措施和有效控制提供参考依据，减少事故造成的损失。

溢油模拟的研究经历了：油膜扩展理论、对流扩散方法和油粒子法三个阶段。油粒子法在溢油模型中得到广泛应用，最早是由Johansen^[3]和Elliot^[4]提出，该方法将溢油离散成大量的油粒子，每一个油粒子代表一定油量，通过拉格朗日法追踪油粒子的运动，模拟油膜的漂移扩散过程。刘伟峰^[5]对油粒子法进行了改进，提出了溢油运动的两阶段模拟法。该方法将溢油运动过程分为自身扩展和紊动扩散两个阶段, 前一阶段根据Fay理论修正模式计算, 后一阶段采用油粒子方法模拟, 通过“油膜粒子化”技术将两阶段进行衔接，比油粒子方法更能真实地反映溢油的运动规律。王永学和郭为军等^[6]基于欧拉-拉格朗日理论, 提出了一种混合溢油模型，用拉格朗日方法来模拟溢油在海洋表面的运动过程，而用欧拉法预测溢油在水柱中的分散问题。吴丹^[2]建立了包括全动力轨迹模型和溢油扩展模型的海洋溢油数学模型, 同时考虑油膜扩展的随机性, 在扩展量上加了一个随机量, 以此来描述油膜漂移轨迹的不规则形状, 并将之应用于渤海海域。袁春光等^[7]将基于欧拉-拉格朗日理论体系的溢油模型应用于大亚湾水域，分别模拟了移动点源和固定点源泄漏。李大鸣等^[8]考虑了风场和流场对溢油油膜的扩展及漂移的影响，并借助水槽试验验证油膜扩展公式的正确性，建立了海上溢油数学模型对栾家口港区在不同动力条件、溢油点、溢油方式及油品种类分别进行模拟。杨红等^[9]基于GNOME溢油模型模拟长江口外锚地发生溢油时的漂移轨迹与扩散情况。国外较为成熟的溢油全动力模型有美国ASA研究所开发的OILMAP模型及Reed等(1999)建立的OSCAR模型^[10]。丹麦水环境研究所开发的溢油模型(MIKE OS)在工程研究中被广泛应用，该模型采用“油粒子”法来模拟海上溢油的行为和归宿，添加了油粒子轨迹模块和岸线作用及围油栏、拦油坝等功能，可以得到油粒子轨迹、油膜范围和厚度等参数，能够考虑多种溢油种类、泄漏体运动形式、以及溢油泄漏形式。齐庆辉等^[11]运用MIKE溢油模型模拟了京杭运河苏北段溢油扩散，许婷等^[12]将之应用于厦门港海域的溢油扩散数值模拟，冯静^[13]采用MIKE溢油模型模拟了滨海近岸溢油风险及对附近环境敏感区的影响。

浙江中澳现代产业园项目主要布局于舟山群岛新区金塘北部围垦区，澳洲活牛进境加工项目作为该产业园的启动项目，其配套码头工程位于金塘岛北部鱼龙山-横档山新建围堤的中部，码头东侧临近西堠门水道，与册子岛隔水相望，项目所处位置重要，运营期间一旦发生溢油事故，对该海域环境及周边产业将造成严重影响。本文采用MIKE21中的HD模块建立金塘港海域的潮流数学模型，并采用实测水文资料对模型进行修正和验证，进而利用MIKE中OS模块建立溢油模型，对溢油扩散进行数值模拟，预测分析溢油扩散对周边海域的影响，为金塘港港区溢油事故的应急措施提供有效的参考依据。

1   材料与方法

1.1   溢油模型概述

油类进入水体中，其行为和归宿受风、浪、流等动力环境因素、非环境动力因素和油类的物理化学特性等多种因素的支配。油类的变化过程主要为油膜的扩展、漂移和风化过程。油膜的扩展包括在其自身重力、表面张力和惯性力作用下的自身扩展运动。油膜的漂移是其在流场和风场作用下的漂移扩散运动。风化过程主要是指油类在蒸发、溶解、乳化、生物降解引起油类的物理和化学特征发生变化。MIKE OS模块，主要采用油粒子法，来模拟溢油在水体中的扩展、漂移和风化过程，从而模拟预测油粒子的漂移轨迹、扩散范围、浓度分布、厚度分布等特征的时空变化。

1.2   溢油模型控制方程

1.2.1   扩展过程

油膜的扩展计算采用修正的Fay重力-粘力公式，公式如下：

式中：A_oil为油膜面积；R_oil为油膜直径；A_oil=πR_oil²；K_a为系数；t为时间；V_oil为油膜体积；h_s为初始油膜厚度，V_oil=R_oil²·π·h_s, h_s=10 cm。

1.2.2   输移过程

输移过程指的是油粒子在风、浪、流等综合作用下的平移运动，包括漂移运动和紊流扩散。

(1) 漂移运动

油粒子的漂移运动受水流作用力和风拽力的支配，决定了油粒子的运动方向和范围分布，油粒子的漂移运动速度方程如下：

式中：U_t为总的漂移运动速度；c_W为风漂移系数，取值范围0.03~0.04；U_W为水面以上10 m处的风速；U_s为表层流速。

在Δt时间内的油粒子的位置变化方程如下：

式中：X₀, Y₀为油粒子的初始时刻位置; u, v为流速分量；θ为风向角。

(2) 紊动扩散

油粒子的紊动扩散是一种湍流弥散过程，具有很大的随机性。假定水平扩散各向同性，油粒子的一个时间步长内方向的随机扩散距离S_α的计算公式如下：

式中：[R]_-1⁺¹为-1~1间的随机数；D_α为α方向的扩散系数；Δt_p为扩散时间。

1.2.3   风化过程

油粒子的风化过程主要包括蒸发、溶解和乳化等，在此过程中油粒子的物理和化学特征发生变化，但水平位置不发生变化。

(1) 蒸发

油粒子的质量传输主要通过蒸发来完成，蒸发过程受油组分、温度、油膜厚度和面积、风速和太阳辐射等因素的影响，其中轻组分蒸发较快，重组分蒸发较困难。油粒子的蒸发率计算公式^[10]如下：

式中：EVAP为蒸发率；K_e为物质输移系数；P_vp为蒸汽压；R气体常数；T为温度；f为油组分的蒸发分数；M为油组分的摩尔重量。

物质输移系数K_e的计算公式如下：

式中：k为蒸发系数; S_c为蒸汽施密特数。

(2) 溶解

油中的可溶性部分在海水中发生溶解，油膜面积相应减少。溶解受油膜面积、风、海面条件、气温以及日照强度等因素的影响。在模型中，将油类分为挥发性组分和重组分两部分分别计算溶解率，计算公式^[10]如下：

式中：k_v, k_h分别为挥发性和重组分的溶解率；M_v, M_h分别为挥发性和重组分油粒子质量；M_oil为油粒子总质量；ρ_v, ρ_h分别为挥发性和重组分的密度；f_d为化学分散剂效果；C_v, C_h分别为挥发性和重组分的溶解度。

(3) 乳化

溢油发生若干个小时后，油和水混合在一起形成油水乳化物，该过程取决于油膜的厚度、溢油本身的密度、黏度及风浪大小等因素。在模型中乳化看作是油包水和水包油两个阶段的动态平衡过程，可用含水率的变化来表示乳化过程，计算公式如下：

式中：y_W为实际含水率；R_u, R_r分别为水的吸收速率和释放速率；k_u, k_r分别为吸收系数和释放系数；y_W^max为最大含水率；Wax为油中的石蜡含量(重量比)；μ_oil为油粘度。

1.3   水动力模拟

1.3.1   模型方程

MIKE21 HD模块采用非结构化网格，运动方程为垂向平均的二维浅水方程。

(1) 水动力模型连续方程如下：

(2) 水动力模型动量方程如下：

式中：ξ为潮位；t为时间；x, y为直角坐标系坐标轴；u, v为x, y方向上的垂线平均流速分量；h为水深；N_x, N_y为x, y方向的水平紊动粘性系数；f为科氏力参量；f_b为底部摩擦系数。

1.3.2   模型建立

工程所处水域位于舟山、金塘、岱山诸大岛间，形成包括深槽、边滩、潮流沙脊等地貌单元的峡谷型潮汐通道, 受岸线及地形影响，潮流基本为往复流形式。给定潮位开边界，使用干湿判别法对水陆交界、防波堤及围堰等进行处理。考虑到边界条件的选取，计算域东西方向长约474 km，南北方向长约452 km，包括整个杭州湾(如图 1所示)。模型采用三角形网格，在金塘岛附近进行局部加密(如图 2所示)。

图  1  计算域范围

Fig.  1  Calculational area range

下载: 全尺寸图片 幻灯片

图  2  网格划分

Fig.  2  Model grids

下载: 全尺寸图片 幻灯片

1.3.3   模型验证

模型验证采用2017年10月的实测数据，包括一个潮位测站(T1)和四个潮流测站(D1~D4)，以上测站的位置如图 3所示(图中白色圆点)。潮位观测期为一个月，观测结果表明该海区潮汐性质属于“规则半日潮”类型，浅海分潮显著，平均潮差为2.55 m，最大潮差为3.70 m，最小潮差为0.82 m。4个潮流观测站，进行了15 d的连续观测。涨落潮流主要呈往复流特性，也有一定的旋转流特征，大、中、小潮的平均流速分别为0.49 m/s、0.46 m/s、0.30 m/s，实测期涨、落潮平均流速分别为0.40 m/s和0.58 m/s，涨潮流速小于落潮流速。

图  3  水文观测站布置(白色)/溢油位置(红色)

Fig.  3  Hydrological observation distribution(white)/oil spill locations(red)

下载: 全尺寸图片 幻灯片

潮位验证如图 4所示，潮位实测值与计算值吻合较好。流向流速验证如图 5、图 6所示，实测值与计算值拟合良好。以上表明二维潮流数学模型能较好地反映工程海域潮流场的时空分布，可以进一步为溢油模拟提供必要的水流动力条件。

图  4  潮位验证

Fig.  4  Tidal level verification

下载: 全尺寸图片 幻灯片

图  5  流速验证

Fig.  5  Current velocity verification

下载: 全尺寸图片 幻灯片

图  6  流向验证

Fig.  6  Current direction verification

下载: 全尺寸图片 幻灯片

2   结果与讨论

海上溢油一般由于船舶间碰撞或触礁导致船体破损，油从破损处泄漏流入海中。本项目为澳牛运输，考虑为3万吨级运输船舶，溢油源强为375 t的燃油，泄漏点为航道处和码头前沿的靠泊平台，具体溢油位置如图 3所示(图中红色圆点)。溢油方式为固定点源瞬时溢油，油品为燃料油，油密度取970 kg/m³。在模拟中，将燃油分为四个组分：挥发性组分、重组分、蜡和沥青，其中蜡和沥青不溶解于水，不发生降解，也不蒸发。

项目所在海区属亚热带季风气候区，同时又属于大陆性气候向海洋性气候的过度区，风向有明显的季节性变化。秋冬两季受大陆性气候影响，冷空气活动频繁，风向以偏北风为主；春夏两季主要受海洋性气候影响，太平洋暖流比较活跃，风向以南到偏南风为主。本文中溢油计算工况考虑静风、冬季主导风、夏季主导风以及不利风向，分为涨潮和落潮时刻发生溢油。根据金塘测风站2003年11月—2004年10月的10 min平均风速资料，选用07—10月风速资料作为夏季主导风，风向为SE，风速为7.07 m/s，11—02月风速资料作为东季主导风，风向为NNW，风速为7.17 m/s。工程区的常、强风向为NW和SW，风速取船舶的最大作业风速6级风(10.8~13.8 m/s)下限值10.8 m/s^[14]。综合上述条件，本次研究选择8个工况进行分析，选择工况如表 1所示。不同工况下，重点关注油膜在24 h内海域的影响范围，预测分析溢油事故发生后油膜在24 h内的扩散输移方向和主要分布，为金塘港港区溢油事故的应急措施提供有效的参考依据。

表  1  溢油事故发生后24 h内扫海面积 (单位：km²)

Tab.  1  Sweeping area within 24 h after oil spills (Unit:km²)

下载: 导出CSV 

| 显示表格

溢油事故发生后，油膜的扩散输移与燃油泄漏地点、溢油发生时刻、流场特征以及风速风向等因素相关。经过计算得到溢油事故发生后24 h的扫海面积见表 1，油膜分布见图 7和图 8。

图  7  航道溢油24 h扫海范围

Fig.  7  Sweeping area within 24 h when oil spills occur on channel

下载: 全尺寸图片 幻灯片

图  8  靠泊平台溢油24 h扫海范围

Fig.  8  Sweeping area within 24 h when oil spills occur on berth platform

下载: 全尺寸图片 幻灯片

(1) 航道和靠泊平台溢油点相距较近，溢油发生后24 h内油膜的扫海面积和漂移路径相似，无明显差异。

(2) 无风时，油膜运动由涨、落潮流主导，输移过程中受涨落潮影响随水流流向出现摆动。涨潮时溢油，油膜最初随涨潮水流向NW向扩散，随后部分油膜随落潮水流向SE向扩散，最终油膜主要分布在溢油点的NW和SE两个方向；落潮时溢油，油膜的扩散范围和输移规律同涨潮相似，其24 h后的扫海面积小于涨潮时溢油。靠泊平台涨潮时溢油，24 h的扫海面积最大为195.74 km²。

(3) SE风，油膜的运动主要受风向控制，长期扩散方向与风向一致，向杭州湾开敞水域输移，不受岸线和岛屿的阻挡，扩散范围更广。涨潮时，水流流向与风向相近，油膜始终往NW向扩散漂移；落潮时，水流流向与风向相反，溢油初期，油膜受水流支配，随落潮水流流向SE向。航道涨潮时溢油，24 h的扫海面积最大为209.71 km²。

(4) NW向风速最大，油膜的扩散速度更快和漂移距离更大，在NW向风的作用下，油膜在金塘岛南侧岸线和周边岛屿沿岸堆积。落潮水流流向与风向相近，油膜扩散范围更广。航道落潮时溢油，24 h的扫海面积最大为131.99 km²。

3   结论

本文利用MIKE21 HD模块建立舟山海域二维水动力模型，模拟该海域的流场特征。在此基础上，采用MIKE OS模块模拟溢油油品在海面的扩散漂移过程，预测分析澳洲活牛进境加工项目运营期间可能的溢油事故对周边海域的影响。研究结果表明，溢油事故发生后，油膜的扩散输移与燃油泄漏地点，溢发生时刻，流场特征以及风速风向相关。静风时，油膜运动由涨、落潮流主导，随水流流动扩散输移，最终油膜主要分布在溢油点的NW和SE两个方向。不同风向条件下，油膜的运动主要受风向控制，长期扩散方向与风向一致。其中，SE向风时，油膜主要向杭州湾开敞水域输移，不受岸线和岛屿的阻挡，扩散范围更广；NW向风速最大，油膜的扩散速度更快和漂移距离更大，在NW向风的作用下，油膜在金塘岛南侧岸线和周边岛屿沿岸堆积。