韩龙喜张防修

(1.扬州市职业大学 资源与环境工程学院，江苏 扬州 225000； 2.河海大学 环境科学学院，南京 210000)

1 研究背景

长江航道是我国重要的水运物资通道，近年来由于航运繁忙，已发生多起因事故引发的溢油污染事故[1]。事故发生后，如何科学预测溢油输移扩散特征，对于提高应急处置水平具有重要意义。溢油数值模拟是当前该领域的研究热点，现有的研究更多集中于海洋、河口区域[2-6]，而针对狭长型内陆河道的研究相对有限。Sayed等[7]、Goeury[8]在河道溢油事故处理中构建的数值模型简便实用。Toz等[9]、Gautama等[10]建立的河道溢油数值模型综合考虑了水流、随机扩展以及溢油风化等因素的影响，并将遥感技术应用于模型的验证中，取得了较好的应用效果。在国内，刘栋等[11]基于机理性实验分析了潮汐作用下河道中溢油的扩展、漂移特性。赵琰鑫等[12]在长江武汉段采用传统经验模型和油粒子模型相结合的方法分析了溢油迁移特征。

溢油在扩展过程中不接触水陆边界是传统经验模型应用于溢油自身扩展分析的前提条件。这种假设在事发水域较宽敞且溢油以中心排放时较易实现，但是如果事发水域为内陆河道或溢油以岸边排放形式进入水体时则很难实现。针对狭窄水域中发生的溢油事故，韩龙喜等[13]提出的基于油粒子模型分析溢油漂移扩散和扩展过程的研究方法能较准确地反映溢油扩散过程，但是该方法单纯通过水平扩散方式模拟油膜自身扩展，导致油粒子扩散速度响应过慢，且未能有效反映溢油量对扩散面积影响。

本文采用水深平均二维潮流模型模拟评价区域水流流场，采用油粒子模型分析溢油输移扩展。为了能准确反映溢油自身扩展的变化过程，采用随机运动模拟了溢油的惯性扩展、黏性扩展和表面张力扩展过程，使油粒子自身扩展更符合Fay提出的三阶段扩展理论[14]。数值模型同时考虑了溢油风化及触岸吸附等过程。将数值模型应用于假设工况的溢油事故分析中，取得了理想预测效果。该方法可为内陆河道突发性溢油事故预报和应急处理提供技术支持。

2 水动力学模型

根据长江河道特点，本文采用水深平均二维潮流模型[15]模拟研究区水流场，水动力数据获得准确性验证后，可以作为溢油水动力场的驱动条件。水动力模型通过连续方程和动量守恒方程来表征流场和潮位的变化特征。

连续方程为

(1)

动量守恒方程为

(2)

式中：Z，H分别表示水位和水深(m)；u，v分别为x，y向的流速分量(m/s)；νt为紊动粘性系数(m2/s)；g为重力加速度；c为谢才系数；f为柯氏系数，f=2ωsinφ，其中ω为地球自转角速，φ为计算地水域的地理纬度。

上述二维潮流模型基本方程中含有非线性混合算子，采用剖开算子法[16]进行离散求解。

3 油粒子预测模型

油粒子模型[17-18]最早由Johansen和Andunson(1982)提出，该模型假设溢油为大量离散的油粒子聚合体，每个油粒子代表一定油量。模型通过综合分析油粒子在水流、风向以及自身随机游走作用下的运动输移，分析溢油随时间迁移及其空间分布特征。油粒子模型将溢油在Δt时间内的运动分成对流、风导漂移和随机游走3个过程。单个油粒子的运动方程可表示为

Xn+1=Xn+ΔXC+ΔXW+ΔXD。

(3)

式中：Xn+1，Xn分别表示单个油粒子在(n+1)Δt，nΔt时刻的空间位置列向量；ΔXC，ΔXW，ΔXD分别表示该油粒子受表层水流对流运动影响、风应力影响以及水体紊动扩散影响而产生的空间位置变化的列向量。

3.1 溢油随流运动模拟

在Δt时段内，受水流对流作用产生的油粒子空间位移可以用确定性方法模拟，即

(4)

式中Un，Un+1分别表示油粒子在nΔt，(n+1)Δt时刻的水流速。

3.2 溢油风导(应力)漂移模拟

在Δt时段内，受表层风应力直接作用，油粒子产生的漂移距离ΔXW可以表示为

ΔXW=αDW10Δt。

(5)

式中：α为风漂移因子，取值范围为0.03～0.04；D为考虑风向偏转角的转换矩阵；W10为水体表面10 m高处的风速向量。

3.3 溢油随机游走模拟

受随机游走影响，油粒子云团的尺度和形状随时间不断变化。本文基于随机步长法[19]计算油粒子在水平方向上产生的随机运动距离ΔXD，即

(6)

式中：Rn表示[-1，1]之间的均匀分布随机数；θ′表示位于[0，π]之间的均匀分布的随机角度；De是由溢油自身扩散引起的扩散系数；DT是油粒子作布朗运动的扩散系数。

DT的计算式可表示为

(7)

式中：nb是曼宁系数；Vc为水流速；h为水深。

采用油粒子模型模拟溢油机械扩展时，存在溢油运动响应慢，且未能反映溢油量对溢油扩展影响的不足。鉴于Fay提出的溢油扩展三阶段理论能较真实地反映溢油机械扩展过程[18]。本文采用Fay公式推导溢油自身扩散系数De，如式(8)所示。该方法的实质是通过溢油随机运动反映溢油惯性扩展、黏性扩展和表面张力扩展过程。

式中：常数K2i=1.14，K2v=1.45，K2t=2.30；Vp是油粒子的体积(m3)；σn为张力系数；ρw，γw分别为水的密度和运动黏性系数。

3.4 溢油风化过程

溢油风化包括蒸发、乳化及溶解等过程。风化过程对溢油组成产生影响，但对溢油水平位置不发生影响[20]。

3.4.1 蒸 发

在溢油事故初期，溢油质量传输主要通过蒸发完成。溢油蒸发率与溢油本身物化特征有关，通常通过溢油蒸发率表示，即

(9)

3.4.2 乳 化

溢油乳化表现为油包水的过程。溢油乳化进程受水域条件影响，但最终乳化量仅取决于溢油中乳化剂的含量。本文采用含水率来表征溢油乳化程度，即

(10)

3.4.3 溶 解

溶解是指油粒子在外界能量的扰动下，均匀进入水体的过程，溶解量和溶解速度主要受溢油自身物化性质、风速和水况条件的影响。通常由溢油溶解率反映溢油溶解程度，即

(11)

3.5 河岸吸附

溢油沿水流进入长江岸边区域时，由于水流速往往与水边线存在一定夹角。溢油粒子极易触岸，根据Shen等[21]的研究结论，触岸油粒子可能会出现被岸边完全吸收、部分吸收、完全反射3种可能。本文作简化分析，只考虑完全吸收和完全反射2种情况。设定一个黏附概率，某油粒子触岸时会获得一个介于0～1之间的随机数，当该随机数小于设定黏附概率时，油粒子被完全吸收，该油粒子坐标被修正为计算网格中距离靠岸点最近的节点坐标。

3.6 油膜厚度计算

油膜厚度变化是溢油事故应急处理中一个重要分析指标。在掌握油粒子空间分布规律后，可通过分析一定水体内油粒子数量、质量、体积等参数计算得到油膜厚度分布。

假设水中油粒子个体大小均等，A表示水面面积，N表示水面油膜中包含的油粒子个数，m表示受风化作用后单个油粒子质量。按液体厚度计算公式，油膜在t时刻的厚度ht可表示为

(12)

式中ρ为溢油密度。

4 溢油事故数值模拟

选择镇江扬州段河道中仪征码头上游12 km至下游13 km范围，作为溢油事故模拟的研究区域(119°13′E—119°08′E，32°22′N—32°13′N)。该江段包含了仪征市水厂的重要取水口，周边分布有码头、润扬大桥、江心洲、边滩等，航道条件较为复杂(如图1)。近年来长江航运业日益繁忙，过往船舶在水上运输或码头作业过程中一旦发生溢油事故，将对该地区工农业生产和群众饮水安全产生严重影响。

图1 模拟江段位置示意图Fig.1 Position of the simulated reach of the Yangtze River

溢油事故产生的影响受水文、气象等多条件综合作用。本文根据水环境敏感目标位置、水文水动力条件以及污染源位置的代表性确定6种计算方案，具体信息见表1。这里按照考虑最不利影响原则，选择方案3(溢油时刻为大潮涨潮时，风向取ESE方向，风速取平均风速2.6 m/s，化学品本底浓度取0 mg/L)具体分析说明。

表1 溢油事故风险预测方案Table 1 Schedules of risk prediction for oil spill accidents

4.1 模拟事故工况

根据以往研究成果[22]确定长江镇扬段河道的糙率系数，长江主槽糙率一般取值0.018～0.022，河道滩地糙率一般取值0.024～0.028。分散系数选用Ex=αxhu确定，αx取为4.0，αy取为0.5。基于GIS平台，采用矩形网格划分水动力模拟计算区域，平面共布置360(纵向)×300(横向)个节点(网格)，网格边长50 m。依据1∶10 000水下地形等值线图，计算提取各节点河底高程。考虑到待评价江段与大通站间支流入流量相对较小，这里以大通站最小月平均流量作为一维水流模拟的上边界条件；用同期的下游潮位站潮位过程作为下边界条件[23]。经一维水动力学数学模型模拟后得到评价区域二维水动力学模拟的上、下游边界水文要素变化过程，并以此作为设计潮流量、潮位边界条件，模拟设计潮流过程的水动力特征。

模拟事故情景为某油船在仪征拟建码头处发生原油泄露，溢油源强取为140 t，泄漏时间为10 min，溢油油品为轻质中原原油，其理化性质如表2所示。

表2 油样主要理化性质Table 2 Main physical and chemical properties of oil samples

注：ρ0为油品密度；μoil为原油黏度

4.2 结果分析

采用二维非稳态水动力模型模拟分析研究区水动力流场特征。水动力模型通过多个测站的潮位潮流资料进行验证。大潮涨急时刻流场见图2。采用Willmott[24]的统计学方法对水动力模型计算结果进行效率评价，结果显示效率评价值为极好。水动力模型能较好地呈现研究区复杂水流特征。

图2 大潮涨急时刻流场Fig.2 Flow field at spring tide

如图3所示，根据模拟计算结果，溢油在大潮涨潮开始时进入长江，在潮流和风力的共同作用下，油粒子大致沿西方向沿上岸向上游漂移。仪征取水口离事发位置较近，溢油事故发生后对取水口的影响如表3所示，约20 min后油膜到达仪征水源地准保护区下边界，油膜在事故发生后150 min离开水源地保护区，对仪征水源地保护区的持续影响时间为130 min，对取水口的持续影响时间为5 min。水源地保护区内油膜最大厚度范围为2.67～42.36 mm，其中取水口附近油膜最大厚度为6.77 mm。在未采取有效应急措施的情况下，溢油事故对该水源地水质造成直接影响。因此，应及时启动水源地应急监测计划，确保供水安全。

图3 油粒子在事故发生后各时刻运动轨迹Fig.3 Trajectories of oil particles at differentinstances after an accident

位置影响目标油粒子中心到达时间/min溢油持续影响时间/min折算油膜最大厚度/mm仪征取水口准保护区下边界202042.36二级保护区下边界401014.60一级保护区下边界50108.45仪征取水口6056.77一级保护区上边界80156.46二级保护区上边界130505.21准保护区上边界150202.67

5 结 论

溢油事故是长江航运过程中的易发事故类型。本文针对长江河道溢油事故特点，采用二维潮流模型模拟评价区域水流流场。在溢油迁移扩散分析中，采用改进型油粒子模型分析溢油迁移扩散特征，使得模型较准确地反映了溢油迁移和扩展过程。该模型同时考虑了溢油风化及岸边吸附对溢油迁移的影响。将建立的溢油事故分析模型应用于镇江—扬州江段中虚拟事故的分析，结果表明在风向取ESE方向、大潮涨潮时发生溢油污染事故，油膜大致沿西方向沿上岸向上游漂移。该模型不仅可以科学预测油膜对仪征水源地保护区和取水口的持续影响时间，以及油膜厚度变化，对于分析溢油迁移扩散和厚度变化也具有较好的指导意义。

目前该主题的研究更多关注溢油在二维平面的迁移变化，在后续研究中，可以尝试开展溢油在纵向立面的三维分析。在研究方法上，可以尝试将水动力模型、水质模型及GIS动态可视化技术集成，以提高分析结果的呈现效果。