长江口区域溢油事故对青草沙水源地影响的数值模拟

殷 健，陈志铮，梁珊珊

(1. 上海市水文总站，上海 200232;2. 上海市排水管理处，上海 200001)

据统计，每年全球平均发生千吨级以上的溢油事故4.4 起，百吨级溢油事故近百起［1］，我国年均发生500 起海洋溢油事故。导致海洋溢油事故的最主要原因是海上运输，每年因航运产生的石油泄露达160 ～200 万t［2］。近年来部分沿海地区的海水含油量已超标2 ～8 倍［3］，海洋石油污染形势十分严峻。

溢油的风险概率与水域状况、天气条件、船只密度、航道管理等因素相关。溢油事故发生后，需要迅速、准确地预测溢油的变化趋势，为科学有效的决策提供依据。

溢油的数值模拟指基于先期计算的溢油水域流场，结合实地实时的溢油品种、溢油量、风、浪、盐度、水温、气温等因素，预测溢油在迁移与风化过程中的状态、组分、性质的改变以及最终归宿，为应急决策的制定、清除手段的选择和溢油损害的评估提供科学依据。

1 长江口区域的溢油风险分析

长江口地处长江东西运输通道与海上南北运输通道的交汇点，是上海港及长江干线的通海咽喉，在内陆航运及海路航运中占据重要地位，是我国通航密度最高的水域之一，年货运吞吐量超过7 亿t，日均船舶量达2 300 多艘，高峰时期更可达6 000 多艘［2］，并以大型船舶为主，发生溢油事故的潜在风险较大。

同时长江口也是上海最重要的水源地青草沙水库［4］所在地，区域内有崇明东滩鸟类自然保护区、中华鲟自然保护区、九段沙湿地自然保护区等全国重要的生态保护区，研究本区域船舶溢油对溢油应急处置具有重要的现实意义。

2 数学模型的动力

2.1 长江口水动力数学模型

建立高效稳定、精确可靠的水动力数学模型是进行海洋溢油行为与归宿数值模拟的重要前提。MIKE21HD 包括了广泛的水力现象，可有效模拟由于各种作用力而产生的水位及水流变化，进而为溢油模拟模型MIKE21SA 提供可靠的水动力计算基础。模型计算区域及水深分布如图1，模型区域包括长江口、杭州湾以及邻近海区;采用非结构网格;模型计算时间步长为30 s，网格节点为26 121 个，单元数为47 174 个，计算时采用高阶求解格式。

图1 模型计算区域及水深分布Fig.1 Model's Region and Depth Distribution

模型边界条件分别对应径流、潮汐、海表面风应力及床底摩擦应力。模型主要参数包括床面阻力系数、紊动粘性系数、引潮势等，其中床面阻力系数反映水流与床面相互作用过程中床面边界的形态及粗糙程度对水流阻力的综合影响［5］，采用实测资料反推获得;紊动粘性系数的选择对岸线变化复杂且有回流产生的岸段十分关键，采用经验公式获取。

于模拟区域内选取不同站点进行水动力率定，如图2 所示。模型能较为准确地模拟区域内潮汐潮流运动，如图3 ～图6 所示。

图2 模型水动力率定站点分布Fig.2 Model's Hydrodynamic Calibration Station

图3 高桥站潮位过程Fig.3 Gaoqiao Station's Tide Process

图4 大戢山站潮位过程Fig.4 Dajishan Station's Tide Process

图5 中浚大潮潮流流向过程Fig.5 Zhongjun Tide's Flow Process

图6 北港大潮潮流流速过程Fig.6 Beigang Tide's Velocity Process

2.2 长江口溢油数学模型

2.2.1 溢油过程的数值模拟

溢油在海洋上的扩散受其物理、化学、生物过程的影响，而这些过程又与溢油的性质、水动力环境、气象环境等因素密切相关［6］。海面溢油发生后，首先受引力、表面张力、惯性等影响，以薄膜形式散开;随后在海流、波浪、风场的共同作用下，产生漂移，对油膜漂移的精确模拟是迅速有效地清除溢油，将污染降至最低程度的重要前提;石油烃的蒸发是溢油质量传输的主要过程，增加油的比重、倾点、粘性及表面张力;油膜在迁移与风化的同时还发生溶解过程和光氧化过程，光氧化产物虽浓度低但毒性强，长期积累危害较大;由于溢油对水面的覆盖，减少了溢油层以下水体与空气的交换和循环，易导致大面积水体缺氧，进而对水体环境及水生生物正常生长产生不利影响。海洋溢油行为与归宿，如图7 所示。海洋溢油过程时间尺度，如图8 所示。

MIKE21SA 溢油模型基于欧拉-拉格朗日理论体系，通过对溢油在水体中的扩展、漂移、蒸发、分散、沉降、乳化、溶解和光氧化等过程的模拟，获取油膜的漂移过程及厚度增减，同时计算溢油的比重、倾点、粘性等性质的变化。

溢油模型基于油粒子模式进行数值模拟，油粒子模式是将溢油视为大量质量不等的粒子所组成的集合，并以粒子的宏观运动来表达溢油在水环境中行为过程的一种模拟方法［7］。油粒子模式将溢油离散为大量油粒子，每个油粒子代表一定油量，借助油粒子的随机运动反映油膜的扩展与漂移，通过油粒子的质量改变体现油膜的蒸发、乳化、分散和沉降等过程，并通过计算设定域内油粒子的数量、体积获取油膜的厚度分布。油粒子模式可更有效地响应海洋水动力条件，较为精确地模拟油膜运动过程中顺风拉长、迎风压缩以及扭曲断裂等变化。

图7 海洋溢油行为与归宿Fig.7 Behavior and Fate of Marine Oil Spill

图8 海洋溢油过程时间尺度Fig.8 Marine Oil Spill Process

溢油模型首先计算各油粒子的组分和位置，分析含水率与体积，随后将运算区域划分网格并于各时间步长末统计域内油粒子的数量和分布，依据油粒子体积及网格面积计算油膜的厚度与浓度，在模拟溢油漂移与风化的同时，借助热量迁移计算油膜的热容量、表面张力、粘度等物理化学性质的变化。

2.2.2 模型构建与参数设置

溢油模型的运算基于水动力模型，两者的计算区域相同。溢油模型网格分辨率取100 m，模型模拟的时间步长小于网格间距与计算域最大流速的比值，模拟时长为30 h。模型考虑的条件场及参数因子包括溢油特征(源通量、油粒子数量、油组分特征、油粘度);水体特征(水温场、盐度场);大气特征(气温场、阴暗度);扩散(横向水平扩散、纵向水平扩散、垂向扩散);对流(风场、风因子、风偏转角、风深度);涡流(对数速度分布);蒸发(蒸发率);溶解夹带(传质系数、油包水界面张力);乳化(油最大含水率、吸收系数、释放系数);热量迁移(油初始温度、油辐射率、水辐射率、大气辐射率、漫射率)。

2.3 长江口船舶溢油模拟

2.3.1 5·18 溢油事故模拟

2012 年5 月18 日夜间，通银6 加油船装载重油于吴淞口东长江口6#锚地进水沉没，事发位置121°36'10.46″E、31°23'31.75″N，距离上游青草沙水库取水泵闸约17 km，存在对青草沙水源地造成影响的可能。本次事故为瞬时溢油，溢出约25 m3重油，通过模拟溢油发生后30 h 内油膜位置、面积、形状、厚度等变化和溢油蒸发、溶解、乳化、沉降等趋势(如图9 ～图12)，与溢油事故监测结果进行对比，初步探讨溢油的迁移与风化规律。

图9 模拟30 h 后油膜分布Fig.9 Film Distribution after 30 Hours'Simulation

图10 各时段油膜迁移范围Fig.10 Scope of Film with Time

图11 油膜平均厚度变化曲线Fig.11 Thickness of Film with Time

图12 油膜10 μm 厚度包络面积变化曲线Fig.12 Area of 10 μm Film Cover with Time

以下3 个模拟成果表明所构建的油粒子模式溢油模型能较好地反映溢油的迁移与风化规律，适用于长江口区域的溢油模拟。

(1)油膜扩展主要受溢油自身物理化学性质、潮流场、风场等因素的影响。在溢油发生后初期，油膜在重力及惯性力作用下，覆盖面积逐步增大，平均厚度迅速减小，这一阶段油膜呈椭球形均匀扩展，随着时间推移，潮流逐渐成为支配油膜扩展的主导因素，油膜在海面被潮流拉成彗星状油膜带，并随潮汐来回漂移。

(2)油粒子模式中，油膜在海面的漂移受表层流场及表面风场驱动，考虑到二维水动力模型计算的潮流为垂向平均值，模型建立对数函数关系通过垂向平均流速推算表层流场，并借助数值积分获取油粒子团的运动轨迹。潮流与风应力的矢量叠加决定了油膜的漂移方向、速度、距离，油膜的漂移随表层流速的增大而加速，随表层流速的减小而减速或转向，在涨、落急以及涨、落憩等潮流特征时段表现得尤为显著，与此同时表面风场的存在有效促使了油膜向下风方向漂移。溢油迁移时段与范围表明，溢油发生后30 h 内油膜未对青草沙水库取水泵闸产生影响，与溢油事故监测结果一致。

(3)溢油在进行扩展、漂移的同时发生蒸发与乳化。前期蒸发随油膜面积的增大而增强，在易蒸发的轻组分快速蒸发后，溢油的蒸发过程逐渐减弱，蒸发速率逐步减缓趋于稳定。在波动产生的混合能量作用下，溢油不断乳化，含水率增加，有效体积及覆盖面积增大，约16 h 后含水率的变化趋缓，并逐渐稳定在80%左右。溢油的蒸发及乳化现象导致溢油的密度变化，油膜密度因蒸发过程造成的轻组分损失以及乳化过程引起的含水率增加而呈上升趋势，并随时间推移与周围水体的密度逐渐接近。

2.3.2 溢油模拟的影响因子

环境是一个开放性系统，自然因素及人类活动都会对环境产生各种复杂的影响，这些影响往往难以精确定量，数值模拟中模型结构、输入条件、参数取值等也存在一定的不确定性［8］。借助影响因子分析建立低灵敏度系统，可有效提高模型运行的准确性，通过比对不同输入条件下的模拟结果，掌握各种状况下油膜输移运动的规律，从而可在事故发生后有针对性地收集模拟预测所需的关键信息。

溢油数值模型主要用于突发性溢油事故的模拟预测，目的是为事故的应急响应提供支持，油膜的运动轨迹、溢油影响范围、油膜到达敏感区域的时间等是模型模拟预测的重点，这些指标主要受溢油位置、溢油类型、溢油方式、溢油量、潮流场、潮汐形态、风场、温度场、盐度场等因素的影响［9］，本文选取溢油量与风场因素进行影响因子分析，通过比对长江口历年重、特大溢油事故案例(如表1)，溢油量因子选择50、100 m3;考虑到长江口最频繁风向为NW-N 及ESE-SSE，频率分别达24%和23%，年平均风速3.7 m/s，风场因子选取长江口最频繁风向以及年平均风速，在此基础上于事故位置建立工况1 ～工况6(如表2)进行溢油模拟运算(如图13 ～图26)，讨论不同工况溢油事故对青草沙水源地的潜在影响。

表1 长江口溢油事故案例Tab.1 Oil Spill Case of Yangtze Estuary

表2 溢油模拟工况Tab.2 Simulated Conditions of Oil Spill

图13 工况1 下30 h 后油膜分布Fig.13 Film Distribution under Condition 1 after 30 Hours'Simulation

图14 工况1 下各时段油膜迁移范围Fig.14 Scope of Film under Condition 1 with Time

图15 工况2 下30 h 后油膜分布Fig.15 Film Distribution under Condition 2 after 30 Hours'Simulation

图16 工况2 下各时段油膜迁移范围Fig.16 Scope of Film under Condition 2 with Time

图17 工况3 下30 h 后油膜分布ig.17 Film Distribution under Condition 3 after 30 Hours'Simulation

图18 工况3 下各时段油膜迁移范围Fig.18 Scope of Film under Condition 3 with Time

图20 工况4 下各时段油膜迁移范围Fig.20 Scope of Film under Condition 4 with Time

图21 工况5 下30 h 后油膜分布ig.21 Film Distribution under Condition 5 after 30 Hours'Simulation

图22 工况5 下各时段油膜迁移范围Fig.22 Scope of Film under Condition 5 with Time

图23 工况6 下30 h 后油膜分布ig.23 Film Distribution under Condition 6 after 30 Hours'Simulation

图24 工况6 下各时段油膜迁移范围Fig.24 Scope of Film under Condition 6 with Time

图25 工况1 ～6 下油膜平均厚度变化曲线Fig.25 Thickness of Film under Condition 1 ～6 with Time

通过对比不同工况的模拟成果，分析溢油量及风场因素对油膜迁移与风化过程的影响，具体分析如下。

图26 工况1 ～6 下油膜10 μm 包络面积变化曲线ig.26 Area of 10 μm Film Cover under Condition 1 ～6 with Time

(1)溢油量的增加提高了油膜厚度与覆盖范围，而风场则在溢油的迁移、风化过程中起显著作用。静风环境下，油膜大体随潮汐涨落方向漂移，风场的出现促进油膜向下风方向移动，且风速与流速的大小不同、风向与流向的夹角不同，对油膜漂移的影响不同［10］。工况3、4 中涨潮阶段风向与流向反向，减缓了油膜向上游漂移的速度，落潮阶段风向与流向同向，油膜加速向下游漂移;工况5、6 中涨潮阶段风向与流向同向，加快了油膜向上游漂移的速度，落潮阶段风向与流向反向，油膜减速向下游漂移，相比静风条件，各时段油膜位置均偏向西北。

(2)油膜厚度变化历经两个阶段，第一阶段油膜厚度迅速减小，持续约3 h。第二阶段油膜厚度随时间推移缓慢减小，与之相应油膜面积则逐步扩大。相比静风环境，风场引起的水流湍动及水面波浪影响驱动油膜的蒸发和乳化过程，导致油膜厚度与面积的变化梯度时缓时陡。在潮流场与风场的共同作用下，在工况3、4 的转潮时段，油膜各部分受到方向不同的潮流作用，发生扭曲及断裂，并进一步相互分离形成油斑，工况5、6 中油膜抵达岸边后受局部地形影响显著变形，模型采用的油粒子模式借助大量独立运动的油粒子，较好地模拟了油膜的扭曲变形及油斑形成。

(3)溢油事故发生后，短时间内油膜在水中的主要行为为扩展及漂移，因此，建议在突发性溢油事故发生后，如若溢油源距离青草沙水库较近，应重点考虑油膜在水中扩展漂移到达取水泵闸的时间［11］。溢油迁移时段与范围表明，溢油发生后30 h 内工况1、2、5 和6 油膜未对青草沙水库取水泵闸产生影响，溢油发生后21 ～24 h 内工况3、4 油膜影响至取水泵闸区域，期间油膜厚度维持在1 ～10 μm。

3 结论

运用MIKE21SA 建立长江口杭州湾二维溢油数值模型，针对5·18 溢油事故进行后评估，与事故监测结果的对比表明，所构建的油粒子模式溢油模型能较好反映溢油的迁移与风化规律，可适用于长江口区域的溢油模拟。事故位置不同工况的模拟结果表明溢油量的增加提高了油膜厚度与覆盖范围，而风场则在溢油的迁移、风化过程中起显著作用;油膜厚度变化历经两个阶段，第一阶段油膜厚度迅速减小，持续约3 h，第二阶段油膜厚度随时间推移缓慢减小，与之相应油膜面积则逐步扩大;如若溢油发生在青草沙水源地附近，应重点以油膜的扩展漂移作用评估其扩散至水源地的时间。

［1］王长海.溢油漂移扩散计算模式初步研究［J］.交通环保，2000，21(2):7-9.

［2］陈莎.长江口船舶溢油模拟及其对环境的风险研究［D］.上海:上海海洋大学，2008.

［3］岳文洁，刘书俊.浅析海洋产业的环境损耗核算［J］. 环境与可持续发展，2008，(1):44-46.

［4］上海青草沙水源地工程简介［J］.净水技术，2009，28(3):78.

［5］许婷.丹麦MIKE21 模型概述及应用实例［J］. 水利科技与经济，2010，16(8):867-869.

［6］刘伟峰，孙英兰. 海上溢油运动数值模拟方法的探讨与改进［J］.华东师范大学学报(自然科学版)，2009，(3):90-97.

［7］刘生宝.基于POM 的长江感潮河段溢油数值模拟研究［D］.南京:河海大学，2007.

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

［9］刘彦呈，殷佩海，林建国，等. 基于GIS 的海上溢油扩散和漂移的预测研究［J］.大连海事大学学报，2002，28(3):41-44.

［10］郭运武，刘栋，钟宝昌，等.风对河道溢油扩展、漂移影响的实验研究［J］.水动力学研究与进展A 辑，2008，23(4):446-452.

［11］龙绍桥，娄安刚，谭海涛，等. 海上溢油粒子追踪预测模型中的两种数值方法比较［J］. 中国海洋大学学报(自然科学版)，2006，36(s1):157-162.