瓯江口航道海域溢油扩散数值模拟

黄毅峰 ，许 婷 ，刘 涛

（1.长沙理工大学水利学院，长沙410076；2.交通运输部天津水运工程科学研究所工程泥沙交通行业重点实验室，天津300456）

近半个世纪以来全球发生了近10起重大的原油泄露事件。这些原油的泄露造成水质污染，海域环境恶化，危及海洋生物，造成生态破坏。人民生活也因此受到影响，企业因此蒙受巨大的损失。今年4月份美国墨西哥湾钻井平台发生爆炸引发漏油事件引起了国际社会的广泛关注。

鉴于溢油事件的频发性和严重性，国内外有很多学者对溢油事件进行了数值模拟研究。国外的比较著名的有美国的OILMAP系统、英国的OSIS系统、挪威的OSCAR系统、荷兰的MS4系统、比利时的MU-SLICK系统［1］。国内学者在不同的海区都进行了溢油数值模拟研究，俞济清、黄立文等在舟山水域建立了溢油预报模型［2-3］；珠江水域、渤海湾水域、辽东湾水域、胶州湾等水域建立了溢油预报模型［4-7］。文中利用MIKE21 HD模块建立了瓯江口二维水动力模型，并用近几年最新实测资料验证模型的可靠性和适用性，模拟结果为瓯江口航道二维溢油模型的建立提供水动力基础数据，再次利用MIKE21 SA模块建立瓯江口航道二维溢油模型，应用“油粒子”模型模拟输移、风化和热量迁移等过程，对瓯江口航道一期工程溢油泄露事故进行影响预测，选取7种具有代表性情形的事故场景，对不同气象条件和不同事故发生地点的工况进行了模拟分析。瓯江口以外的海域存在着水产养殖、风景旅游、生态保护等敏感区域，一旦发生溢油事故，将对周围海洋环境产生严重的影响，应采取溢油应急措施，控制污染影响程度，同时追踪油膜污染带，并通知港航监督部门。

1 工程海域自然动力条件

1.1 潮汐类型

瓯江口是浙江省第二大河流，为强潮河口，据2005年6月水文全潮测验统计，各潮位站（HK1+H01）/HM2=0.24～0.28，均小于0.5，因此潮汐属正规半日潮类型。一昼夜两个潮，潮高不等现象较为明显［8］。

1.2 潮流特征

据1999年10月～2005年7月瓯江河口及其附近海域水文全潮各站潮流速观测资料调和分析和计算，各站主要全日分潮流与主要半日分潮流的比值F=（W01+WK1）/WM2均小于0.5，各站浅水分潮流比值G==0.27～0.48，综合该两项指标，表明该水域潮流属正规浅海半日潮流类型；由于涨落潮流受水域地形限制，基本呈往复流运动。

1.3 风况

本地区全年常风向为N-NE向，频率为54.7%，全年平均风速为3.8 m/s；强风向为SSW向，最大风速为32.0 m/s（1975年8月12日）。全年中夏季多为SW向大风，春秋季节多为偏S向或偏N向大风，以偏N向大风为主，冬季盛行N-NE向大风。根据其多年风资料统计，多年平均≥6级风日数为37 d，多年平均≥7级风日数为8.5 d。

1.4 波浪特征

据洞头岛NE端的甲米礁东岸波浪观测站（1990年12月～1992年9月）资料统计，该海域常浪向为NNE向，频率为55%，次常浪向为S-SW向，频率为20%；涌浪向为ENE向，频率为26.3%，风浪常浪向主要为N-NE向，频率55%，次为S-SSW向，频率为16.5%。最强浪向为NNE-ENE向，年平均H4%波高大于0.6 m，对应年平均周期为0.5～4.2 s，多发生在8、9月台风侵袭时期，实测最大波高4.3 m（E向），由9216号台风造成。

据温州浅滩围涂工程1#（北堤处）和2#（温州浅滩南口水道附近）波浪观测站资料统计，2001～2004年1#站测得波高2.1～3.0 m出现24次（2001年12月～2002年5月14次、2002年6～7月10次），波高大于3.0 m仅出现2次（2001年12月）；2001～2003年2#站波高2.1～3.0 m出现3次（2001年6月2次、2002年9月1次），两站波高为0.0～0.5 m的出现年频率分别为95.0%和91.9%。

2 溢油预测模型概述

丹麦水环境研究所开发的MIKE 21/3溢油分析（SA）模块，是基于欧拉-拉格朗日理论体系，通过对油膜在水体中的扩展、传输（水流和风场作用）、紊动扩散、分散（夹带）、蒸发、乳化、溶解等各种过程的模拟，MIKE 21/3 SA能提供油膜随时间变化的漂移位置、厚度，以及漂移过程中粘度、油膜表面温度、倾点等属性的变化。另外，MIKE 21/3对复杂水体（如各种水工建筑物）的水动力能进行精确模拟，强大的前后处理功能适用于决策系统的开发。MIKE系统是国际上同类型软件中最先进的模拟工具之一，在世界各地已有许多成功的应用案例［9］。

2.1 输移过程

油粒子的输移包括扩展、漂移、扩散等过程［9］，这些过程是油粒子位置发生变化的主要原因，而油粒子的组分在这些过程中不发生变化。

（1）扩展运动。采用修正的Fay重力-粘力公式计算油膜扩展

式中：Aoil为油膜面积，Aoil=πR2oil；Roil为油膜直径；Ka为系数；t为时间；油膜体积为

初始油膜厚度

（2）漂移运动。油粒子漂移的作用力是水流和风拽力，油粒子总漂移速度由以下权重公式计算

式中：Uw为水面以上10 m处的风速；Us为表面流速；cw为风漂移系数，一般在0.03～0.04。风场数据从气象部门获得，而流场从二维水动力模型计算结果获得。

（3）紊动扩散。假定水平扩散各向同性，一个时间步长内α方向上的可能扩散距离Sα可表示为

式中：［ ］

R1-1为-1～1的随机数；Dα为α方向上的扩散系数。

2.2 风化过程

油粒子的风化包括蒸发、溶解和形成乳化物等过程，在这些过程中油粒子的组成发生改变，但油粒子水平位置没有变化。

（1）蒸发。油膜蒸发受油分、气温和水温、溢油面积、风速、太阳辐射和油膜厚度等因素的影响。假定在油膜内部扩散不受限制（气温高于0℃以及油膜厚度低于5～10 cm时基本如此），油膜完全混合，油组分在大气中的分压与蒸气压相比可忽略不计。

蒸发率可由下式表示

式中：N为蒸发率；kei为物质输移系数；为蒸汽压；R为气体常数；T为温度；M为分子量；ρ为油组分的密度；i为各种油组分。kei由下式估算

式中：k为蒸发系数；Sci为组分i的蒸气Schmidts数。

（2）乳化。油向水体中的运动机理包括溶解、扩散、沉淀等。扩散是溢油发生后最初几星期内最重要的过程。扩散是一种机械过程，水流的紊动能量将油膜撕裂成油滴，形成水包油的乳化。这些乳化物可以被表面活性剂稳定，防止油滴返回到油膜。在恶劣天气状况下最主要的扩散作用力是波浪破碎，而在平静的天气状况下最主要的扩散作用力是油膜的伸展压缩运动。

从油膜扩散到水体中的油分损失量计算

式中：Da为进入到水体的分量；Db为进入到水体后没有返回的分量

油滴返回油膜的速率为

油中含水率变化可由以下平衡方程表示

式中：R1和R2分别为水的吸收速率和释出速率。

（3）溶解。溶解率用下式表示

2.3 热量迁移

蒸气压与粘度受温度影响，而且观察发现通常油膜的温度要高于周围的大气和水体。图1为油膜的热平衡示意图。

图1中1为大气与油膜之间的传热过程；2为大气与油膜之间热辐射过程；3为太阳辐射；4为蒸发热损失；5为油膜与水体之间的热量迁移；6为油膜与水体之间散发和接受的热辐射。

3 水动力条件模拟

溢油扩散模拟结果是建立在水动力条件模拟的基础之上，文中采用MIKE21 HD模块计算水动力条件。

3.1 平面二维水流控制方程

连续方程

运动方程

式中：x、y为直角坐标系坐标；t为时间；h为水深（基准面到床面的距离）；ζ为潮位（基准面到自由水面的距离）；u、v分别为x、y方向的垂线平均流速分量；f为科氏系数；g为重力加速度；Ex、Ey分别为x、y方向的水平紊动粘性系数；τx、τy分别为波流共同作用下床面剪切力在x、y方向的分量；Sxx、Sxy、Syy分别为各方向的波浪辐射应力。

3.2 计算域的确定及网格剖分

该海域岛屿众多，潮波运动受地形控制显著、径流影响明显的特点，要复演该水域的潮流运动，必须模拟瓯江口整体海域的潮波传播过程，建立足够大的计算区域。因此，此次数学模型的计算区域西边界取至120°30′经度线，东边界至121°45′经度线，北至乐清湾，南至南麂岛南侧27°10′纬度线。东西向宽约120 km，南北向为150 km。模型采用正四边形网格，计算网格步长为Dx=Dy=50 m，并就复杂的陆地和岛屿岸线进行二次精细光滑处理，以保证计算精度。

3.3 模型验证

为了验证模型的合理性，采用2006年10月瓯江口水文全潮测验（每次测验均为大、中、小3个潮型）中的大潮资料进行了验证。水文全潮共布置了13个潮位站、14条水文垂线（图2）。其中1#～7# 垂线施测大、中、小潮，8#～14# 垂线仅施测大潮，测验内容为潮位、流速、流向等。该次水文全潮规模宏大、内容丰富、测点涵盖了瓯江口大范围海域。部分测站大潮潮位及流速流向验证曲线见图3～图4。

通过潮位、流速、流向的验证，模型中无论计算的位相、量值还是过程线，均与原型吻合良好，符合《海岸与河口潮流泥沙模拟技术规程》［10］的要求，可以用来进行溢油扩散数值模拟计算研究。

3.4 工程海域潮流场模拟

采用已验证的潮流数学模型对工程海域潮流场进行模拟计算。从工程海域的流场看，瓯江口及口外的海域涨落潮流主要受岸线和地形控制及工程方案布置的影响，总体呈现往复流运动；由于江口外的岛屿众多，导致附近区域的流态比较复杂。涨潮的时候潮流从外海区域向瓯江口门涌入，受小门岛、大门岛和洞头列岛的状元岙岛和霓屿岛影响分成两股主要的涨潮流，其中大门岛的南北两侧形成了瓯江口潮流的主流和口门的分流区，涨潮流在流经小门岛被分为一南一北两支，北支流在流向沿岸时受岸线约束又分为南北2个流向，向南涨潮流跟大门岛和小门岛之间的支流汇合成口门分流。落潮时基本呈现出与涨潮相反的运动过程。当落潮时，外海潮位降低，大小门岛的西侧及南北侧的落潮流向外海运动，落潮流经小门岛分为南北两支流，汇合小门岛东侧海域，而后跟大门岛南侧的主流汇合向外海方向运动，至低平潮沿岸形成一定范围的露滩。工程海域涨、落急流场分布见图5。

4 溢油影响预测

在水动力条件模拟的基础上，利用MIKE21 SA模块建立瓯江口航道二维溢油模型，应用“油粒子”模型模拟输移、风化和热量迁移等过程，对瓯江口航道一期工程溢油泄露事故进行影响预测，预测时刻选择大潮的4个潮时（高潮时、落急时、低潮时、涨急时），考虑瓯江口地区全年常风向为N－NE向，频率为54.7%，全年平均风速为3.8 m/s；强风向为SSW向，最大风速为32.0 m/s。溢油计算的结果受水动力、风向、溢油事故发生位置、发生时间等影响。综合考虑以上各影响因素，共计36个计算工况。限于篇幅，这里选取7种代表工况进行结果展示，其扫海范围及漂移轨迹见图6～图12，图中所示为时间总长2个潮周期（即24 h），时间分别取2 h、6 h、12 h和24 h的结果。不同时间段过后其溢油扫海面积见表1。

表1 事故溢油扫海面积预测结果Tab.1 Sweeping area of oil spill accident prediction

从图6～图12和表1中均可以看出，溢油事故发生后其扫海范围及迁移路径与溢油发生地点、发生时间、风向情况密切相关，不同情况下，其扫海范围及迁移路径差别较大。根据施工期和营运期溢油风险预测可得：（1）不同风向的影响时。风对溢油油膜漂移轨迹的影响比较大。静风条件下，油膜的漂移路径主要受控制于潮流动力，在NE向风的作用下，油膜会发生明显向西北移动现象，而在SSW风向下，油膜同样会发生明显向东北至东南该区间方向偏转；（2）不同潮时的影响。溢油的漂移轨迹受溢油时刻的影响极大，在潮周期的不同阶段发生溢油事故时，溢油油膜漂移轨迹可能完全不同。可以通过扫海范围及油膜漂移轨迹图明显看出不同潮时所带来的不同溢油事故后果。若发生溢油事故的地点不同，不同潮时所带来的影响差别大小也不尽相同。

5 结论

利用MIKE21 HD模块，建立瓯江口水域水动力数学模型模拟潮流特性。在此基础上利用MIKE 21/3溢油分析（SA）模块，对瓯江口航道二期工程溢油泄露事故进行影响预测。计算结果表明：（1）不同风向的影响。风对溢油油膜漂移轨迹的影响比较大。静风条件下，油膜的漂移路径主要受控制于潮流动力，在NE向风的作用下，油膜会发生明显向西北移动现象，而在SSW风向下，油膜同样会发生明显向东北至东南区间方向偏转。（2）不同潮时的影响。溢油的漂移轨迹受溢油时刻的影响极大，在潮周期的不同阶段发生溢油事故时，溢油油膜漂移轨迹可能完全不同。可以通过扫海范围及油膜漂移轨迹图明显看出不同潮时所带来的不同溢油事故后果。若发生溢油事故的地点不同，不同潮时所带来的影响差别大小也不尽相同。

［1］徐艳东.海上溢油风化过程及其预测模型研究［D］.青岛：中国海洋大学，2006.

［2］俞济清，魏敏，黄立文，等.中国舟山港溢油模拟信息系统研制［J］.武汉理工大学学报，2002，26（4）：463-465.YU J Q，WEI M，HUANG L W，et al.The spilled oil simulation information system designed for Zhoushan harbor［J］.Journal of Wuhan University of Technology，2002，26（4）：463-465.

［3］黄立文，宋文铜，孙俊，等.实时流场预报及在海面溢油轨迹预测中的应用［J］.中国航海，2000（2）：36-41.HUANG L W，SONG W T，SUN J，et al.Real time current field forecast and its application to marine oil-spill trajectory prediction［J］.Navigation of China，2000（2）：36-41.

［4］熊德琪，杨建立，严世强.珠江口区域海上溢油应急预报信息系统的开发研究［J］.海洋环境科学，2005（2）：63-66.XIONG D Q，YANG J L，YAN S Q.R&D on oil spill contingency forecast information system for Pearl River estuary［J］.Marine Environmental Science，2005（2）：63-66.

［5］张波，吴冠.渤海溢油应急预报系统［J］.海洋环境科学，1998（3）：72-75.ZHANG B，WU G.Contingent predicting system of oil spill for Bohai Sea［J］.Marine Environment Science，1998（3）：72-75.

［6］窦振兴，宫成，张存智，等.海上溢油动态预测应急系统［J］.中国海洋平台，1995（5）：31-33.DOU Z X，GONG C，ZHANG C Z，et al.The Dynamic Predication and Emergency System of Oil Spill at Sea［J］.China Offshore Platform，1995（5）：31-33.

［7］娄安刚，王学昌，孙长青，等.胶州湾海面溢油轨迹的数值模拟［J］.黄渤海海洋，2001（1）：1-8.LOU A G，WANG X C，SUN C Q，et al.A simulation of Oil spill Trajectory on Sea Surface of Jiaozhou Bay［J］.Journal of Oceanograpgy of Huanghai&Bohai Seas，2001（1）：1-8.

［8］赵洪波，许婷.瓯江口航道二期治理潜堤工程潮流泥沙盐度数学模型研究［R］.天津：交通运输部天津水运工程科学研究所，2009.

［9］姜卫星.黄浦江溢油事故的数值模拟研究［D］.上海：同济大学，2007.

［10］JTJ/T233-98，海岸与河口潮流泥沙模拟技术规程［S］.