国外对于海上溢油预测研究始于20世纪60年代,Fay通过实验得到了油类在水中分别受惯性力、粘性力及表面张力进行扩展运动的公式^[1]。ehr和Elliott通过海上现场实验,考虑了风及水平、垂直扩散系数,对Fay理论及公式进行了修正^[2-3]。80年代中后期,由Johansen和Andunson提出的基于拉格朗日追踪法的"油粒子"模型越来越受到学者们的青睐,被广泛应用于溢油模拟中。

我国海上溢油研究起步较晚,始于20世纪90年代。赵文谦等^[4-5]在溢油模型中加入风场来研究溢油过程。杨庆霄^[6]、赵云英等^[7]在溢油模型中考虑了油的溶解过程和乳化过程,更全面的对溢油进行了研究。管卫兵^[8]、杨红^[9]先后基于拉格朗日粒子追踪法及自由走动理论建立了长江口溢油模型,分别就风场对溢油影响、水源地污染等问题开展了研究。宋泽坤^[10]、李曰嵩等^[11]基于"油粒子"模型,对油膜厚度变化、漂移轨迹及扫海面积等参数进行了模拟。渤海湾是溢油研究热点,黄娟^[12-13]、娄安刚^[14-15]等对石油钻井平台进行了溢油预报。除此之外,黄成^[16]、吴亚楠^[17]及顾恩慧^[18]等人分别对广西近海、蓬莱海庙港、福建乐清湾等地区进行了溢油研究。国内对于溢油研究研究区域主要分布在沿海河口及港湾、海上钻井平台和海上事故发生地,对于内湾溢油研究较少。

福建霞浦东吾洋地区主要从事海水养殖、内海捕捞和贸易等渔业生产和服务,随着当地海洋渔业经济的迅速发展,渔业生产规模逐年扩大,水产品的卸港量和养殖饵料的吞吐量也大幅度增加,故建设配套的渔港码头是很有必要的。关门地处福建霞浦县东吾洋北部,位于东安岛西侧,紧邻关门水道,自古以来都是进出东吾洋的必经关口。渔港码头在建设和运营的时候存在着较大的溢油风险,一旦发生溢油事故油类会随涨潮流进入东吾洋湾内,东吾洋为半封闭海湾,水体交换能力相对较弱,不利于油类的清除降解;若油类随落潮流进入三沙湾,该区水动力比较大,若不及时采取有效措施,溢油将有可能扩散至外海造成更大影响。为了能在溢油发生后第一时间迅速而有效地做出溢油事故应急反应,本文基于MIKE21软件对关门渔港附近海域潮流进行数值模拟,在此基础上通过"油粒子"模型对该区进行溢油模拟,并分析不同情况下溢油油膜的漂移轨迹、扫海面积以及到岸时间等。

1   材料与方法

1.1   水动力模型控制方程

MIKE21 HD是一个可用于河口、海湾以及近海潮位及潮流的综合数值模型,对于二维流场模拟已经被应用于许多研究之中。本文使用该模块对福建东吾洋进行了水动力模拟实验,得到潮流数据验证良好后输入溢油模块,作为溢油水动力场的驱动条件。水动力模型用连续方程和动量守恒方程来描述流场和潮位的变化:

连续方程:

动量守恒方程:

式中:x,y为直角坐标系坐标;t为时间;h为水深(基准面到底床的距离);ξ为潮位(基准面到自由水面的距离);u,v分别为x,y方向上的垂线平均流速分量;f为科氏力系数;g为重力加速度;E_x,E_y分别为x,y方向上的水平涡粘系数;τ_bx,τ_by分别为波流共同作用下床面剪切力在x,y方向上的分量;S_xx,S_xy,S_yy分别为各向的波浪辐射应力。

1.2   溢油控制方程

1.2.1   水平运动

水平运动是水表面油膜输运的主要过程,其主要动力为潮流及风力,因此在模拟溢油之前,水动力场的精确模拟是十分必要的。在潮流和风力共同的作用下,油粒子的中心位置经Δt 时间后漂移到新的位置,油粒子坐标位置变化由如下公式计算:

式中:X₀和Y₀为油膜初始时刻坐标; u、v为流速; θ为风向角。油粒子漂移速度由下式权重公式计算:

式中:U_W为水面以上10 m处风速; U_S为表面流速;C_w 为风漂移系数,一般取0.03~0.04,本文模型中取其平均值为0.035。由于不同时刻流场及风场不同,故在同一地点不同时刻发生的溢油其漂移轨迹及漂移速度会存在差异。

1.2.2   紊动扩散

本文用随机走动方法进行模拟该过程,假定油膜在水平各方向上扩散同性,随机扩散距离可用如下公式计算:

式中:[R]_-1¹为-1到1的随机数;S_α为扩散距离;D_α为α方向上扩散系数;Δt_p为扩散时间。

1.2.3   扩展运动

溢油在进入水体后受重力、惯性力、粘性力和表面张力作用迅速向四周扩展,形成大面积油膜,此过程中油膜厚度不断减小。对该运动的模拟采用修正的Fay重力－粘力公式:

重力扩展阶段:

粘性力扩展阶段:

表面张力扩展阶段:

式中:g 为重力加速度;V为溢油总体积;ρ_w 为水的密度;v_w 为水的运动粘性系数;σ为空气和水、空气和油、油和水表面张力系数之差;K₁、K₂、K₃为待定系数。

1.2.4   蒸发

油类蒸发指液态的石油烃变成气态造成油膜与空气间进行物质交换的过程,油类的蒸发受油品本身及环境因素影响,油品本身性质决定着油类在水中的蒸发率,油膜厚度、风场及温度、水温、风速计太阳辐射等因素也影响着油类的蒸发。蒸发率由下式表达:

式中:N为蒸发率;k_e为物质输移系数;P^SAT为蒸汽压;R为气体常数;T为温度;M_i为油类的各组分分子量;ρ_ι为油类各组分密度。k_ei由下式估算:

其中:k为蒸发系数;Sc_i为各组分的蒸汽Schmidts数。

1.2.5   乳化

油膜的乳化包括形成水包油型及油包水型两个物化过程。在溢油发生最初的时刻,扩散是最主要的运动方式,由于其是一种物理过程,油膜在自身伸展压缩运动及水流和风浪的扰动下破碎,将油膜分散成各种粒径的油滴,油滴进入水中,形成水包油型乳化。从油膜中扩散到水中的油量损失可由下式计算:

式中:W₁为油类进入到水中的分量;W₂为进入到水中后没有返回的分量;μ_oil 为油类的粘度;γ_ow水油界面的张力。

当水进入油类中,形成油包水型乳化,会增加油类的含水率。油类含水率变化由下式得出:

式中: R₁为吸收速率;R₂释放速率;y_w为实际含水率;A_s为油类中沥青的质量分数;W_ax 为油类中石蜡的质量分数;K₁为吸收系数;K₂为释放系数。

1.2.6   溶解

本文中油类的溶解速率由下式计算:

式中:C_i为各组分的溶解度;X_i 为各组分的摩尔质量分数;M_i 为各组分的摩尔质量;Ks_i为溶解传质系数。

油类在水中运动受油品本身、水动力环境、气象环境及生物环境影响,各因素相互作用下发生复杂的物理、化学及生物反应,其运动呈现出较大的随机性及差异性。除上述过程外,还存在油膜与大气、油膜与水体间的热量迁移以及太阳辐射、反射、蒸发热损失等一系列热量迁移过程^[19],这些因素对于溢油初期油膜运动影响并不大,故本文中中设定油膜与外界环境不存在热量迁移。

1.3   水动力计算参数及网格

模型计算区域为26.3°~26.9°N,119.5°~120.1°E。将其划分成非结构三角网格,并在东吾洋内进行了加密处理,在关门渔港附近分辨率最大达到10 m,整个计算区域包括三角网格16735个,网格节点数27854个(见图 1)。开边界设置在三沙湾口,共计节点12个,在图 1中用红点表示。开边界上使用MIKE21模型自带全球潮位预测值进行驱动,闭边界采用干湿边界处理,取h_dry=0.005 m,h_flood=0.05 m,h_wet=0.1 m。模型使用冷启动,时间步长为3600 s,每小时输出潮流场数据,计算45 d,取后30 d稳定流场进行验证分析,得到该区水动力场。

图  1  模型设计网格及关门渔港所在区域放大

Fig.  1  The grid of model andenlarged view of the fishing port

下载: 全尺寸图片 幻灯片

1.4   溢油模型参数设置

当地渔民进出东吾洋依靠关门水道及大门水道,关门水道来往船只较为密集,加之渔港的建设,使关门成为渔船避风及卸货补给的集散地,使渔船发生碰撞的概率大大增加,该区可能发生的溢油事件是到港船舶碰撞后燃料油泄漏。不同船型、不同事故情况的溢油量都不同,具有较大的随机性。此次模拟的溢油油品设定为燃料油,密度为920 kg/m³,最大含水率为0.85,溢油量取渔港靠泊最大船型携带的单舱燃油量:5 t,溢油类型为固定点源瞬时溢油,油粒子数为1000。在本工程营运时,溢油最大风险可能发生位置位于渔港西侧近关门水道处,因此在渔港西侧关门水道附近设置溢油点,坐标为:119.90°E,26.68°N。由于东吾洋四周环山,风速受山体摩擦及阻挡而减小,该区夏季盛行东南风,冬季盛行西北风,根据现场实测风速及霞浦县志记载^[20],模型中夏季风速取2.5 m/s,冬季风速取3.5 m/s。本文模拟了夏季盛行风、冬季盛行风及静风状态下涨急、落急、涨憩、落憩时刻分别发生溢油的共计12种情景,具体溢油风况及潮时组合见表 1。利用水动力结果作为溢油模型驱动,考虑风况及溢油潮时组合,溢油模型时间步长为3600 s,每小时输出溢油结果,模拟两个潮周期(24 h)内,油膜漂移运动。

表  1  溢油风况及潮时组合

Tab.  1  Conditions of wind regimeand tide

下载: 导出CSV 

| 显示表格

2   结果与讨论

2.1   模型验证

潮位验证数据为福建罗源迹头站及帮门两个验潮站的实测潮位,小潮验证数据为两站2014年9月18日至9月19日实测潮位,大潮验证数据为两站2014年中9月24日至9月25日实测潮位,流速及流向验证数据为浒玙澳及仙东2014年9月24日的现场实测资料,潮流验证站及模型水深分布见图 2。根据验证结果可知,该模型模拟所得潮位及流速流向与实测数据验证较为吻合,可较好地反映当地水动力形态,计算所得水动力场可用于溢油模型(图 3~图 6)。

图  2  潮流验证站及水深分布

Fig.  2  Distribution of tidal and current verification station and water depth

下载: 全尺寸图片 幻灯片

图  3  帮门大潮、小潮计算水位与实测水位的比较

Fig.  3  Comparison of computed and observed water levels at Bangmen station

下载: 全尺寸图片 幻灯片

图  4  罗源迹头大潮、小潮计算水位与实测水位的比较

Fig.  4  Comparison of computed and observed water levels at uoyuanjitou station

下载: 全尺寸图片 幻灯片

图  5  浒玙澳流速与流向计算值与实测值的比较

Fig.  5  Comparison of computed and observed current at Huyuao station

下载: 全尺寸图片 幻灯片

图  6  仙东流速与流向计算值与实测值的比较

Fig.  6  Comparison of computed and observed current at Xiandong station

下载: 全尺寸图片 幻灯片

2.2   水动力模拟结果

油类在海面漂移,表层流场为其主要的驱动力,图 7为东安岛附近大潮期涨落潮表层流场分布图,东吾洋内岛屿较少,潮波在东吾洋内传播较少出现变形,潮流运动方式单一,主要为往复流。涨潮时,外海潮流由关门水道进入东吾洋湾内;落潮时潮流由东吾洋流向外海,受东安岛地形岸线[]阻挡,东安岛附近潮流分流、漩涡较明显,且涨潮流稍大于落潮流。

图  7  关门附近大潮涨急(左)、大潮落急(右)表层流场分布

Fig.  7  The surface currentof Guanmen in spring tide at the moment of up (left) and fall (right) anxious

下载: 全尺寸图片 幻灯片

2.3   溢油模拟结果

2.3.1   静风状态下,不同潮时溢油

选取静风状态下,落急、落憩、涨急、涨憩时刻(情景1、情景4、情景7、情景10)溢油进行对比分析,由图 8可知,不同时刻溢油,其漂移轨迹存在一定差异。由于风场为静风,故油膜仅在潮流的作用下进行漂移运动,其轨迹与潮流流向较为一致。

图  8  静风状态下各潮时溢油轨迹

Fig.  8  Oil particle tracks for each tide under clam condition

下载: 全尺寸图片 幻灯片

静风状态下各潮时到岸时间及溢油漂移路程、扫海面积随时间变化见表 2。到岸时间从小到大为:落急＜涨急＜落憩＜涨憩,涨落急时油膜刚开始漂移速度较快,因关门水道较为狭长,两侧岸线距离较近,油膜容易随潮流运动到岸边,并吸附于岸线之上。溢油事故发生后,油膜的漂移路径与溢油发生时刻及风况有密切关联。静风条件下溢油,油膜运动主要受潮流影响,漂移轨迹与潮流运动方向较为一致,流速变化对于油膜漂移速度变化起了决定作用。油膜经过一定时间的漂移运动后,有较大可能到达东安岛西北部岸线,并对当地的生态环境产生一定的影响。扫海面积随漂移路程增大而增大,静风状态下各潮时溢油扫海面积从小到大依此为:涨憩＜落急＜涨急＜落憩,落憩时刻油膜随落潮流最远漂移至东吾洋与三沙湾交界处,其漂移距离较其他三种情况更大,故扫海面积最大。

表  2  静风状态下各潮时到岸时间(h)及溢油漂移路程(km)、扫海面积(km²)随时间变化

Tab.  2  The reaching shoreline time(h),drift distance(km),swept area(km²) changing with time for each tide under clam condition

下载: 导出CSV 

| 显示表格

2.3.2   相同潮时,不同风况下溢油

各潮时在不同风况下溢油情况大致相同,油膜受风况及潮流共同驱使,油膜漂移方向为两者矢量和方向,若潮流与风向一致,油膜漂移速度较快;若相反,其漂移速度较慢。落憩时刻油膜漂移受风况及潮流作用较明显,3种不同风况下油膜漂移路径差异最大,因此本文选取落憩时刻3种风况溢油(情景4、情景5、情景6)作为相同潮时不同风况溢油的典型进行比较,漂移轨迹如图 9所示。

图  9  落憩时不同风况溢油轨迹

Fig.  9  Oil particle tracks for each windregimeat ebb moment

下载: 全尺寸图片 幻灯片

落憩时不同风况溢油到岸时间及漂移路程、扫海面积随时间变化见表 3。油膜在潮流及风场的叠加作用下进行漂移运动,油膜漂移轨迹为流场与风场的合力方向,当流场作用大于风场时,油膜漂移轨迹偏向于流场一侧;反之,则偏向风场方向。落憩冬季盛行风情况下,油膜漂移轨迹受风力影响漂移轨迹较静风状态下偏向东南方向。落憩夏季盛行风下,油膜受夏季盛行风影响,在关门水道西北方岸线附近来回移动,到岸时间较其他[]两种情况短。落憩时溢油扫海面积从小到大依次为:夏季风＜冬季风＜静风,由于夏季风对油膜的挤压作用,油膜漂移路程较短,故其扫海面积也相对较小。油膜在冬季风作用下到达岸线的时间短于静风状态下溢油,故其漂移路程与扫海面积小于静风状态溢油。

表  3  落憩时不同风况溢油到岸时间(h)及漂移路程(km)、扫海面积(km²)随时间变化

Tab.  3  The reaching shoreline time(h),drift distance(km),swept area(km²) changing with time for each wind condition at ebb moment

下载: 导出CSV 

| 显示表格

3   结 论

利用MIKE21 HD模块建立了福建霞浦东吾洋海域水动力模型,较好的模拟了该海域的潮流情况。在此基础上,利用拉格朗日"油粒子"模型对关门渔港附近溢油情景进行了模拟,由数值模拟结果可得到以下结论:

(1) 油类在水中运动受多种因素共同影响,不同潮时、不同风况下溢油的漂移轨迹、漂移路程及扫海面积都存在较大差异。

(2) 关门渔港附近溢油漂移扩散主要受潮流和风场影响,其中风场起主导作用。24 h内油膜最大扫海面积及漂移路程均出现在落憩静风时溢油,分别为6.66 km²及17.63 km。关门渔港附近溢油,油膜大都在关门水道附近来回移动,并有较大可能到达东安岛西南侧弥勒湾内,对该区生态环境会造成一定影响。

(3) 静风状态下溢油,油膜受潮流驱动,漂移轨迹与潮流流向较一致;不同风况下溢油,油膜受风及潮流共同驱使,油膜漂移方向与两者矢量和方向大致相同。