广西是我国西南地区唯一的沿海省份，也是我国西南出海口的重要通道。近年来，随着中国—东盟贸易合作的加强以及广西经济的快速发展，广西沿海港口船舶进出港更为繁忙，船舶碰撞、搁浅等事故发生的可能性大大提高，溢油事故发生的可能性不断加大。同时，广西沿海工业发展迅速，一些大型石油化工落户广西，岸上溢油进入海洋造成污染的可能性加大。广西北部湾是一个半封闭式海湾，水交换能力差，拥有着丰富的海洋资源，拥有红树林、珊瑚礁、海草床等典型生态系统，拥有中华白海豚、儒艮等珍稀濒危物种，一旦溢油事故发生，极有可能对海洋生态、环境造成极大的影响。开展广西近海海上溢油漂移扩散特性研究，构建便捷的溢油应急预报系统，对广西近海海上溢油事故防范、应急处理，保护广西海洋生态环境具有重要意义。海上溢油漂移扩散规律的研究始于20世纪60年代。Blokker[1]通过单组分油建立溢油风化模型，提出油膜扩展理论，但其只考虑重力和溢油量的影响，结果不能真实地反应溢油变化情况。Fay[2]在已有研究基础上作出改进，提出油膜扩展同时受惯性力、重力、黏性力及表面张力4个因素的作用，并将溢油扩展分为惯性力扩展、黏性力扩展和表面张力扩展3个阶段，形成目前最为经典的三阶理论。Mackay等[3]在Fay理论的黏性力扩展阶段中考虑风作用，建立了油膜扩展与风之间的经验公式，该理论后来成为溢油风化研究的基础。Liu等[4]结合Fay溢油三阶段理论，将油膜扩展与水体运动结合，提出油膜扩展与时间之间的变化公式[4]。Lehr等[5]修正Fay理论，建立了油膜变化的椭圆模型，将油膜的扩展在风向上拉长。Johansen等[6]根据拉格朗日跟踪法对“油粒子”模型进行改进，使之可直观地反映油膜在风、波浪等作用下破碎、分离的真实情况。杨红等[7]结合拉格朗日追踪法和随机运动理论模拟长江口南港水道溢油情况，模拟结果与实际情况吻合较好。类似的还有黄娟等[8]和李怀明等[9]学者的研究。目前，溢油漂移扩散理论的研究相对成熟，各界学者也提出较多的溢油模型，如MIKE SA模型、IKU风化模型、美国的ADIOS风化模块及OILMAP模型系统、英国的OSIS风化模块等[10-12]。国内关于溢油扩散机理的研究也取得了较多成果，在渤海、长江口水域、珠江口、北部湾等海域均有较多的应用[13-14]。然而，以往的研究过程中，针对广西近海的海上溢油漂移扩散研究鲜有报道，且目前所采用的主流溢油模型计算复杂，并不适用于应急处理业务化应用。本研究试图构建适用于广西近海的海上溢油应急预报系统，并以钦州湾为例对系统进行应用，为北部湾海上溢油事故防范、应急处理、海洋生态环境保护等提供科学依据。
1 系统结构系统结构主要有模型计算程序和可视化界面两个部分(图 1)。模型计算程序包括模型默认的参数、用户自定义参数，以及核心计算程序代码组成。模型默认的参数主要为模型一些经验参数，如风漂移系数、溢油扩展系数、质量转移系数等，这些经验参数简化了客户端参数输入，使得用户更容易使用系统，同时保证了计算的准确性。用户自定义参数，如溢油时间、溢油量等，根据用户需要输入到系统中，使计算更为接近实际。核心计算程序根据模型默认的参数和用户自定义参数，完成溢油漂移、蒸发、乳化、溶解等过程的计算。
|
| 图 1 系统结构示意图 Fig. 1 Schematic diagram of system structure |
模型可视化界面基于Windows操作系统以及客户端浏览器，主要用于用户输入自定义参数、地图展示以及计算结果输出等。用户将自定义参数输入后，数据将会保存至服务器内，用户发出计算指令后，模型计算程序将根据输入的数据进行计算。计算完成后，服务器将计算结果返回至用户端的模型界面，供用户查看、导出数据。
2 系统主要功能和特点 2.1 系统功能(1) 通过有访问权限的客户端浏览器，用户可方便地将溢油开始时间、结束时间、溢油量等参数输入至系统内。同时，参数输入界面右侧地图展示界面将显示溢油的计算区域，并且在地图上标识出溢油源的位置，方便用户检查溢油位置是否准确。
(2) 根据用户输入的计算参数，模型对假定的溢油情况进行计算。计算完成后，系统将计算结果返回至图形界面，用户可在图形界面中查看漂移轨迹、扩散范围、扩散面积，以及油膜厚度、乳化率等溢油信息，其结果有助于准确寻找油膜的位置和动态，及时将溢油情况提供给应急救援队伍，以利于开展应急救援行动或发出溢油预警和应急部署。
2.2 技术特点与国内其他同类研究和应用成果相比较，本研究系统主要有以下几个方面的技术特点：
(1) 系统内置配置有默认的地形、岸线数据以及大部分模型参数，从而能在保证一定计算精度的同时，降低模型的操作难度，这使得即使是非专业人员也能使用本系统，极大提高了系统的便捷性。
(2) 系统布置于服务器上，将数据处理、计算过程移至高配置的服务器上，用户通过指定网络上具有访问权限的PC端浏览器即可进行操作计算，因此，系统对客户端PC要求较低，无需安装、无需布置，只要打开浏览器操作即可，极大减轻了客户端的压力，同时大大提高了计算速度。
(3) 图形界面采用图像叠加的方式进行图像处理，将底图设为固定的图像预先保存至服务器内，而系统计算时只计算油膜涉及的范围。系统出图时对地图进行调用，同时将油膜与底图行进叠加展示，从而避免了模型计算范围过大而影响计算速度的弊端，极大地提高了系统的计算效率。
3 计算方法 3.1 水动力模型水动力模型主要为溢油模型提供潮流数据，其计算结果的准确性对溢油模型计算结果有显著的影响。系统采用FVCOM海洋数值模式为溢油模型提供水动力场，该模式由美国麻省大学海洋科学技术学院陈长胜博士领导的研究小组与美国伍兹霍尔海洋学协会的罗伯特C比尔兹利博士合作开发，其原始方程主要包含动量方程、质量连续方程以及温度、盐度和密度方程。在物理和数学上用Mellor-Yamada 2.5阶垂向湍流闭合模型及Smagorinsky水平湍流闭合模型对方程组进行闭合。在垂向上适用σ坐标系或通用垂向s坐标系对不规则底部地形进行拟合，在水平上利用无结构三角形网格对水平计算区域进行空间离散。
模型结合了有限元方法的自由几何拟合特性和有限差分方法的离散结构简单及计算高效的特性，从而能综合这两种方法的优点。利用有限体积积分格式，能更好地保证复杂几何结构的河口海湾及海洋计算中的质量、动量、盐度、温度及热量的守恒性。模型输入输出基于标准化NETCDF格式，使得该模型计算结果能方便地被其他计算程序调用，具有较好兼容性，在数值处理方法和岸线地形拟合上的优势使其在河口海岸地带得到了广泛应用。模型求解方程组及具体操作方法见文献[15-16]。
3.2 溢油模型 3.2.1 概述溢油模型采用经典的“油粒子”模型进行计算。“油粒子”模型采用拉格朗日法确定平流过程，并采用随机走动法模拟扩散过程，实现了确定性和随机性的结合。通过追踪油粒子随水体平流和湍流扩散运动的过程，模拟油膜的形状变化和漂移轨迹，同时通过扩展、扩散、溶解、乳化和蒸发来计算溢油在水体表面的扩展和风化过程。
3.2.2 漂移油粒子漂移是潮流运动及风拽力共同作用的结果，其漂移速度与潮流流速及风速有关。油粒子漂移速度可用以下权重公式计算[17]：
| $U_{ {oil}}=V_{s}+C_{ {wind}} \cdot U_{ {wind}}, $ | (1) |
式中，Uwind为水面以上10 m处的风速，Cwind为风漂移系数，一般取值在0.03~0.04，Vs为表面流速度。
3.2.3 扩展油粒子扩展是重力、惯性力、黏性力和表面张力等共同作用的结果，根据修正的Fay重力-黏力公式，油膜扩散计算式为[17]
| $\left( {\frac{{d{A_{oil}}}}{{{d_i}}}} \right) = {K_a} \cdot {A_{oil}}^{1/3} \cdot {\left( {\frac{{{V_{oil}}}}{{{A_{oil}}}}} \right)^{1/3}}, $ | (2) |
其中，Aoil=πRoil2，Voil=πRoil2hs。式中，Aoil为油膜面积，Roil为油膜半径，Voil为油膜体积，t为时间，Ka为扩展系数，hs为油膜厚度。
3.2.4 扩散油粒子扩散与水流紊动有关，根据随机扩散理论，假设各方向扩散具有同一性，一个时间步长内油粒子在a方向的扩散距离可表示为[7]
| ${D_a} = R\sqrt {6 \cdot {C_a} \cdot \mathit{\Delta }t} , $ | (3) |
式中，Da为时间步长内溢油在a方向上的扩散距离，R为-1~1的随机数，Ca为a方向上的扩散系数，Δt为时间步长。
3.2.5 蒸发油膜蒸发速率与油分气温、水温、溢油面积、风速、太阳辐射和油膜厚度等因素有关，油膜蒸发速率可用以下公式表示[18]：
| $\frac{\mathrm{d} Q}{\mathrm{d} t}=-k_{E} A_{ {oil}} X M P / R T, $ | (4) |
其中，kE=-k·Aoil0.045·Sc-2/3·Uwind0.78。式中，
溢油在水流紊动作用下油膜的影响逐渐增大，逐渐被紊动能量撕裂成油滴，形成水包油的粒子，从而形成一种油水相交乳化物。油膜乳化可用油膜含水率来表示，其计算式为[7]
| $Y_{w}=\frac{1}{K_{B}}\left(1-e^{-K_{A} K_{B}\left(1+U_{w i n d}\right)^{2}} t\right) K_{A}, $ | (5) |
其中，KB=1/YWF，YWF为油膜最终含水率，Yw为乳化物的含水量，KA为经验参数，取值约为4.5×10-6。
3.2.7 溶解溢油的溶解主要发生在低碳轻质油组分，溶解量通常很小，其溶解率用下式表示[18]：
| $\frac{\mathrm{d} V_{d s}}{\mathrm{d} t}=K s \cdot C_{d s} \cdot X \cdot \frac{M}{\rho} \cdot A_{oil}, $ | (6) |
式中，
溢油时间选择涨潮和落潮时两个时刻。地点选择钦州港附近海域，溢油事故通常发生于港口航道等航运较繁忙区域。溢油量根据《船舶油污染事故等级标准》中重大溢油事故的划分标准，取为10 t，并假设为瞬时溢出。气象条件根据钦州龙门站历史统计资料，选择年平均风速5 m/s，风向为常方向N向。横向扩散系数取0.10 m/s2，纵向扩散系数取0.25 m/s2。模拟时长设为24 h，并假定油膜最终扩展厚度为0.01 mm[19]。各溢油计算工况及参数如表 1所示，根据上述模型及参数，对钦州湾海上溢油漂移扩散进行模拟。
| 项目Item | 参数设定Parameter setting |
| 溢油量Spill volume | 10 t |
| 溢油位置Oil spill location | 108°37′40″, 21°41′45″, |
| 溢油时刻Oil spill time | 涨潮、落潮Flood tide, ebb tide |
| 溢油方式Oil spill way | 瞬时Instantaneous |
| 最终扩展厚度 Final extended thickness |
0.01 mm |
| 模拟时长Simulation duration | 24 h |
| 风向Wind direction | N、E、S、W、静风 N，E，S，W，static wind |
| 风速Wind velocity | 5 m/s |
| 气温Air temperature | 25℃ |
| 水温Water temperature | 25℃ |
| 纵向扩散系数 Longitudinal dispersion coefficient |
0.25 m/s2 |
| 横向扩散系数 Transverse dispersion coefficient |
0.10 m/s2 |
4.2 溢油轨迹分析
静风情况下，当涨潮时刻发生溢油时，北向涨潮流控制溢油运动，油膜随着潮流向西北偏北方向移动，由钦州湾湾颈中部进入茅尾海。涨潮期间，钦州湾湾颈中部流速较大，油膜漂移至龙门附近时移动速度可超过0.5 m/s，进入茅尾海后潮流运动减弱，油膜漂移速度随之下降，约为0.2~0.3 m/s。由于“油粒子”漂移具有随机性的特点，多次计算发现油膜飘移至钦州湾湾颈时还可能从七十二泾进入茅尾海，使得七十二泾一带岛屿受不同程度的影响。整个计算时段内，溢油向茅尾海方向漂移的距离可达17.5 km之远。带潮流反转后，油膜基本沿涨潮的路径从钦州湾湾颈中部退出茅尾海(图 2)。
|
| 图 2 不同时刻溢油漂移轨迹 Fig. 2 Drift track of spilled oil at different time |
当落潮时刻发生溢油时，油膜受南向潮流的影响向钦州湾外湾方向移动，大致方向沿着深水航道，影响区域主要为深水航道附近，其原因可能与航道处水动力条件好，潮流流速大，有利于物质的输运有关。一个落潮周期内，油膜可向钦州湾外湾方向漂移最远距离约17.9 km，随后潮流涨潮时又从南向北移动(图 2)。
4.3 溢油扩散特征分析溢油进入水体后，油膜在潮流及自身重力的作用下快速扩展，油膜面积扩散速度较快，油膜面积曲线斜率较高，但此时油膜面积并不大。随着时间的推移，油膜面积不断扩大，但油膜面积曲线斜率不断下降，扩散速度越来越慢。当油膜扩散至一定面积，厚度达到设定的下限(0.01 mm)后，油膜面积不再扩大，此时，油膜面积曲线不再上升。
相比于落潮时刻，涨潮时刻发生溢油时相同时间内油膜扩散面积相对要小(表 2)，原因在于涨潮时溢油发生后，油膜向茅尾海方向漂移，经过的七十二泾区域海域面积较小，水流集中，不利于油膜扩散；进入茅尾海后，海域面积扩大，水流分散，油膜扩散得以加强，油膜扩散面积增加较快。落潮时发生溢油，油膜扩散面积不断增大，且比涨潮时稍大(表 2)，这与溢油发生时油膜向更开阔的海域移动，更有利于扩散有关。扫海面积表现为开始时落潮大于涨潮，随着时间推移，逐渐变为涨潮大于落潮(表 2)，原因在于涨潮时发生溢油，油膜经过的七十二泾油膜面积较小，扫海面积并不大，进入茅尾海后，油膜面积扩展加强，从而扫海面积增加比较明显；而落潮时发生溢油，虽然油膜扩散面积较大，到移动路径基本沿着深水航道，所以后期扫海面积并没有涨潮时大。
| 潮况 Tide condition |
扩散面积 Diffusion area (km2) |
扫海面积 Sweeping area (km2) |
|||||||
| 6 h | 12 h | 18 h | 24 h | 6 h | 12 h | 18 h | 24 h | ||
| 涨潮 Flood tide |
2.88 | 3.92 | 6.87 | 9.64 | 8.58 | 22.39 | 62.47 | 70.46 | |
| 落潮 Ebb tide |
3.52 | 4.81 | 7.76 | 10.24 | 10.52 | 19.36 | 37.69 | 55.85 |
|
4.4 风对溢油漂移扩散的影响
有风情况下，油膜进入水体后，其漂移路径、扩散特征受潮流及风的共同影响。风作用对钦州湾溢油漂移扩散有显著的影响，不同风向下的油膜漂移方向、影响区域有明显的区别(图 3)。N风情况下，溢油将沿着钦州湾颈中部进入茅尾海，受西北潮流及北风的影响，溢油几乎没有触及七十二泾一带的岛屿；E风情况下，油膜受风的影响表现得尤为明显，东风推动下，油膜明显偏向西向移动，影响至龙门一带海域及岸线；S风情况下，溢油受南风和西北向流的影响，移动路径更偏向北向，使得七十二泾一带的岛屿受到明显的影响，同时，因为有南风的作用，进入茅尾海后，油膜移动更偏向茅尾海东部。W风作用下，油膜几乎紧贴钦州港附近岸线移动，进入七十二泾东部区域，受岸线影响明显，油膜移动距离很短，影响面积较小。
|
| 图 3 不同风向下涨潮时24 h溢油影响范围 Fig. 3 Influence scope of spilled oil affected by different wind direction at flood tide within 24 h |
可见，风对溢油的影响通常比较显著，一般5 m/s左右的风即可对油膜漂移扩散产生明显的影响。风对溢油的影响主要通过表面风生流的形式出现。在风力的作用下，表层水受风摩擦的影响产生一定的水流动，而这部分的风生流将带动漂浮于水面的薄油运动，从而对溢油结果产生显著影响。
5 结论本研究利用计算机代码，结合Windows可操作界面，构建广西近海海上溢油应急预报系统。系统自带了风漂移系数、溢油扩展系数、质量转移系数等经验参数，同时允许设定溢油时间、溢油量等自定义参数，简化计算的同时，保证了一定的计算精度。通过预先计算的水动力模型计算结果作为背景场，驱动溢油模型计算溢油漂移扩散过程。以钦州湾为例，涨潮时刻发生溢油时，溢油将向西北偏北方向移动，由钦州湾湾颈中部进入茅尾海；而落潮时刻发生溢油时，溢油主要沿着深水航道向钦州湾外湾方向移动。溢油在进入水体后，同时进行扩展及漂移过程。涨潮时刻发生溢油时，受钦州湾湾颈水流集中影响，油膜扩散面积相对较小，进入茅尾海后则明显增加。涨潮时发生溢油，扫海面积相对落潮时先小后大，主要原因在于涨潮时溢油开始阶段油膜面积较小，进入茅尾海后油膜面积增加较快，且进出茅尾海油膜移动路径少有重合。风对钦州湾溢油漂移扩散结果有显著的影响，不同风向下油膜漂移方向、影响区域有明显的区别。
系统计算结果基本能反映溢油进入水体后的漂移扩散特征。系统具有方便、准确以及良好的人机互交等特点，可为广西近海海上溢油事故防范、预警预测、应急处理、生态环境损害评估等提供科学依据，也可以为广西近海海上溢油应急预警预测提供技术支持。
| [1] |
BLOKKER P C.Spreading and evaporation of petroleum products on water[C]//Proceedings of the 4th International Harbour Conference.Antwerp, the Netherlands: [s.n.], 1964: 911-919.
|
| [2] |
FAY J A. The spread of oil slicks on a calm sea[M]. Oil on the Sea, New York: Plenum Press, 1969: 53-63.
|
| [3] |
MACKAY D, MATSUGU R S. Evaporation rate of liquid hydrocarbon spills on land and water[J]. Canadian Journal of Chemical Engineering, 1973, 51(4): 434-439. DOI:10.1002/cjce.5450510407 |
| [4] |
LIU D S K, LEENDERTSE J J.A 3-D oil spill model with and without ice cover[C]//Proceeding of Mechanics of Oil Slicks.Paris, France: Association Amical des Ingenieurs, 1981.
|
| [5] |
LEHR W J, CEKIRGE H M, FRAGA R J. Empirical studies of the spreading of oil spills[J]. Oil and Petrochemical Pollution, 1984, 2(1): 7-11. DOI:10.1016/S0143-7127(84)90637-9 |
| [6] |
JOHANSEN O.The Halten Bank experiment-observations and model studies of drift and fate of oil in the marine environment[C].Environment Canada. Proceedings, 11th Arctic Marine Oil Spill Program Technical Seminar.Ottawa: Environment Canada, 1984: 18-36.
|
| [7] |
杨红, 刘成秀, 李曰嵩, 等. 长江口南港水道溢油数值模拟[J]. 海洋通报, 2013, 32(3): 345-351. |
| [8] |
黄娟, 曹雅静, 高松, 等. 渤海海上溢油漂移扩散数值模拟研究[J]. 海洋科学, 2014, 38(11): 100-107. DOI:10.11759/hykx20130910002 |
| [9] |
李怀明, 娄安刚, 王璟, 等. 蓬莱19-3油田事故溢油数值模拟[J]. 海洋科学, 2014, 38(6): 70-77. |
| [10] |
严志宇, 熊德琪, 殷佩海. 海上溢油风化模型评述[J]. 大连海事大学学报, 2001(4): 36-39. DOI:10.3969/j.issn.1006-7736.2001.04.009 |
| [11] |
李涛, 李筠, 陈俊峰. 基于OILMAP模型的洋浦港溢油风险研究[J]. 交通节能与环保, 2012, 8(3): 61-67. DOI:10.3969/j.issn.1673-6478.2012.03.014 |
| [12] |
吴亚楠, 董胜, 高俊国, 等. 海庙港溢油扩散数值模拟研究[J]. 海洋湖沼通报, 2015(2): 177-184. |
| [13] |
刘钦政, 张存智, 刘煜, 等. 渤海溢油数值预报研究[J]. 海洋预报, 2005, 22(1): 70-76. |
| [14] |
宋泽坤, 程和琴, 刘昌兴, 等. 长江口溢油数值模拟及对水源地影响[J]. 长江流域资源与环境, 2013, 22(8): 1055-1063. |
| [15] |
尼建军, 王新怡, 张凤烨, 等. 基于FVCOM的渤海潮波数值模拟[J]. 海洋科学, 2013, 37(2): 89-94. |
| [16] |
杨万康, 尹宝树, 杨德周, 等. 基于FVCOM的南海北部海域潮汐潮流数值模拟[J]. 海洋科学, 2013(9): 10-19. |
| [17] |
王翠, 郭洲华, 李青生, 等. 基于MIKE SA模型的厦门西港海域溢油影响的数值模拟研究[J]. 应用海洋学学报, 2014, 33(2): 229-235. DOI:10.3969/J.ISSN.2095-4972.2014.02.011 |
| [18] |
姜卫星.黄浦江溢油事故的数值模拟研究[D].上海: 同济大学, 2007.
|
| [19] |
吴晓丹, 宋金明, 李学刚, 等. 海洋溢油油膜厚度影响因素理论模型的构建[J]. 海洋科学, 2010, 34(2): 68-74. |