湛江湾海域溢油漂移扩散数值模拟研究

孙 琰，刘大召,3，许源兴，李 卓，韩泽文

(1.广东海洋大学电子与信息工程学院，广东 湛江 524088；2.中海油信息科技有限公司天津分公司，天津 300450；3.广东省海洋遥感与信息技术工程技术中心，广东 湛江 524088)

引 言

近年来随着我国经济不断发展，对能源的需求日益增长，我国海域发生海洋溢油事故的风险显著增加[1]。海上溢油事故不仅对近海海洋生态环境造成严重破坏，而且会直接危害我国近海社会经济正常发展[2]。

我国是世界最大的原油进口国，湛江港作为中国大陆通往东南亚、欧美、非洲等国家和地区航程最短的天然深水良港与全国最大的陆岸原油专业化码头之一，近年来湛江湾海域原油储运船舶数量上升使得该海域发生溢油事故的风险不断加大[3]，此外湛江湾内生态保护区与养殖区分布广泛，溢油事故会对湛江湾内生态、经济产生严重破坏，因此有必要建立该海域的溢油数学模型进行事故风险模拟与环境敏感区影响分析。对海上溢油的数值模拟预报研究开始于20世纪60年代[4]，经历了油膜扩展理论、对流扩散方法与油粒子法三个阶段，其中由Johansen[5]与EIliot[6]提出的拉格朗日油粒子追踪法目前在溢油模型中得到广泛的应用。国内学者采用数值模拟方式在长江口[7]、胶州湾[8]、渤海湾[9]、江苏滨海[10]、大亚湾[11]、北部湾[12]、福建东吾洋湾[13]、福宁湾[14]等典型区域做了大量关于溢油事故的预报研究工作，表明利用数学模型的方法对溢油事故模拟高效准确，满足防范和治理突发海洋污染事故的要求。

本研究采用Delft3D-FLOW水动力模块建立湛江湾潮流数学模型，进而与Delft3D-PART溢油模块耦合对湛江湾内溢油漂移扩散进行数值模拟，预测分析溢油的漂移扩散、扫海面积、到达环境敏感区的时间、岸线吸附，为海洋环境保护与溢油事故应急处理提供科学依据。

1 材料与方法

1.1 Delft3D-FLOW水动力模型控制方程

Delft3D-FLOW是一个可用于河口、港湾以及海洋的水动力数值模拟模块，对于流场的模拟已经被应用于许多研究之中。采用Delft3D-FLOW模块对湛江湾海域进行水动力环境模拟试验，经过实测潮流数据验证后输入溢油模块作为溢油水动力的驱动条件。Delft3D-FLOW水动力模块在Navier-Stokes方程的基础上，基于静水压强假定和Boussinesq假定，采用循环隐式进程ADI(Alternating Direction Implicit)法离散求解非线性浅水方程。在正交曲线坐标系(ξ,η)下，沿水深积分的连续性方程为：

(1)

模型在ξ和η方向的动量方程为：

(2)

(3)

式中，f为科氏力；ρ0为水的密度；Pξ和Pη为ξ和η方向水压力梯度；Fξ和Fη为ξ和η方向紊动量通量；Mξ和Mη为源(汇)项的动量分量。

1.2 Delft3D-PART溢油控制方程

Delft3D-PART溢油模块耦合FLOW模块的计算结果，通过拉格朗日粒子追踪法，将溢油分成若干个代表一定质量的油粒子来模拟油膜在水中的输移扩散，此外模型也将油粒子的蒸发、乳化、岸线吸附等变化过程也纳入计算。

1.2.1 油的溢出

根据Fay[15]的公式计算溢油泄漏半径：

(4)

式中，k1和k2为Fay常数；V0为溢油体积；ρw为水体密度；ρ0为油的密度；g为重力加速度；VW为水体的运动粘滞系数。

1.2.2 漂移

漂移又称溢油平流，是油膜在受到风、表层平流的作用下产生的移动，是溢油运动转移的主要因素之一。溢油由风引起漂移过程，受当地风向风速以及偏转因子与偏转角的影响。漂移速度通常在当地风速的 2.5～4.4%，风对溢油的漂移影响关系为：

V=Cwd(Vw-Vf)

(5)

式中，Cwd为风拖曳力系数；Vw为风速；Vf为流速。

1.2.3 蒸发

蒸发是较轻的石油烃组分从液态转为气态向大气进行质量传输的过程[16]。Delft3D-PART溢油模块通过将溢油挥发组分的蒸发设为一阶衰减过程来实现蒸发过程，相应的计算公式如下：

(6)

式中，Fvol为挥发组分；Fv为溢油蒸发系数；K为蒸发速率常数，当Fvol-Fv<0时，无蒸发。

1.2.4 乳化

乳化是溢油形成油包水和水包油的平衡过程，该过程是一个不可逆转的过程，其将液体转化为高粘度的半固体物质[17]。乳化过程本身是比较快速的，实验室0.1—3 h内就能完成乳化过程。Delft3D-PART模块根据Mackay等[18]提出的算法计算乳化过程，相应计算公式如下：

(7)

式中，F′wc为吸水率；Uw为风速；Fwc为含水量；C1为乳化率取值 2×10-6；C2为最大含水量。

1.2.5 分散(夹带)

浮油在水中的分散或夹带是通过零阶衰减来实现的，分散的速率与浮油浓度无关，只取决于耗散的波能与油品。分散速率相应的计算公式由Delvigne[19-20]给出：

(8)

(9)

式中，Q为分散速率；Q(d)是直径为d的液滴单位直径的分散速率；C′′为夹带常数，取决于油的种类；D为单位面积耗散的波能；Fwc为每个波周期破波数；N(d)为油粒度分布函数。

1.2.6 岸线吸附

油粒子在漂移到岸边时会受岸线吸附作用不再进入水体，模型计算时对每个油粒子设定一个0～1之间随机数，当油粒子靠岸且随机数小于模型中给定的岸线吸附概率时，该油粒子将粘附在岸上或水底不再进入计算，湛江湾内岸线吸附概率选取为25%[21]。

1.3 水动力模型网格与参数

基于Delft3D-FLOW模块建立水动力模型为Delft3D-PART溢油模块提供水动力基础。模型计算域包括湛江湾与其近海水域约5.65×103km2，采用正交曲线网格划分并在湛江湾内重点关注水域进行网格加密处理，整个计算域包含10246 个网格单元与6406 个网格节点，网格最大分辨率达到30 m。水深数据采用中国人民解放军海军海道测量局高精度电子海图，外海开边界采用俄勒冈州立大学中国海潮汐模型所提供的八个主要分潮(M2，S2，N2，K2，K1，O1，P1，Q1)调和常数。模型初始水位、流速设为0，计算时间步长为60 s。模型的计算域网格划分以及水下地形如图1所示。

图1 模型计算域网格和地形

1.4 溢油模块参数及工况设置

作为中国十大港口之一的湛江港近年船舶密度逐年增大，其40 万t级航道所在海域海况复杂，航线多且繁忙，另湛江湾航道深水区涨落潮流速大，船舶在航道航行过程中易受横流作用偏离航道中心线，从而有船舶碰撞的风险。本文模拟的溢油点选在船舶往来频繁且船速较快的航道上，溢油点(21.09°N，110.44°E)的位置如图2所示。中石化、中石油、中国航油集团等在湛江建有中大型原油库，湛江湾海域原油运输船舶数量多，故本次溢油模拟泄漏油品为密度890 kg/m3的原油，泄漏量为100 t，其中原油最大含水率为0.75，蒸发系数为0.1。在PART模块中将溢油分为40000个油粒子，风的拖曳力系数设为3 %[22]。

风对油膜的转移、蒸发、乳化等过程都有影响，湛江湾海域地属季风区，风向与风速随季节变化明显，根据湛江市多年气象资料与环境敏感区位置(图2)选取五种溢油发生时的风况条件：静风；夏季常风；冬季常风；不利风况1；不利风况2，另选取两种潮时(涨潮、落潮)作为泄漏发生时的潮流条件。以上两种潮时，五种风况组合为10 种典型工况，具体见表1。

图2 溢油点、实测站点与环境敏感区位置

表1 溢油模拟情形组合

2 结果与讨论

2.1 模型验证

为保证Delft3D-PART模块对溢油漂移扩散模拟的准确性，对Delft3D-FLOW模块水动力结果进行潮位、潮流验证，验证采用在一个全潮期间(2016年3月8日至3月16日)研究区内L1—L3三个潮位站点，S1—S4四个潮流站点实测数据(点位具体分布见图2)与模型计算值相比较，篇幅限制本文具体仅列出L1、L2大、小潮的潮位验证(图3)与大潮期间(2016年3月8日12时至3月9日12时)S1、S2潮流验证(图4)结果。通过比较得出大小潮期间最低潮位误差小于10 cm，误差小于10%；流速、流向计算值与实测值误差在±12%以内，能较好模拟出潮流运动特征，综上所建立潮流数学模型能较准确地反映湛江湾海域潮流环境情况，可以为耦合溢油模块提供水动力基础。

图3 (a) L1大潮潮位验证 (b) L2大潮潮位验证 (c) L1小潮潮位验证 (d) L2小潮潮位验证

图4 大潮流速流向验证图(S1、S2)

2.2 水动力模拟结果分析

图5为模型计算域的涨急流场与落急流场的分布图，通过计算分析得出湛江湾及其近海海域属于不规则半日潮，涨、落潮方向大致与东海岛、南三岛东侧岸线垂直，涨潮时来自外海的潮波通过模型开边界进入计算域，海水大致从东南向西北方向传播，落潮时潮流沿原路返回，流向外海。

图5 模型计算域涨急(左图)与落急(右图)流场分布

图6为湛江湾内涨潮与落潮时刻的流场分布图，湛江湾内海域海流受潮汐作用较为明显，在湛江湾独特狭长地形影响下呈现明显往复流特征，涨、落潮期间海水从湛江湾口航道汇入或流出，涨潮时海流从湛江湾口以及南三岛北侧水道汇入湾内，海流在湾口到湾顶大致由东南向西北流，落潮时湛江湾内北部水道与西侧海流在东头山岛北侧汇集后一并向东南流向湾外；流速方面，湾口区域为流速高值区，其流速最大可达1.6 m/s，海流进入湾内后，过水断面增大，使得流速下降，但基本为航道深水区平均流速大于其它海域，浅滩、岛屿以及近岸地区流速较小，涨潮平均流速略大于落潮平均流速。

图6 湛江湾内涨急流场(左图)与落急(右图)流场分布

2.3 溢油扩散影响预测

2.3.1 油膜24 h漂移扩散

由于所模拟溢油事故是一个短时间内发生溢油泄漏情形，为对突发溢油事故油膜高值区快速定位追踪，对十种工况在24 h内的计算结果进行后处理，根据《中国海上船舶溢油应急计划》中影响评价指标，以轻度影响所对应的油膜厚度(1 μm)为起算值进行统计，绘制油膜漂移扩散分布图，能直观地反应不同工况下溢油的运动轨迹与影响范围。图7、图8为十种典型溢油工况情形下，油膜高值区在24 h内漂移扩散分布图(油膜厚度单位：μm，油膜厚度绘制时间间隔：30 min)。

图7 涨潮时发生溢油24 h内油膜漂移扩散分布

图8 落潮时发生溢油24 h内油膜漂移扩散分布

通过图7、图8油膜漂移扩散分布图可以得出不同风况、潮流条件都对油膜的漂移扩散产生影响。

在涨潮时刻发生溢油，各种风况下油膜都会随着涨潮流向西北方向的湾内漂移，这是由于事故刚发生后，油膜厚度集中、质量较大且溢油事故发生在流速较大的航道深水区，此时油膜漂扩散方向由潮流方向主导。涨潮-静风条件下(工况1)，油膜漂移扩散只受潮流作用影响，油粒子随着涨、落潮在湾内振荡，岸线吸附、夹带、乳化作用较小，油粒子随着涨、落潮在湾内振荡，此种情况下厚度大于1 μm的油膜分布影响范围最广。涨潮-夏季常风条件下(工况2)，涨潮流方向与夏季常风风向相近，使得油膜3 h内到达湛江近岸，油膜在湛江近岸海域受涨、落潮与夏季常风(SE风)共同作用沿岸线方向振荡，油粒子在此岸线受到大量吸附，此工况下溢油事故只对海滨公园保护区产生污染影响。涨潮-冬季常风条件下(工况3)，事故发生后油膜向西北方向漂移扩散至东头山岛西北侧，4 h后跟随落潮流沿东头山岛岸线向东南漂移，在此期间部分油粒子被东头山岛岸滩吸附，油膜漂移过东头山岛后在冬季常风(N风)作用下继续向南偏移漂移扩散至东海岛北岸工业与城镇用海区。在第一种不利风(SW风)影响下(工况4)，溢油事故发生后2.5 h油膜漂移扩散至特呈岛海洋保护区，由分布图可知此时油膜厚度较大，部分油粒子被吸附在特呈岛南岸，4 h后落潮使得潮流方向发生变化，在SW风与落潮流共同作用下油膜向东经过南三岛养殖区后在第6 h漂移至南三岛南侧岸滩，对该处红树林保护区产生影响。在第二种不利风(S风)作用下(工况5)，油粒子向上游水道漂移扩散，发生溢油后2.5 h油粒子在特呈岛海洋保护区西岸部分被吸附，5 h内潮流方向与风向相同，两驱动因素叠加使油膜转移比静风时速度更快且对湛江海滨公园水域造成污染。

在落潮时发生溢油，各风况油粒子在事故发生2 h内都会随着落潮流大致向东南方向漂移扩散。落潮-静风条件下(工况6)，油膜在湾内随往复流振荡，最远向东到达湛江湾口海域。落潮-夏季常风条件下(工况7)油膜先随落潮流向东南漂移至东海岛北侧工业城镇用海区界线北侧，随后在涨潮流与夏季常风(SE风)共同作用下，油膜漂移扩散至湛江近岸，在此油粒子受到岸线吸附。在落潮-冬季常风条件下(工况8)下油膜的漂移相较涨潮时刻溢油，溢油不再对东头山岛西北岸线产生影响，对东海岛工业城镇用海区东侧岸线产生了污染影响。两种不利风况(工况9、工况10)条件下，油粒子在不利风1(SW风)与不利风2(S风)作用下向东漂移扩散后向北到达南三岛近岸，对南三岛红树林保护区产生污染影响，而不利风2(S风)工况下，油膜对红树林保护区与南三岛养殖区都产生了不同程度的污染。

对各工况下溢油到达敏感区的时间做出统计如表2所示，可以得出溢油污染对南三岛红树林保护区、南三岛渔业养殖区、东海岛工业与城镇用海区、特呈岛海洋保护区与湛江海滨公园等主要环境敏感区有污染影响，油膜最快到达各环境敏感区的时间分别为5，6，1.5，2.5，4 h，其中只有不利风情况下会对南三岛红树林保护区与南三岛渔业养殖区产生影响，特呈岛海洋保护区由于距离溢油点较近最易受到污染威胁，有6种工况下的溢油事故对其产生污染影响。

表2 事故发生24 h内油膜到达敏感区时间

2.3.2 溢油扫海面积统计

溢油事故发生24 h内各工况条件下油膜累积扫海面积变化统计如图9，图9(a)为涨潮时发生溢油的五种不同风况累积扫海面积统计，图9(b)为落潮时发生溢油的五种不同风况累积扫海面积统计，分析得出溢油发生后油膜的扫海面积在一定时间内随时间增长而扩大，且在24h内大于1um厚度油膜扫海面积达到最大。在涨潮时发生溢油，24 h后油膜扫海面积：静风>不利风2>冬季常风>不利风1>夏季常风；在落潮时发生溢油，24 h后油膜扫海面积：静风>夏季常风>不利风2>不利风1>冬季常风。十种工况中涨潮-静风工况的扫海面积最大，24 h后面积达到57.89 km2，涨潮夏季常风工况扫海面积最小为9.06 km2。

图9 各工况下油膜累计扫海面积(a:工况1～工况5；b:工况6～工况10)

溢油事故在同一潮时不同风况下的油膜扫海面积变化不同(如图9a)，这表明风况对油膜的漂移扩散有着重要影响；同一风况条件下不同潮时发生事故也对油膜的扫海面积产生不同影响，如在夏季常风条件下，落潮时发生溢油的累计扫海面积约为涨潮时发生溢油的六倍，涨潮与落潮时发生事故的油膜扫海面积变化趋势存在较大差异，是由于事故发生后涨潮流带动油粒子加速到达岸边，大量油粒子受岸线吸附作用不再转移，最终溢油污染影响范围变小，表明在湛江湾海域潮流场对油膜的转移也存在一定影响。

2.3.3 岸线吸附分布

岸线吸附是溢油粘附于岸边不再进入水体的过程，粘附于岸滩的溢油对海岸生态环境也会产生破坏，为追溯溢油被岸线吸附部分，绘制涨、落潮时发生事故72 h后被岸线吸附的溢油分布图(图10、图11)，涨潮时溢油影响岸线长度：静风>冬常风>夏常风>不利风2>不利风1，落潮时溢油影响岸线长度：静风>冬常风>夏常风>不利风1>不利风2。产生吸附的岸线与风向有明显相关性，静风时油膜在只受潮流作用下在湾内振荡，影响岸线范围最大，在有风情况下有油膜吸附的岸线相对较少，其中湛江湾顶岸线最容易受到溢油污染，有五种工况下的溢油会在此受到吸附。通过岸线吸附分布图可进一步明确溢油归宿，对溢油事故的应急处理与防范有指导意义。

图10 涨潮时发生溢油72 h后岸线吸附分布 (单位：kg/m2)

图11 落潮时发生溢油72 h后岸线吸附分布 (单位：kg/m2)

3 结论

基于Delft3D-FLOW模块建立湛江湾海域水动力模型，较好地反映了湛江湾海域流场特征，在此基础上采用Delft3D-PART模块耦合水动力结果模拟湛江湾海域航道发生100 t原油泄漏的事故，通过数值模拟结果分析得到以下主要结论:

(1) 湛江湾海域溢油的漂移受多种因素的共同影响，发生溢油时的不同潮时、不同风况都对溢油漂移轨迹、油膜扫海面积以及到环境达敏感区时间产生影响。静风时油膜漂移只受涨、落潮流影响，有风时受潮流与不同风况共同作用，其中油膜的长期扩散方向与风向一致。

(2) 绘制十种不同工况下油膜24 漂移扩散分布图，直观显示对各区域污染影响的油膜厚度。分析得出溢油对南三岛红树林保护区、南三岛渔业养殖区、东海岛工业与城镇用海区、特呈岛海洋保护区与湛江海滨公园等主要环境敏感区产生污染影响，油膜最快到达各敏感区的时间分别为5 h、6 h、1.5 h、2.5 h、4 h，其中南三岛红树林保护区与渔业养殖区只有在不利风况下会被溢油污染影响，特呈岛海洋保护区最容易受污染。

(3) 对溢油扫海累积面积进行统计，静风条件下油膜只受潮流作用影响沿涨落潮方向在湛江湾内振荡，岸线与岛屿的吸附、夹带、乳化等作用减小，使得静风条件下溢油影响范围更广，其中涨潮静风时溢油扫海面积最大，事故发生后24 h油膜累积扫海面积达到57.89 km2。

(4) 追溯溢油岸线吸附部分，绘制溢油岸滩吸附部分分布图，得出静风状况下产生吸附的岸线最多，且岸线吸附受与风况呈明显相关性，为岸线保护与事故防范提供依据。

综上，本次模拟对湛江湾内航道上潜在溢油事故风险的应急处理有积极意义，为海洋环境影响评价与制定相应的防护措施提供科学依据。