渤海海上溢油漂移扩散数值模拟研究

黄 娟, 曹雅静, 高 松, 徐江玲, 刘桂艳

(国家海洋局 北海预报中心, 山东 青岛 266061)

近年来, 国内外海上溢油事件频发, 造成海洋生态环境受到污染损害。2006年2月份, 山东长岛、龙口、蓬莱岸段及其附近海域发生大面积原油污染,当地水产养殖业和渔业损失惨重。2004年12月7日,巴拿马籍集装箱船和德国籍集装箱船在珠江口发生碰撞, 约1 200 t燃油溢漏, 海上形成了长9海里、宽200 m的油带, 造成我国近年来较大的一次海洋污染事故。2002年11月23日, 马耳他籍油轮“塔斯曼海”与中国籍“顺凯1号”轮相撞, 约200 t原油泄漏, 溢油污染影响渤海湾海域面积达359.6 km2。2010年4月 20日, 位于墨西哥湾的“深水地平线”钻井平台发生爆炸并引发大火, 大量漏油进入墨西哥湾, 5月原油漂浮带长200 km, 宽100 km, 且还在进一步扩散,此次漏油事件演变成美国历来最严重的油污大灾难,造成了巨大的环境和经济损失。2010年 7月, 在中国大连新港发生石油储备库输油管道发生爆炸事件,部分石油泄漏入海。这次溢油事故在我国历史上溢油量较大, 影响范围较广。2011年6月, 渤海中部蓬莱 19-3平台发生“井涌”事故, 对渤海海洋环境造成巨大影响, 此次溢油事件持续事件长, 影响面积大,是对渤海海域影响最大的一次溢油。

本文将业务化气象数值模型(Weather Research and Forecast, WRF)、海流数值模型(Regional Ocean Model System, ROMS)的数值预报结果作为海洋环境驱动场, 采用“油粒子”的海上溢油形态数值模拟方法, 结合油粒子的随机走动算法, 对海上溢油扩散范围等要素进行模拟预测。本文以蓬莱19-3平台溢油事件为例, 设计不同模型分辨率、不同风、流系数对溢油漂移扩散模型进行敏感试验, 获得适合于渤海的溢油数值模型参数。利用优化的溢油漂移扩散模型对“南海一号”平台溢油情况进行数值模拟预测, 并与飞机实测资料进行对比, 得到的预测结果与实际观测结果基本一致。

1 渤海海域海洋环境要素预报模型

在北海预报中心业务化运行的北海区海洋、气象环境动力预报模式的基础上, 采用嵌套技术研发渤海精细化海洋、气象环境动力预报模式。气象预报采用WRF中尺度气象模式, 同时利用3维变分同化技术, 进一步提高模式预报准确度。海流预报采用ROMS区域海洋模式, 进行海流(潮流和环流)、水位、温度和盐度等预测。

1.1 渤海海域气象场预报模型

1.1.1 WRF气象模型简介

北海预报中心建立的WRF气象模型是以美国大气研究中心(NCAR)开发的 WRF-ARW(AdvancedResearch WRF)版本为基础。WRF模式程序具有操作的可移植性、可维护性、扩展性、易读性、运行结构性和互用性等特点。为提高大气分析场的精度,NMFC-WRF采用三维变分同化方法, 将观测资料与数值预报产品结合。

1.1.2 WRF气象模型参数设置

模型采用美国国家环境预报中心(NCEP)提供的FNL(Final)资料为初始场和边界条件, 资料空间分辨率为1°×1°, 时间间隔为6 h。模式结果每小时输出1次。利用WRF的网格嵌套功能, 在溢油事发海域建立高分辨率的气象强迫场, 更加准确地模拟溢油海区的风场情况。其中, 大区为北海区(图1所示范围),空间分辨率为15 km; 小区为渤海海域(图1中黑色框范围), 空间分辨率为5 km。图2为利用WRF模型模拟渤海风场的预测图。

图1 WRF气象场模拟范围示意图Fig.1 Simulation range of WRF model

1.1.3 风场模拟结果验证

为检验模拟预报结果的准确度, 将风场模型预报结果与渤海中部浮标实测数据进行对比验证, 浮标位置见图3, 其中BQF103, BQF109, BQF110分别表示浮标编号, 19-3C为蓬莱19-3C油田位置。图4和图5分别是BQF109和BQF110两个浮标2011年7月1~4日(72 h), 每小时1次的观测数据与模拟结果的对比图。可以看出, 模拟和实测的风速、风向的变化趋势较为一致, 且大小也较为符合, 模拟结果较好地反映了渤海风场特征。

图2 WRF模型模拟渤海风场预测图Fig.2 The forecast of the Bohai Sea wind field with WRF model

图3 渤海中部浮标的站位图Fig.3 The bitmap of buoy stations in the Bohai Sea

1.2 渤海海域海流场预报模型

1.2.1 ROMS海流模型简介

渤海海流业务化预报采用ROMS区域海洋海流模式。该模式是自由表面的、基于流体力学原始方程的静力平衡模式, 在现在海洋学界应用广泛[1-2]。Roms由S-coordinate Rutgers University Model (SCRUM)发展而来[3], 在数值方法[4]和模式架构上都有所改进, 以方便于更高效地完成从单线程到多线程的计算任务。它包括了一系列新的特性: 可选的高阶平流格式, 更准确的压强梯度算法, 更多的次网格尺度参数化方法, 大气、海洋、深海边界条件, 生态模块和数据同化模块。

图4 数值模拟与浮标 BQF109观测的 2011年 7月 1~4日(72 h)风速风向对比Fig.4 The comparison of wind speed and direction between simulation results and BQF109 data from July 1 to 4 in 2011

图5 数值模拟与浮标 BQF110观测的 2011年 7月 1~4日(72 h)风速风向对比Fig.5 The comparison of wind speed and direction between simulation results and BQF110 data from July 1 to 4 in 2011

1.2.2 ROMS海流模型参数设置

根据渤海所在位置和高精度、长时效的业务化需求, 建立2个区域模式: 大区域是整个东中国海海域 (117°34′~127°54′E, 31°54′~41°N), 水 平 分 辨 率(1/30)°×(1/30)°, 垂向 6 层; 小区域是渤海海域(117°30′~122°18′E, 37°~41°N), 水平分辨率(1/240)°×(1/240)°, 垂向 6层。大区主要采用环流场和温盐场驱动。利用美国 HYCOM(Hybrid Coordinate Ocean Model )业务化预报的水位、流速和温盐场作为大区的初值和边界值。小区采用业务化预报北海区WRF风场和热通量场驱动模式。另外, 模型考虑了渤海沿岸4个主要河流(黄河、辽河、海河和滦河)月平均径流量, 小区 ROMS模式外海开边界条件水位采用M2、S2、K1、O1等8个分潮来驱动。

1.2.3 流场模拟结果验证

为检验模拟预报结果的准确度, 将流场模型预报结果与渤海浮标实测数据(浮标位置见图3)进行对比。图6为ROMS模型模拟渤海海流场。图7~图9分别是3个浮标(BQF103、BQF109和BQF110)连续72 h(2011年7月1~4日)每个小时1次的观测数据与模拟结果的对比图。可以看出, 模拟的表层流速和实际观测的数据较符合, 海流场模拟结果较好地反映了渤海海流特征。

2 溢油漂移扩散预测结果分析

2.1 溢油漂移扩散数值模型

石油进入海洋后, 不仅会随着海洋环流和海表面的风进行空间位置的输运, 还存在着极其复杂的物理、化学和生物过程。如石油在海面上的分散、运移、蒸发、溶解、光分解、生物降解、乳化、悬浮物的吸附和沉积过程等[5-6]。溢油在海表面的过程可以分为动力学过程和非动力学过程 2种, 其中溢油动力学过程的驱动场由海流和气象模型提供, 非动力学过程则由蒸发、乳化、溶解等物理、化学变化组成。根据前人的研究成果, 溢油形态学的数值模拟方法是利用不同的数学模型进行计算, 其计算理论可分为油膜扩展模式、对流扩散模式和“油粒子”模式3类[7,8]。

溢油模块采用国际上广泛应用的“油粒子”模型[9-12], 该模型可以很好地模拟溢油进入水体后产生的复杂现象, 包括扩展、漂移、扩散等油膜组分保持恒定的输移过程, 以及输移和风化过程中黏度、表面张力等油膜属性随油膜组分和温度变化而发生的变化。“油粒子”模型基于拉格朗日输运概念, 将溢油离散为大量的油粒子, 每个油粒子代表一定的油量,油膜就是由这些大量的油粒子所组成的“云团”, 模型首先计算每个油粒子的空间位置变化及化学过程变化, 然后统计各网格上的油粒子数, 从而模拟出油膜浓度的时空分布和油量变化。

“油粒子”的物理过程是指由于平流流动和湍流波动引起的粒子运动。平流运动是指每个粒子在特定的流场条件下发生的平移, 适宜用拉格朗日法模拟, 湍流波动是指由于剪切流和湍流引起的扩散运动, 适宜用随机游动法模拟。湍流可视为随机流场,而每个模型粒子在湍流场中的运动则类似于分子的布朗运动。每个粒子的随机运动导致了整个云团在水体中的扩散过程。这种方法实际上是确定性方法和随机方法的结合, 即采用确定性方法模拟平流运动, 采用随机性方法模拟扩散过程[11-12]。

图6 渤海海流预报结果Fig.6 The current forecasting for the Bohai Sea

图7 BQF103浮标表层实测流速和流向与模式预报结果的比较Fig.7 The comparison of surface current speed and direction between BQF103 buoy data and simulation results

图8 BQF109浮标表层实测流速和流向与模式预报结果的比较Fig.8 The comparison of surface current speed and direction between BQF109 buoy data and simulation results

图9 BQF110浮标表层实测流速和流向与模式预报结果的比较Fig.9 The comparison of surface current speed and direction between BQF109 buoy data and simulation results

另外, 模型还考虑了影响油膜变化的非动力过程: 蒸发是由于石油烃类从液态向气态的相变而造成的油膜与空气之间的物质交换, 蒸发导致溢油中的低沸点烃类迅速减少; 乳化是一个油包水的过程,乳化物具有较高的密度和黏性, 影响溢油的扩散过程。这些溢油的风化过程对海洋生态环境系统带来了重大的损害, 造成长期的影响[13-15]。

2.2 溢油数值模型敏感试验

2.2.1 敏感试验设计

近年来, 溢油被看成是由大量“油粒子”组成的,每个油粒子代表一定的油量, 在表层海流和风场的作用下漂移运动。溢油漂移过程常指油膜在水面流场、波浪、风力等外界因素的共同作用下, 所发生的迁移过程。多年来有关溢油漂移的研究普遍认为: 油膜在水体表面漂移的主要驱动力来自于水体表层流场与表面风场。因此, 海洋环境背景场的预测差异对溢油漂移过程的影响较大。为获得渤海海域的溢油数值模型的参数, 提高溢油漂移预测的模拟精度,本文设计了 2个敏感试验, 即溢油数值模型分辨率和风流系数对油膜预测结果的影响, 以优化溢油漂移扩散模型。验证资料为蓬莱19-3平台溢油事故中,2011年7月24日18: 33和8月26日5: 08的2幅SAR(Synthetic Aperture Radar)卫星影像解译的油膜监测数据。

针对溢油数值模型分辨率敏感试验, 设计了模型水平网格分辨率为400 m×400 m和20 m×20 m 2种情况, 优化溢油漂移扩展范围的预测。

根据本文 2.1部分和 2.2部分中对业务化 WRF风场和 ROMS流场的验证结果, 设计不同风流系数对风场和流场进行修正, 优化现有溢油漂移扩散模型的预测结果, 获得适用于渤海的溢油参数。设计了风(流)系数为0.8(0.5)、0.3(0.5)、0.3(0.2) 3种情况。

2.2.2 敏感试验结果分析

图10是水平分辨率为400 m×400 m和20 m×20 m时, 溢油数值模型模拟的2011年7月24日和8月26日的结果对比图。图10中阴影部分表示卫星遥感观测溢油范围, 红线表示分辨率为400 m×400 m的模型模拟溢油范围, 蓝线表示分辨率为20 m×20 m的模型模拟溢油范围。从图中可以看出, 模型分辨率为400 m×400 m时, 对油膜影响范围、形状等的模拟效果较差, 很难反映油膜在海水中的整体状态, 特别是对少量溢油, 模拟误差将大于实际溢油的扩散状态; 而模型分辨率为20 m×20 m时, 能够较好地反映油膜在海水中的分布, 与实际监测信息基本吻合。

图11是风(流)系数分别为 0.8(0.5)、0.3(0.5)和0.3(0.2)时, 溢油数值模型模拟的2011年7月 24日和8月26日的结果对比图, 网格分辨率选择20 m×20 m。图中阴影部分实线表示卫星遥感观测溢油范围, 蓝线表示风(流)系数为 0.3(0.5)的溢油模型模拟结果, 红线表示风(流)系数为 0.8(0.5)的溢油模型模拟结果, 绿线表示风(流)系数为 0.3(0.2)的溢油模型模拟结果。从图中可以看出, 在该海区, 选择风系数为0.3、流系数为0.5的溢油模型模拟结果能更好地反映溢油在海面上的状态, 与实际监测资料吻合较好;另外2种系数下模拟的溢油范围都有不同的偏差。

综上所述, 根据模型分辨率的敏感性试验结果,优化现有模型, 采用区域嵌套的方法, 根据不同溢油位置, 选择较小区域进行网格加密, 模型最小分辨率可达20 m。风(流)系数敏感性试验的结果表明, 环境背景场的差异对溢油漂移方向、扩散范围有很重要的作用。根据上述试验, 我们选择风系数 0.3, 流系数 0.5,修正和优化现有模型, 提高溢油漂移扩散预报的精度。

2.3 溢油数值模型预测结果验证

2.3.1 2010 年5月“南海一号”平台溢油漂移预测验证

2010年5月13 日, 中海油天津分公司报告“南海一号”钻井平台在作业时钻透相邻采油平台的套管, 发生溢油事故。本文根据5月13日、5月14日海监3807飞机监测的海面溢油分布情况和5月15日5时29分过境的中国遥感卫星雷达监测的溢油情况, 利用上文得到的优化的渤海溢油数值模型, 对平台溢油的漂移方向、主要路径和影响范围进行模拟预测。

将5月13日海监3807飞机监测的海面溢油情况作为模型初始场, 结合当时海洋环境要素场, 进行 48h溢油漂移轨迹和扩散面积的模拟预测, 预测的结果与实测对比如图12所示。图中实线分别表示5月13日和5月14日海监飞机监测溢油范围, 虚线表示模拟预测24 h和48 h溢油可能位置, 十字线表示模拟预测溢油中心位置的 48h的漂移轨迹, 阴影部分表示模拟溢油48 h内可能扩散的影响范围。从图12中可以看出, 预计13~14日海上油膜将向北偏西方向漂移, 中心位置向西北方向移动约9 km, 预测漂移方向和中心位置与14日海监3807飞机现场监测到的溢油情况完全吻合; 模拟预测溢油可能影响面积与海监3807飞机14日现场监测影响范围基本吻合, 溢油影响面积大约为10 km2, 溢油基本向北偏西方向漂移。从图12所示对比结果可以看出, 溢油模拟24 h预测油膜中心点误差约为3 km, 漂移方向偏差小于10°。

图12 “南海一号”溢油漂移轨迹和面积模拟预测结果与监测情况比较Fig.12 The comparison between simulation results of oil spill drift and diffusion area and monitoring for the oil spill event of “No.1 of Nahai ” platform

2.3.2 2011 年“蓬莱19-3”平台溢油油膜形态短期模拟验证

2011年6 月渤海蓬莱19-3油田发生严重的溢油泄漏事件, 2011年6月4日蓬莱19-3B平台油藏持续发生渗油事故, 6月17日蓬莱19-3C平台发生井涌,大量溢油暴发入海。此次溢油事故发生地点在海底,大量原油从海底喷出, 造成水体中大量的油污, 而且持续时间较长, 直至10月份, 监测仍然发现19-3平台附近有油膜出现。

本文选取2011年7月24日18时33分SAR卫星影像解译的油膜监测数据作为验证资料。该 SAR卫星影像特点是油膜由平台连续溢油形成, 溢油在海流和海面风的共同作用下, 呈条带状分布, 油带一端与平台相连。直至卫星成像时, 已经持续溢油3~4 h。以平台位置作为溢油模拟的起始位置, 结合当时海洋环境要素场, 模拟4 h内油膜在海表面的扩散情况。图13为溢油模拟结果与卫星监测结果的比较图, 其中实线表示卫星遥感溢油范围, 阴影部分表示模拟预测溢油扩散范围。可以看出, 溢油模拟结果和卫星解译结果在油膜的走向和形态上都非常的一致, 模拟结果较好。但油膜影响面积的模拟结果比监测结果稍大, 溢油扩展面积预报相对误差小于20%。

图13 溢油模拟结果与卫星监测结果比较图Fig.13 The comparison chart of oil spill simulation result with the satellite monitoring

3 结论

本文将业务化气象数值模型WRF、海流数值模型 ROMS的数值预报结果作为海洋环境驱动场, 以蓬莱 19-3平台溢油事件为例, 设计不同溢油敏感试验, 优化和完善溢油漂移扩散预测模型, 对渤海海上溢油的漂移路径、扩展面积和扫海区域等要素进行预测。提高了溢油漂移扩散预报的精度, 为海洋溢油应急处置和防灾减灾提供技术支持。得到的结论如下所述:

通过模型分辨率敏感性试验, 优化现有模型,采用区域嵌套的方法, 根据不同溢油位置, 选择较小区域进行网格加密, 模型最小分辨率可达20 m。

风(流)系数敏感性试验结果表明, 环境背景场的差异对溢油漂移方向、扩散范围有很重要的作用。根据上述试验, 我们选择能较好反映溢油漂移扩散结果的风系数0.3, 流系数0.5。

利用优化的溢油漂移扩散模型对“南海一号”平台溢油情况进行模拟预测, 并与飞机实测资料进行对比验证, 得到预测结果与实际观测结果吻合。

[1]Haidvogel D B, Arango H G, Hedstrom K, et al.Model evaluation experiments in the North Atlantic Basin:Simulations in nonlinear terrain-following coordinates[J].Dyn Atmos Oceans, 2000, 32: 239-281.

[2]Ezer T, Arango H G, Shchepetkin A F.Developments in terrain-following ocean models: Intercomparisons of numerical aspects[J].Ocean Model, 2002, 4: 249-267.

[3]Song Y, Haidvogel D.A semi-implicit ocean circulation model using a generalized topography-following coordinate system[J].J Comput Phys, 1994, 115: 228-244.

[4]Shchepetkin A F, Mcwilliams J C.The regional oceanic modeling system (ROMS): A split-explicit, free-surface,topography-followingcoordinate oceanic model[J].Ocean Model, 2005, 9: 347-404.

[5]Fay J A.Physical processes in the spread of oil on a water surface[C]// Conf Prevention and Control of oil Spill.Washington D C: American Petroleum Institute,1971: 463-467.

[6]杨庆霄, 徐俊英, 李文森, 等.海上溢油溶解过程的研究[J], 海洋学报, 1994, 16(3): 50-56.

[7]武周虎, 赵文谦.海面溢油扩展、离散和迁移的组合模型[J].海洋环境科学, 1992, 11(3): 34-40.

[8]李冰绯.海上溢油的行为和归宿数学模型基本理论与建立方法的研究[D].天津: 天津大学, 2004.

[9]Chen H Z, Li D M, Li X.Mathematical modeling of oil spill on the sea and application of the modeling in Daya bay[J].Journal of Hydrodynamics, 2007, 19(3): 282-291.

[10]Reed M, Johansen O, Brandvik P J, et al.Oil spill modeling towards the close of the 20th century :over-view of the state of the art [J].Spill Science &Technology Bulletin, 1999, 5(1): 3-16.

[11]娄安刚, 王学昌, 孙长青, 等.胶州湾海面溢油轨迹的数值模拟[J].黄渤海海洋, 2001, 19(1): 1-8.

[12]刘伟峰, 孙英兰.海上溢油运动数值模拟方法的探讨与改进[J].华东师范大学学报, 2009, 3: 96-103

[13]张存智, 窦振兴, 韩康, 等.三维溢油动态预报模型[J].海洋环境科学, 1997, 16(1): 22-29.

[14]Elhakeem A A, Elshorbagy W, Chebbi R.Oil spill simulation and validation in the Arabian (Persian).Gulf with special reference to the UAE coast[J].Water Air Soil Pollution, 2007, 184: 243-254.

[15]Wang S D, Shen Y M, Guo Y K, et al.Threedimensional numerical simulation for transport of oil spill in seas[J].Ocean Engineering, 2008, 35: 503-510.