基于MIKE21的海洋溢油数值模拟在长江口区域的应用

1、基于MIKE21的海洋溢油数值模拟在长江口区域的应用殷健，李丕学，何金林（上海市水文总站 上海 200232）摘要：本文介绍了溢油数值模拟的基本原理，利用MIKE21建立长江口杭州湾溢油模型，针对518溢油事故进行后评估，并通过6种不同工况的模拟研究溢油扩散的影响因子，讨论长江口区航道内船舶溢油事故对长江口青草沙水源地的潜在影响。关键词：长江口；二维数值模拟；溢油模型；油粒子Application of ocean oil spills numerical simulation based on MIKE21 model in the Yangtze EstuaryLI Pi-xue1, HE

2、 Jin-lin1, YIN Jian1，（1.Shanghai Hydrology Administration, Shanghai 200232, China）Abstract：This paper introduces the basic principles of the oil spills simulation. Building the Yangtze Estuary Oil Model using MIKE21, we studied the situation of the 518 oil spill accident, and the influencing factors

3、 of the oil diffusion in the sea in six cases. In the end, we discussed potential impact of oil spill accidents in Yangtze Estuary on Qing Cao Sha water source area.Keywords：the Yangtze Estuary; two-dimensional numerical simulation; oil spills model; oil particle1. 引言据统计，每年全球平均发生千吨级以上的溢油事故4.4起，百吨级溢油

4、事故近百起(王长海，2000)，我国则年均发生500起海洋溢油事故。导致海洋溢油事故的最主要原因是海上运输，每年因航运产生的石油泄露达160至200万吨(陈莎，2008)。近年来部分沿海地区的海水含油量已超标2至8倍(岳文浩 等，2008)，海洋石油污染形势十分严峻。溢油的风险概率与水域状况、天气条件、船只密度、航道管理等因素相关。溢油事故发生后，需要迅速、准确地预测溢油的变化趋势，为科学、有效的决策提供依据。溢油的数值模拟指基于先期计算的溢油水域流场，结合实地实时的溢油品种、溢油量、风、浪、盐度、水温、气温等因素，预测溢油在迁移与风化过程中的状态、组分、性质的改变以及最终归宿，为应急决策的制

5、定、清除手段的选择和溢油损害的评估提供科学依据。2. 长江口区域的溢油风险长江口地处长江东西运输通道与海上南北运输通道的交汇点，是上海港及长江干线的通海咽喉，在内陆航运及海路航运中占据重要地位，是我国通航密度最高的水域之一，年货运吞吐量超过7亿吨，日年均船舶量达2300多艘，高峰时期更可达6000多艘(陈莎，2008)，并以大型船舶为主，发生溢油事故的潜在风险较大。同时长江口也是上海最重要的水源地-青草沙水库所在地，区域内有崇明东滩鸟类自然保护区、中华鲟自然保护区、九段沙湿地自然保护区等全国重要的生态保护区，研究本区域船舶溢油对溢油应急处置具有重要的现实意义。本文主要讨论长江口区航道内船舶溢油

6、事故对长江口青草沙水源地的潜在影响。3. 长江口水动力数学模型建立高效稳定、精确可靠的水动力数学模型是进行海洋溢油行为与归宿数值模拟的重要前提。MIKE21HD包括了广泛的水力现象，可有效模拟由于各种作用力而产生的水位及水流变化，进而为溢油模拟模型MIKE21SA提供可靠的水动力计算基础。本模型区域包括长江口、杭州湾以及邻近海区；采用非结构网格；模型计算时间步长为30 s，网格节点26121个，单元数47174个，计算时采用高阶求解格式(图1)。模型边界条件分别对应径流、潮汐、海表面风应力及床底摩擦应力。模型主要参数包括床面阻力系数、紊动粘性系数、引潮势等，其中床面阻力系数反映水流与床面相互作

7、用过程中床面边界的形态及粗糙程度对水流阻力的综合影响(许婷，2010)，采用实测资料反推获得；紊动粘性系数的选择对岸线变化复杂且有回流产生的岸段十分关键，采用经验公式获取。图1 模型计算区域及水深分布选取徐六泾、白茆、高桥、小衢山、大戢山等站点率定模型潮位；选择徐六泾、东风西沙、北港、连兴港、中浚等位置率定模型潮流，水动力模型的率定验证结果表明模型能较为准确地模拟区域潮汐潮流运动(图3-图9)。图2 白茆站潮位过程图3 高桥站潮位过程图4 小衢山站潮位过程图5 大戢山站潮位过程图6 北港大潮潮流流向过程图7 北港大潮潮流流速过程图8 连兴港小潮潮流流向过程图9 连兴港小潮潮流流速过程4. 长江

8、口溢油数学模型4.1 溢油过程的数值模拟溢油扩散受其物理、化学、生物过程的影响，而这些过程又与溢油的性质、水动力环境、气象环境等因素密切相关(刘伟峰 等，2009)。海面溢油发生后，首先受引力、表面张力、惯性等影响，以薄膜形式散开；随后在海流、波浪、风场的共同作用下，产生漂移，对油膜漂移的精确模拟是迅速有效地清除溢油，将污染降至最低程度的重要前提。同时，石油烃的蒸发是溢油质量传输的主要过程，其增加油的比重、倾点、粘性及表面张力。油膜在迁移与风化的同时还发生溶解过程和光氧化过程，氧化产物虽浓度低但毒性强，长期积累对水生生物及水体环境危害较大(图10-图11)。图10 海洋溢油行为与归宿图11 海

9、洋溢油过程时间尺度MIKE21SA溢油模型基于欧拉-拉格朗日理论体系，通过对溢油在水体中的扩展、漂移、蒸发、分散、沉降、乳化、溶解、光氧化等过程的模拟，获取油膜的漂移过程及厚度增减，同时计算溢油的比重、倾点、粘性等性质的变化。溢油模型基于油粒子模式进行数值模拟，油粒子模式是将溢油视为大量质量不等的粒子所组成的集合，并以粒子的宏观运动来表达溢油在水环境中行为过程的一种模拟方法(刘生宝，2007)。油粒子模式将溢油离散为大量油粒子，每个油粒子代表一定油量，借助油粒子的随机运动反映油膜的扩展与漂移，通过油粒子的质量改变体现油膜的蒸发、乳化、分散、沉降等过程，并通过计算设定域内油粒子的数量、体积获取油

10、膜的厚度分布。油粒子模式可更有效地响应海洋水动力条件，较为精确的模拟油膜运动过程中顺风拉长、迎风压缩以及扭曲断裂等变化。溢油模型首先计算各油粒子的组分和位置，分析含水率与体积，随后将运算区域划分网格并于各时间步长末统计域内油粒子的数量和分布，依据油粒子体积及网格面积计算油膜的厚度与浓度，在模拟溢油输移与风化的同时，借助热量迁移计算油膜的热容量、表面张力、粘度等物理化学性质的变化。4.2 模型构建与参数设置溢油模型的运算基于水动力模型，两者的计算区域相同。溢油模型网格分辨率取100 m，模型模拟的时间步长小于网格间距与计算域最大流速的比值，模拟时长30 h。模型考虑的条件场及参数因子包括溢油特征

11、(源通量、油粒子数量、油组分特征、油粘度)；水体特征(水温场、盐度场)；大气特征(气温场、阴暗度)；扩散(横向水平扩散、纵向水平扩散、垂向扩散)；对流(风场、风因子、风偏转角、风深度)；涡流(对数速度分布)；蒸发(蒸发率)；溶解夹带(传质系数、油包水界面张力)；乳化(油最大含水率、吸收系数、释放系数)；热量迁移(油初始温度、油辐射率、水辐射率、大气辐射率、漫射率)。5. 长江口船舶溢油模拟5.1 518溢油事故模拟2012年5月18日夜间，通银6加油船装载重油于吴淞口东长江口6#锚地进水沉没，事发位置1213610.46E、312331.75N，距离上游青草沙水库取水泵闸约17 km，可能对青

12、草沙水源地造成影响。本文设定本次事故为瞬时溢油，溢出25 m3重油，在此基础上模拟溢油发生后30 h内油膜位置、面积、形状、厚度等变化以及溢油蒸发、溶解、乳化、沉降等趋势，与本次溢油监测结果进行对比，初步探讨溢油的迁移与风化规律(图12-图16)。图12 518溢油事故位置图13 30 h后油膜分布图14 油膜迁移时段与范围图15 油膜平均厚度变化曲线图16 油膜10 m厚度包络面积变化曲线模拟成果表明，所构建的油粒子模式溢油模型能较好地反映溢油的迁移与风化规律，适用于长江口区域的溢油模拟：(1) 油膜扩展主要受溢油自身物理化学性质、潮流场、风场等因素的影响。在溢油发生后初期，油膜在重力及惯性

13、力作用下，覆盖面积逐步增大，平均厚度迅速减小，这一阶段油膜呈椭球型均匀扩展，随着时间推移，潮流逐渐成为支配油膜扩展的主导因素，油膜在海面被潮流拉成彗星状油膜带，并随潮汐来回漂移。(2) 油粒子模式中，油膜在海面的漂移受表层流场及表面风场驱动，考虑到二维水动力模型计算的潮流为垂向平均值，模型建立对数函数关系通过垂向平均流速推算表层流场，并借助数值积分获取油粒子团的运动轨迹。潮流与风应力的矢量叠加决定了油膜的漂移方向、速度、距离，油膜的漂移随表层流速的增大而加速，随表层流速的减小而减速或转向，在涨、落急以及涨、落憩等潮流特征时段表现的尤为显著，与此同时表面风场的存在有效促使了油膜向下风方向漂移。溢

14、油迁移时段与范围表明，溢油发生后30 h内油膜未对青草沙水库取水泵闸产生影响。(3) 溢油在进行扩展、漂移的同时发生蒸发与乳化。前期蒸发随油膜面积的增大而增强，在易蒸发的轻组分快速蒸发后，溢油的蒸发过程逐渐减弱，蒸发速率逐步减缓趋于稳定。在波动产生的混合能量作用下，溢油不断乳化，含水率增加，有效体积及覆盖面积增大，约16 h后含水率的变化趋缓，并逐渐稳定在80%左右。溢油的蒸发及乳化现象导致溢油的密度变化，油膜密度因蒸发过程造成的轻组分损失以及乳化过程引起的含水率增加而呈上升趋势，并随时间推移与周围水体的密度逐渐接近。5.2 溢油模拟的影响因子环境是一个开放性系统，自然因素及人类活动都会对环境

15、产生各种复杂的影响，这些影响往往难以精确定量，数值模拟中模型结构、输入条件、参数取值等也存在一定的不确定性(姜卫星，2007)。借助影响因子分析建立低灵敏度系统，可有效提高模型运行的准确性，通过比对不同输入条件下的模拟结果，掌握各种状况下油膜输移运动的规律，从而可在事故发生后有针对性地收集模拟预测所需的关键信息。溢油数值模型主要用于突发性溢油事故的模拟预测，目的是为事故的应急响应提供支持，油膜的运动轨迹、溢油影响范围、油膜到达敏感区域的时间等是模型模拟预测的重点，这些指标主要受溢油位置、溢油类型、溢油方式、溢油量、潮流场、潮汐形态、风场、温度场、盐度场等因素的影响，本文选取溢油量与风场因素进行

16、影响因子分析，通过比对长江口历年重、特大溢油事故案例，溢油量因子选择50 m3、100 m3；考虑到长江口最频繁风向为NW-N及ESE-SSE，频率分别达24%和23%，年平均风速3.7 m/s，风场因子选取长江口最频繁风向以及年平均风速，在此基础上建立工况1-工况6进行溢油模拟运算(表1-表2，图17-图24)。表1 长江口溢油事故案例序号事故时间事故位置事故溢油量1高桥20 t2吴淞口32#灯浮100 t3吴淞口182 t4吴淞口锚地500 t5崇明新华14 t6长兴岛130 t7深水航道W-2灯浮30 t8长江口3号锚地20 t表2 溢油模拟工况工况溢油量风场150 m3静风2100 m

17、3350 m3ESE-SSE 3.70 m/s4100 m3550 m3NW-N 3.70 m/s6100 m3图17 工况4 30 h后油膜分布图18 工况4油膜迁移时段与范围图19 工况4油膜平均厚度变化曲线图20 工况4油膜10 m包络面积变化曲线图21 工况5 30 h后油膜分布图22 工况5油膜迁移时段与范围图23 工况5油膜平均厚度变化曲线图24工况5油膜10 m包络面积变化曲线通过对比不同工况的模拟成果，进而可分析溢油量及风场因素对油膜迁移与风化过程的影响：(1) 溢油量的增加提高了油膜厚度与覆盖范围，而风场则在溢油的迁移、风化过程中起显著作用。静风环境下，油膜大体随潮汐涨落方向

18、漂移，风场的出现促进油膜向下风方向移动，且风速与流速的大小不同、风向与流向的夹角不同，对油膜漂移的影响不同。工况3、4中涨潮阶段风向与流向反向，减缓了油膜向上游漂移的速度，落潮阶段风向与流向同向，油膜加速向下游漂移；工况5、6中涨潮阶段风向与流向同向，加快了油膜向上游漂移的速度，落潮阶段风向与流向反向，油膜减速向下游漂移，相比静风条件，各时段油膜位置均偏向西北。(2) 油膜厚度变化历经两个阶段，第一阶段油膜厚度迅速减小，持续约3 h，第二阶段油膜厚度随时间推移缓慢增减，与之相应油膜面积则逐步扩大。相比静风环境，风场引起的水流湍动及水面波浪影响驱动油膜的蒸发和乳化过程，导致油膜厚度与面积的变化梯

19、度时缓时陡。在潮流场与风场的共同作用下，在工况3、4的转潮时段，油膜各部分受到方向不同的潮流作用，发生扭曲及断裂，并进一步相互分离形成油斑，工况5、6中油膜抵达岸边后受局部地形影响显著变形，模型采用的油粒子模式借助大量独立运动的油粒子，较好地模拟了油膜的扭曲变形及油斑形成。(3) 溢油事故发生后，短时间内油膜在水中的主要行为为扩展及漂移，在突发性溢油事故模拟中，如若溢油源距离保护目标较近，应重点考虑油膜在水中扩展漂移到达保护目标的时间。溢油迁移时段与范围表明，溢油发生后30 h内工况1、2、5、6油膜未对青草沙水库取水泵闸产生影响，溢油发生后21-24 h内工况3、4油膜影响至取水泵闸区域，期间油膜厚度维持在1-10 m区间。6. 结语随着石油海上开采与水上运输的不断增加，近年来海洋溢油事故的次数及数量均日益增多，正确、快