危险化学品海上漂移扩散数值模拟研究综述

郭凯旋

（国家海洋环境预报中心，北京100081）

1 引言

随着我国沿海经济的快速发展，临海工业及港口物流规模不断增加，海上危险化学品运输需求量越来越大，危险品的种类和数量也逐年增加。截至2018年12月31日，我国沿海省际运输化学品船（含油品和化学品两用船，下同）共计288艘、112.9万载重吨，同比增加16艘、6.72万载重吨，增幅6.33%[1]；全国万吨级及以上泊位中，专业化液体化工泊位217个，比上年增加12个[2]。

自1990年以来我国发生了多起海上危化品泄漏事故：1997年，Blue Sky No.2散化船在杭州以东200 km洋面沉没，船上载有的988 t酞酸二辛酯泄漏入海；2001年，韩国籍散化船“大勇”轮在长江口外海域和香港散货船“大望”轮相撞，造成630 t苯乙烯泄漏，给海洋生态环境造成了严重损害；2011年，渤海蓬莱19-3油田发生溢油事故，造成油田周边及其西北部面积约6 200 km2的海域海水污染(超第一类海水水质标准)，其中870 km2海水受到严重污染(超第四类海水水质标准)；2012年，韩国籍化学品船“雅典娜”轮在广东汕尾海域沉没，船上7 000 t浓硫酸及140 t燃料油发生泄漏；2018年，巴拿马籍油船“桑吉”轮与香港散装货船“长峰水晶”轮在长江口东约160 n mile处发生碰撞事故，“桑吉”轮爆燃沉没，船上的十多万吨凝析油和大量燃油发生泄漏，据国家海洋局的消息，沉海区域有长约12 km，宽约9 km的油污带，面积约58 km2，严重影响了东海海洋生态环境；2018年，宁波舟山“天桐1号”油轮靠泊泉港东港石化公司码头拟接运工业用裂解碳九，由于操作员违规操作，造成69.1 t碳九泄漏，虽通过吸油毡回收了约40 t泄漏物，但是剩余的大部分泄漏入海，给当地海洋生态环境及海水网箱养殖造成影响。

海上危化品泄漏事故与其他海洋环境污染事故相比具有突发性强、影响范围广、影响时间长和泄漏后变化复杂等特点[3]。突发性强表现为海上危化品泄漏大部分都是在没有任何前兆的情况下突然发生的，瞬间排放大量的污染物，给海洋环境造成巨大影响。影响范围广表现为少量的危化品发生泄漏，其产生的影响范围将非常巨大，例如1 t氯的泄漏能够影响4.8 km2的范围[3]。影响时间长表现为危化品泄漏对海洋生态环境的影响不仅仅发生在泄漏时，部分危化品对海洋生态环境的污染极难消除，例如福岛核泄漏事故污染物CS-137的半衰期长达30.07 a[4]。泄漏后变化复杂表现为危化品种类繁多且多样的行为方式造成危化品泄漏入海后物理和化学形态变化十分复杂。

海上危化品泄漏事故中溢油事故最为频繁，因此英、美等发达国家很早就开展了海上溢油事故相关理论研究，并通过大量的试验验证，建立了众多预测模型[5-8]。与溢油相比，其他类型的海上危化品泄漏漂移扩散数值模拟研究相对较少[9-12]。本文将对近年来国内外海上危化品泄漏在海洋水体环境中漂移扩散数值模拟方面的研究进行整理、分类和介绍。

2 危化品泄漏数值模拟水动力模式

水动力数值模拟是研究海上危化品泄漏漂移扩散的基础。三维模型由于能够全面而立体地反映实际河口海岸和近海环境中海水的运动状况和基本规律，因此在海洋水动力数值模拟研究中得到了广泛研究和应用。伴随着计算机的飞速发展以及数值模拟技术的不断完善，目前三维模型已经成为水动力数值模拟的主要发展方向[13-14]，其中研究应用较多的模型有POM（Princeton Ocean Model)、FVCOM（Finite-Volume Coastal Ocean Model）和ROMS（Regional Ocean Modeling System）。

2.1 POM

POM是由美国普林斯顿大学的Blumberg和Mellor于1977年建立起来的一个三维斜压原始方程数值海洋模式[15-17]。POM采用的蛙跳有限差分格式使得模式具有很高的垂向分辨率[18-20]。POM采用时间分裂算法，模式的外模方程是二维的，内模方程是三维的，内外模分离技术比完全三维计算节省很大计算量[18-20]。POM在水平方向采用正交曲线网格，变量空间配置使用“Arakawac C”网格，能更好地拟合岸线边界，减少“锯齿”效应；在垂向上采用σ坐标变换，可体现不规则海底地形的变化特点；通过干湿网格动边界技术，不仅能很好地处理复杂地形水域的模拟问题，而且可更好地解决三维水动力环境中大量浅滩的“干出”与“淹没”等难点问题[15，19-20]。于海亮[18]以POM为基础并基于Lagrange粒子追踪方法开发了模拟海上溢油的三维输运模型，成功模拟了渤海海域1990年的一次溢油事故。李连峰[19]在POM模式的基础上建立了三维海上类油型化学品及沉降型化学品输运模型，该模型采用Langevin方程与Fokker-Planck方程相结合来模拟油粒子的输运过程，有效地解决了具有随机性的粒子运动过程。通过在宁波-舟山港海域的应用和模拟案例的计算校验，该模型显示出良好的效果。

2.2 FVCOM

FVCOM是由陈长胜（美国麻州大学海洋科技学院）于2000年成功建立的一套三维、有限体积、自由表面和非结构网格的海洋环流与生态数值模型[20-24]。FVCOM在水平方向上采用非结构化三角形网格，能够更精确地拟合复杂曲率的岸界。FVCOM借鉴了许多POM模型的优点，例如：垂直方向采用σ坐标用于拟合复杂的海底地形；使用干湿网格判别法处理近岸滩涂演变问题；采用模式分裂法求解方程大大提高计算效率[22，25]。国内学者藏士文[22]、徐国怀[25]在模拟海上危化品泄漏漂移扩散的研究中分别采用FVCOM建立了渤海和宁波-舟山海域的水动力模型，将实际观测资料与模拟结果进行对比验证，发现潮位、流速和流向的模拟结果与实际观测情况吻合性较高、趋势一致且误差范围较小。

2.3 ROMS

ROMS是一个三维自由表面非线性原始方程近海区域模式，是在垂向静压近似和Boussinesq假定下求解自由表面下Reynolds平均的Navier-Stokes方程[14，26-28]。ROMS水平方向采用正交曲线网格,垂直方向采用σ坐标系统[14，27]。为了提高计算效率节约计算时间，ROMS使用经过正压（快）和斜压（慢）模式之间的特殊处理和耦合的分离显式时间步方案求解动量流体静力原始方程[28-29]。ROMS包含非线性计算内核算法等多种算法，并有MY-2.5混合参数方案等多种垂直混合参数化方案供选择[28-30]。由于ROMS具有集合预报功能，可与多种海洋及气象模式进行耦合使用，因此，近年来被广泛应用于近海海洋环流、海洋生态环境以及海冰等领域的模拟研究[27，29]。崔可夫[14]基于ROMS模式建立了渤海湾三维潮流模型，数值模拟结果在水位、流速和流向方面均与实际观测资料吻合良好，为渤海湾水体交换和污染物迁移扩散研究提供了可靠的水动力参数。陈超[27]运用ROMS模型，对大连湾危化品输运过程进行精细化模拟，模拟结果和实测值吻合较好。

3 海上危化品漂移扩散数值模型

在海洋水动力数值模式基础上，通过引入危化品运动扩散方程，发展形成了适用于不同危化品类型的海上危化品漂移扩散预测数值模型。海上危险化学品一旦发生泄漏，危化品将在海流和风的作用下进行平流、扩散、挥发、溶解和沉降（或上升）等理化过程[31]。海上危险化学品泄漏后在海中的输移和扩散模式主要受两方面影响。一方面是外部环境，如地理环境、水文因素、气象因素和其他因素[16]。危化品在海流和风的影响下将进行平流扩散和湍流扩散：平流扩散指海水的整体运动引起污染物有规律的平移；湍流扩散指因风和潮流引起的随机无序的流体运动,运动状态只能通过概率统计方法加以描述。扩散类型由海流类型决定,并受潮流、风生环流、密度流以及水深和岸线状况等诸多因素制约[32-33]。在海洋这样一个大环境中,污染物一旦泄漏，不仅是单一的平流扩散或湍流扩散，而是两种扩散形式同时作用于污染扩散的整个过程，只不过作用的程度不同[33]。另一方面是化学品的主要物理化学性质，其中最主要的3个因素是溶解度、挥发性和密度。溶解度是决定危化品扩散运动形式的主要因素，是对其他性质进行分析的前提；挥发性决定了危化品入海后是否以蒸气形式向空气中扩散，危化品的挥发性用20℃饱和蒸气压衡量，通常以3 kPa作为挥发与不挥发的界限；密度（与水的相对密度）则决定危化品入海后是在水面运动还是进入水体甚至沉降至海底[34]。根据危化品泄漏入海后的运动形式，我们将海上危化品分为溶解扩散型、海面漂移型、沉降型和挥发型4种类型。

3.1 溶解扩散型

溶解扩散型危化品溶于海水，并在海水中以三维形式进行扩散。对于溶解型危化品在海水中的迁移扩散研究，目前主要从数值优化求解污染物对流扩散方程这一方向来解决。国内外研究求解对流扩散方程的数值方法大致可分两类：一类是欧拉法，将研究海域的空间点作为研究对象，利用解析解的方法或者数值求解的方法，求解危化品在该海域的连续性方程、运动方程和输运方程，从而得出危化品浓度随时间的变化特征；另一类是拉格朗日法，该方法将研究海域中的危化品概化为具有一定数量和质量的质点颗粒，追踪每个质点粒子在研究海域流场中的运动轨迹，得到每个时刻和每个质点颗粒所处的空间位置[35]。

对流扩散方程为：

式中：C为危化品浓度；ux、uy和uz分别为X、Y和Z方向的流速；Dx、Dy和Dz分别为X、Y和Z方向的扩散系数；S0为危化品的源和汇。

在水平运动尺度远远大于垂向尺度和危化品垂向混合比较均匀时，为求解方便往往忽略垂向浓度变化，将对流扩散方程简化为二维模型，以此提高计算效率。杜海涛[33]在对大连湾溶解、轻质和保守液体化学品溢漏后归宿行为的研究中使用了简化后的二维模型，得出的模拟结果符合大连湾实际情况。伴随着计算机性能的提高，越来越多的学者使用三维模型进行研究，于海亮[15]在研究海上液体危化品泄漏三维污染扩散的过程中建立了σ坐标系下的溶解型危化品三维输运模型，并对模型中的对流项分别采用中心差分和一阶Smolarkiewicz迎风格式离散进行优化求解，结果表明一阶Smolarkiewicz迎风格式和改进的二维有限体积法相比中心差分格式具有一定的优势。

3.2 海面漂移型（类油型）

海面漂移型危化品不溶于海水，在海面主要以二维形式漂移扩散。人们对于这类危化品的研究主要集中在海上溢油，因此海面漂移型又可简称“类油型”。海上溢油的漂移扩展过程主要分惯性力扩展、粘性力扩展和表面张力扩展3个阶段进行研究，早期具有代表性的是Fay公式[15-16，36-37]。

惯性力扩展阶段：

粘性力扩展阶段：

表面张力扩展阶段：

式中：D为油膜扩散直径（单位：m）；g为重力加速度（单位：m/s2）；V为溢油总体积（单位：m3）；t为从泄漏开始计算的时间（单位：s）；β=1-ρ0/ρw，ρ0和ρw分别为危化品密度和海水密度（单位：103kg/m3）；Vw为海水运动粘性系数（单位：m2/s）；δ为扩展系数；K1、K2和K3分别为3个阶段的比例系数。

泄漏后的油品在海水中形成油膜并在海水运动的作用下产生漂移，同时油膜扩散范围逐渐扩大,油膜漂移距离的长短通常用油膜等效圆中心位移来判断[15]。其位移s为：

Fay理论及其修正模型对于在开阔海域瞬时泄漏情况下的类油型污染物模拟取得了很好的结果，适合类油型危化品在泄漏初期阶段的模拟预测[38]。连续泄漏油膜在距离泄漏源一定距离后，侧向扩散起主导作用，而对于水下泄漏，浮力的作用至关重要，由于Fay模式忽略这些作用，所以它并不适用于连续泄漏和水下泄漏的模拟[39]。

Johanseen等[39-43]提出的“油粒子”模型是当今世界主流的溢油模式。该模型把溢油分散成有限个油粒子来模拟实际溢油的漂移扩散过程。“油粒子”模型不仅解决了溢油在重力扩展停止后的扩散现象，还实现了海上油膜破碎分离现象的模拟，对溢油在海洋环境中受到的海洋动力因素的影响模拟得更加准确[42-44]。国内学者Yang等[45]基于国家海洋环境预报中心自主研发的溢油预报模型，采用Lagrange追踪法模拟油粒子的三维时空运动，建立了渤海湾三维溢油预报模型，较好地重现了2011年6—8月渤海蓬莱19-3溢油事故的发展过程。Li等[46]基于“油粒子”模型重现了“桑吉”轮事故溢油的漂移扩散过程，模拟结果与实际观测相吻合。虽然“油粒子”模型代表了当今溢油模拟的方向，但是其缺点也很明显：所需的模拟粒子数过于庞大，导致计算量非常大；忽略了溢油自身的重力扩展作用。黄娟等[47]和刘伟峰等[43]研究发现在溢油大规模泄漏初期，由于油膜面积在短时间内急剧扩大，溢油自身的扩展效应显著大于湍流扩散效应，如果仅通过水平扩散的方式来模拟油膜扩展过程，将导致“油粒子”扩散速度响应过慢，不能真实反映溢油量对扩散面积的影响。

3.3 沉降型

沉降型化学品主要指密度比海水大且不溶或微溶的化学品，入水后分散成小液滴在浪、潮和流的作用下随水体运动，其运动形式基本可分为扩散和漂移、沉降、沉积和再悬浮[16，48-49]。

海洋水动力要素是影响危化品在海水中扩散和漂移过程的主要因素。扩散是指危化品泄漏入海后由于浓度差和湍流等物理因素而形成污染物向整个水体迁移的混合过程；漂移指泄漏入海的危化品在风、表层和次表层流以及波浪作用下的拉格朗日漂移过程，漂移运动主要由平流条件决定[16，48]。陈协明[48]在研究化学品的溢漏扩散模式时指出，沉降型化学品在海面上的漂移主要受风的切应力、表层/次表层流和余流（波生余流和潮余流）控制。

沉降过程一般可以分为3个阶段：沉降射底、触底和底浪的传播[16]。危化品入水后发生破碎并迅速下降，形成了充分成长的高速射流。在这一过程中，周围海水与危化品发生混合，沉降射流的体积不断扩大；沉降射流触底时将形成以碰撞点为中心向四周扩展的高密度底浪；开始时底浪的传播速度很快，随着传播距离的增加，底浪的厚度和传播速度迅速减少，底浪中携带的危化品迅速下降并沉积[16，48]。

危化品沉降到海底后会形成沉积，在流与浪产生的剪切力的共同作用下，沉积的危化品会发生再悬浮，并回到水体中重新进行扩散、漂移和沉降过程[16，48-49]。沉积物的再悬浮过程十分复杂，目前还没有成熟的预测方法，国外学者McNeil等[50-52]使用水槽装置模拟水体底泥的再悬浮过程，国内学者陈协明[48]和李连峰[19]使用经验参数（沉积物的5%）计算再悬浮通量。

3.4 易挥发型

根据挥发性危化品泄漏后混合气云与环境大气的密度差异，将混合气云分为浮性气云（密度比空气小）、中性气云（密度与空气相近）和重质气云（密度比空气大）3类。研究学者往往将浮性气云和中性气云归为一类进行研究[53]，简称为“非重气云”。非重气云团扩散模型中最具代表性的是高斯（Gauss）模型，其理论基础是湍流扩散梯度输送理论，忽略重力对气体扩散的影响，认为扩散主要由空气的湍流决定，污染物在扩散截面的浓度呈正太分布，扩散系数K为常数[5，54]。高斯模型由于所需参数少、运算量小和计算快捷等优点，被广泛应用于小尺度的污染物模拟。重气云团的扩散比中性和浮性气云复杂的多，主要经历闪蒸和空气夹带、密度差作用下的重量沉降、重力沉降的地表作用和被动扩散4个阶段[54]。重气云团的扩散研究不仅要考虑气体重力因素还要考虑重气云流动扩散阻力等因素的影响。孙莉等[6]通过收集大量重气云扩散实验结果加以分析整理，解析出能较好反应重气云瞬时和连续释放规律的方程式，并成功建立了唯象模型（BM模型）。BM模型是最简单的根据经验判断的重气云扩散模型，可以根据重气云扩散曲线图和方程式快速得出相对应的扩散数据，因此被广泛应用。

计算机仿真模拟是研究易挥发型危险化学品泄漏扩散的又一热点领域。计算机仿真模型考虑的因素较多，模拟时间较长，因此预测精确度很高。此类模型最具代表性的就是CFD（Computational Fluid Dynamics）和ALOHA（Areal Locations Of Hazardous Atmospheres）[55]。CFD具有可视化能强和计算方法成熟等特点[43],在实际中被广泛应用。Ikealumba等[56]在应用CFD研究海上液化天然气的泄漏扩散规律时发现，海洋波动带来的不稳定性比风场更容易促进液化天然气的扩散。ALOHA模型包含两种不同的扩散模型：高斯模型和重气云模型，能根据危化品理化性质和泄漏量等参数自动选择扩散模型。ALOHA不仅能够模拟挥发型危化品的扩散过程，还能模拟危化品源强的变化规律及扩散区域的浓度变化规律[5]。王志宪[5]在应用ALOHA模型分析液氯泄漏灾害的影响因素时指出，ALOHA对气体危化品泄漏后的扩散范围和危险浓度阈值的估算预测具有较高的准确度。

4 总结与展望

本文简要概述了前人关于海上危险化学品泄漏扩散数值模拟方面的研究成果，并将危化品分为漂浮型、溶解型、沉降型和挥发型4大类，分别介绍其泄漏后在海洋环境中的扩散形式及归宿特点；通过最广泛的数值模型阐述其理论方法，并对其优缺点和特性进行讨论。

将危化品划分为4种类型是一种十分理想的状态。数值模型中只考虑了最主要的扩散形式，忽略了次要扩散过程，因此可以将复杂的运动过程简化为单一的运动形式，虽然大大减小了计算工作量，但也会严重影响模拟结果的准确性。现实生活中危化品种类繁多，泄漏到海洋环境中的扩散过程可能会同时发生挥发、溶解、沉降和漂浮等复杂的运动过程。为了更真实地模拟海上危化品泄漏扩散的实际情况，未来应将危化品的泄漏扩散过程作为一个整体，同时考虑沉降、溶解和挥发等运动过程，构建相互耦合的数值模型。