长江口宝钢码头溢油事故油膜漂移扩散数值模拟

顾靖华, 朱建荣, 金智

河口海岸学国家重点实验室(华东师范大学), 上海 200241

码头、港口、航道和桥墩等邻近水域是溢油事故多发之处。溢油事故发生后, 油膜的漂移扩散会对周围水域和敏感目标造成污染。因此, 研究溢油事故后的油膜漂移扩散, 可为溢油事故的处理提供油膜面积、到达敏感目标时间和油膜厚度变化等数据, 对于事故应对和环境保护具有实际意义。

研究溢油事故发生后油膜运动的主要方法有现场实验、水槽试验、理论分析和数值模拟。溢油事故后油膜漂浮于水面, 其运动主要是由水流的平流,以及风引起的漂移和随机扩散决定, 还受蒸发、分散、乳化、沉降及生物降解等过程的影响。张和庆 等(2001)在珠江口通过现场溢油实验, 分析了油类入海后漂流和扩散过程。郭运武 等(2008)通过水槽实验, 分析了在不同风速和水流速度条件下油膜的扩展、漂移运动状况, 结果表明风在溢油油膜扩展、漂移中起着重要的作用。刘栋 等(2006)通过潮汐水槽模拟实验, 观察和记录油膜在潮汐水流中两个半日潮潮周期内的扩展、漂移情况, 分析了溢油油膜在潮汐水流中的扩展、漂移过程。娄厦 等(2008)理论分析了油在水中的行为和归宿, 概括了多个溢油模型的理论依据、基本假设和应用范围。Violeau 等(2007)通过室内试验和应用平滑粒子水动力拉格朗日方法模拟溢油的漂移扩散。在溢油运动的研究和预测中, 数值模拟是一种应用最为广泛的研究方法。较早的一些模式应用垂向平均的二维模式较多,例如娄安刚 等(2001)建立了胶州湾内的二维预报潮流场, 再考虑风场对海面溢油的影响, 模拟了溢油事故后溢油的油膜漂移轨迹。孙长青 等(2003)建立渤海湾二维溢油飘移数值模型, 计算在不同风场和潮流场的共同作用下, 溢油飘移的输移轨迹、可能抵岸的时间、地点及残留量等。李大鸣 等(2008)建立了在潮汐、风浪作用下油膜运动二维数学模型,对大亚湾海区油码头瞬时溢油事故后的影响范围进行了模拟。Periáñez 等(2008)建立了二维水动力模式,模拟直布罗陀海峡潮流和余流, 并质点跟踪技术模拟潮流、混合和风对溢油运动的影响。垂向平均的二维数值模式不能真实反映风生流和密度流等在垂向上的变化, 从而影响模拟油膜漂移扩散的精度。因此, 将三维模式运用到油膜漂移扩散的模拟中,能较为真实地模拟风生流和密度流等的垂向变化,并更为准确地模拟表层流场对溢油漂移扩散的作用。例如, Wang 等(2008)采用粒子方法在二维溢油模式的基础上发展三维溢油模式, 应用到渤海溢油事故的数值模拟。应用三维数值模式, 陈义中 等(2006)模拟了渤海原油码头溢油漂移扩散, 傅利辉等(2008)模拟了长江口北支崇启大桥建成后流场变化及溢油漂移扩散, 陈士谦 等(2012)模拟了鸭绿江公路大桥溢油的漂移扩散。

长江口为三级分汊四口入海的大型潮汐河口,环境敏感保护目标众多(表1)。环境敏感保护目标位于溢油事故发生点的不同方向, 不同的方向对应不同的敏感目标。数值预测油膜运动需要考虑主导风和不利风, 其对油膜的运动起着十分重要的作用。对溢油事故附近气象站实测风作统计分析, 可得出各个方向的频率和平均风速, 第一频率的风取为主导风。不利风向取决于溢油事故发生点和敏感目标的连线, 与指向敏感目标最接近的风向为不利风向。主导风和不利风可从实测风的统计结果给出。本文采用崇明东滩气象站实测的风资料(测站见图1), 给出主导风和不利风向。在此基础上, 应用长期改进和应用的河口海岸水动力三维数值模式

表1 长江口环境敏感保护目标Tab. 1 Environmental sensitive protection targets in the Changjiang Estuary

ECOM-si(semi-implicit estuarine, costal and ocean model), 耦合溢油漂移扩散模块, 模拟并分析上海港宝钢原料区域码头船舶发生溢油事故后油膜的漂移扩散(溢油事故点位置见图1), 以及对环境敏感保护目标(图2)的影响。

图1 长江口形-势图及上海港宝钢原料区域码头溢油事故泄漏点Fig. 1 Sketch map of the Changjiang Estuary, and the site of oil leakage accident at the Baosteel raw material regional wharf of the Shanghai Port. The green dot marks the weather station at the Chongming eastern shoal; the red dot represents the oil leakage point, black dots stand for the water intakes of the reservoirs

图2 长江口环境敏感保护目标位置Fig. 2 Positions of the environmental sensitive protection targets in the Changjiang Estuary

1 数值模式设置

采用作者课题组长期改进和应用的河口海岸水动力三维数值模式ECOM-si, 该模式在长江口水动力过程和盐水入侵等方面研究中经过了大量的率定和验证, 得到了广泛的应用, 取得了大量研究成果

(Lyu et al, 2018; Chen et al, 2019; Zhu et al, 2020; Ma et al, 2022) 。在ECOM-si 的基础上, 耦合油膜漂移扩散模块。该模块将油分成有限个粒子, 每个粒子的位置按以下公式计算:

式中X0、Y0为某粒子的上一时刻坐标;X、Y为某质点的下一时刻坐标;U、V分别为X、Y方向的流速分量, 由水动力数值模式计算得出;W10为海面上10m 处的风速;A为风向;α为风漂移系数。在风的影响下, 油膜漂移速度的增加量约为风速的2%～3%,漂移方向与风向成0°～40°夹角。由于油膜中心的漂移速度和方向是表面水流和风所引起的流速矢量之和, 在模拟中使用风漂移系数α为风与海流的角度修正系数, 一般在0.02～0.03 之间, 本文取0.03;r为随机扩散项,r=RE,R为0～1 之间的随机数,E为扩散系数, 计算公式为E=, 其中Aα为水平扩散系数, 由水动力模型提供;B为随机扩散方向, 计算公式为B=2πR。

在油膜的漂移扩散过程中, 还需考虑油膜的蒸发。油膜蒸发受溢油种类、气温、溢油面积、风速、太阳辐射和油膜厚度等因素的影响。本模块考虑溢油的蒸发公式为(Stiver et al,1984)

其中,Fv为溢油蒸发体积分数;T是环境温度(空气温度), 本文长江口冬季取283K;θ为无因次时间(蒸发暴露时间);A1、B1、T0、TG为各种依赖于油品的常数,A1=6.3,B1=10.3;T0为油在Fv=0 时的初始沸点温度(单位: K);TG为蒸馏曲线的坡度(单位: K);T0、TG根据下列经验公式计算:

其中API 为比重指数, 对重油、重油和轻油可分别取值为16.9162、26.8458 和36.1882。蒸发公式中的θ由下式计算:

其中K是空气端传质系数,A2是溢油面积,t是时间,V是溢油泄漏的体积。针对公式实际使用中蒸发过快的现象, Fingas(1995)对K值进行了调整, 为了解决无风速情况下传质系数公式无法计算的问题, 孙宝楠 等(2010)假设在传质系数上加一个小量R, 从而传质系数可以表示为

其中U10为海表面 10m 风速。R1的范围在(0.1～10)×10–4m·s–1的量级(本文K值取10–4),t为溢油发生后的时间, 单位为h。每个油粒子初始时刻的体积(总溢油体积/划分总粒子数)减去t时刻的溢油蒸发量, 即为t时刻每个油粒子的体积。

溢油事故发生点为上海港宝钢原料区域码头石灰石及砂岩卸船码头, 坐标31°28′6″N, 121°25′33″E,泄漏量1000t, 0#轻质柴油。在径流作用下, 油膜向海方向输运, 落憩时刻之后为涨潮, 相对于涨憩时刻,油膜会有更长时间滞留在口门内, 为不利情况, 故假设溢油事故发生在大潮期间的落憩时刻。水动力模式先运行10d, 水位和流场准稳定后, 再释放油粒子。只要质点数足够多, 油膜漂移的模拟就越精确, 因此本文设置了一千万个油膜质点。每个质点代表一个恒定的体积大小, 每步根据网格的水平面积大小和所含质点数, 可以求出油膜在该网格的厚度值: 厚度=质点数×t时刻单个质点代表的体积/网格的水平面积。

模型计算范围包括整个长江河口、杭州湾和邻近海域(图3a), 东至125°E 附近, 北至接近34°N, 南至27°30′N 左右。模型水平方向上采用曲线非正交网格, 较好地拟合了长江河口的岸线, 并对南北支分汊口、深水航道等区域加密(图3b、c)。长江口内网格分辨率最高可至约100m, 长江口外网格分辨率最高可至10km 左右。模型在垂直方向上采用σ坐标, 由表至底均匀分成10 层。计算时间步长取为60s。由于长江河口潮滩分布较广, 模型采用干湿网格判断实现潮滩动边界, 最小临界水深阈值设置为0.2m。岸线和水深数据为2017 年长江河口实测资料,水深数据以黄海85 高程为基面, 见图1。

图3 模型计算范围和网格(a)以及局部放大的南北支分汊口附近网格(b)和南支口门附近网格(c)Fig. 3 Model calculation domain and grids (a), and the enlarged view of grids near the bifurcation between the South Branch and North Branch (b) and near the mouth of the South Branch. The red lines are the guide dikes and groins of the Deepwater Channel

本文模拟和分析冬季1 月径流量和主导风、不利风作用下油膜的漂移扩散。选取大通作为上游边界, 模式的径流边界条件可直接采用大通水文站的月平均径流量资料。1950—2020 年多年月平均长江入海径流量1 月为11700m3·s–1。外海开边界潮位考虑由16 个分潮(M2, S2, N2, K2, K1, O1, P1, Q1, MU2,NU2, T2, L2, 2N2, J1, M1和OO1) 的调和常数合成得到, 资料来自全球潮汐模式NAOTIDE(http://www.miz.nao.ac.jp)计算结果。

对主导风和不利风向, 本文采用长江口崇明东滩气象站(测站位置见图1)2007 年至2021 年14 年的实测风速风向资料, 统计1 月8 个罗盘方向风的频率和平均风速。1 月第一至第八风向频率分别为26.9%、26.1%、16.5%、8.3%、7.7%、6.3%、5.4%、2.8%, 风向分别是N、NW、NE、E、SE、W、S、SW, 风速分别是5.6、4.8、5.0、4.1、4.1、3.6、4.0、3.4m·s–1。图4 为1 月风玫瑰图, 形象给出了8 个方向风的频率和平均风速, 以及各个方向上不同风速(间隔2m·s–1)及其对应的频率。

图4 崇明东滩气象站2007—2021 年统计得出的1 月风玫瑰图Fig. 4 Wind roses in January calculated from 2007 to 2021 at weather station on Chongming eastern shoal

主导风取为第一频率的风。1 月第一频率(26.9%)风向为北风, 风速5.6m·s–1。不利风由溢油事故发生点和环境敏感保护目标位置确定。由于长江口环境敏感保护目标众多(图2), 需要给出多个不利风向。上海港宝钢原料区域码头船舶发生溢油事故发生点的上游敏感目标为陈行水库、太仓水库、东风西沙水库、刀鲚核心区和实验区, 不利风向为东南风, 为第5 频率(7.7%), 风速4.0m·s–1。对下游东南向敏感目标为青草沙水库、刀鲚实验区、吴淞湿地和九段沙湿地, 不利风向为西北风, 为第2 频率(26.1%),风速4.8m·s–1。对下游东北向敏感目标为青草沙水库、刀鲚实验区、崇明东滩保护区和中华鲟保护区,不利风向为西南风, 为第 8 频率(2.8%), 风速3.2m·s–1。

2 油膜漂移扩散的数值模拟

本文模拟和分析上海港宝钢原料区域码头综合改造工程码头溢油事故发生后油膜的输运扩散及其对敏感目标的影响, 对模拟结果绘制油膜分布和水源地取水口油膜厚度随时间变化过程图进行分析。

2.1 主导风北风5.6m·s–1 风速

油膜的运动主要受水流的平流、风的漂移和随机扩散作用。本文假设宝钢码头的溢油事故发生在大潮落憩时刻, 长江河口为半日潮, 事故发生后3h内为涨潮流, 在潮流的平流作用下向上游运动(图5)。油膜中心位于陈行水库水源地保护区, 厚度超过2μm。之后在径流的作用下, 油膜在随涨落潮流上下移动过程中向下游输运, 至12h 油膜已进入南港,中心位于外高桥附近(图6)。至72h, 油膜继续向下游漂移扩散, 大部分位于南港中下段和南槽中上段的南侧, 小部分进入北槽南侧(图7)。在油膜的运动过程中, 绝大部分位于长江口南岸, 在主导风北风的作用下向南漂移起着重要的作用。

图5 溢油事故后3 小时油膜分布Fig. 5 Distribution of oil film thickness after 3 hours of the oil spill accident

图6 溢油事故后24 小时油膜分布Fig. 6 Distribution of oil film thickness after 24 hours of the oil spill accident

图7 溢油事故后72 小时油膜分布Fig. 7 Distribution of oil film thickness after 72 hours of the oil spill accident

从长江口水库取水口油膜厚度随时间变化过程看(图 8), 因受径流和北风的作用, 溢油事故对南支上端北侧的东风西沙水库取水口无影响。太仓水库位置南支中段南侧, 离溢油事故发生点较远, 同样受径流和北风的作用, 油膜最大厚度小于 1μm, 影响微小。青草沙水库取水口位于溢油点东北侧不远处, 受北风作用, 油膜最大厚度小于 1μm, 影响微小。陈行水库取水口位于溢油点上游不远处, 油膜2h 随涨潮流到达取水口, 之后随着潮流周期性起伏, 最大厚度约为 68μm。31h后, 油膜的影响消除, 期间最长持续影响10h, 总影响时间为29h。

图8 油膜厚度在东风西沙水库取水口(a)、太仓水库取水口(b)、陈行水库取水口(c)和青草沙水库取水口(d)随时间变化情况黑线: 主导风(北风); 红线: 不利风向(东南风); 绿线: 不利风向(西北风); 紫线: 不利风向(西南风)Fig. 8 Temporal variation of oil film thickness at the water intakes of the Dongfengxisha reservoir (a), Taicang reservoir (b),Chenhang reservoir (c) and Qingcaosha reservoir. Black line: prevailing wind (northerly wind); red line: unfavorable wind direction (southeasterly wind); green line: unfavorable wind direction (northwesterly wind); purple line: unfavorable wind direction (southwesterly wind)

对青草沙饮用水水源保护区, 溢油事故发生后11h 油膜到达, 36h 后离开, 持续影响25h, 期间最大厚度为1.42μm。对陈行饮用水水源保护区, 溢油事故发生后 1.5h 油膜到达, 53h 后离开, 持续影响51.5h, 期间最大厚度为151.31μm。对长江太仓浏河饮用水水源保护区, 溢油事故发生后3.5h 油膜到达,5.5h 后离开, 持续影响 2h, 期间最大厚度为2.99μm。溢油事故对东风西沙饮用水水源保护区都无影响。

对长江刀鲚国家级水产种质资源核心保护区,溢油事故发生后2.5h 油膜到达, 29h 后离开, 持续影响26.5h, 期间最大厚度为47.84μm。对长江刀鲚国家级水产种质资源实验保护区, 溢油事故发生后0.5h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为52.99μm。对上海吴淞炮台湾国家湿地公园, 溢油事故发生后8h 油膜到达, 持续影响时间超过第72h, 期间最大厚度为32.70μm。对上海九段沙湿地国家级自然保护区, 溢油事故发生后44.5h 油膜到达, 影响时间超过溢油发生后的72h,期间最大厚度为1.47μm。溢油事故对上海崇明东滩鸟类国家级自然保护区、上海市长江口中华鲟省级自然保护区无影响, 缘于它们远离溢油点, 位于东北的北港口门北侧。

2.2 不利风向东南风4.0m·s–1 风速

在不利风向东南风4.0m·s–1作用下, 油膜朝西北方向漂移, 聚集于南支北岸, 同时受径流作用向下游输运,在溢油事故后72h 油膜位于南支的上下段的中北侧和北港上段的中北侧(图9)。这与主导风北风作用下油膜的分布十分不同, 表明风在油膜漂移扩散中起着十分重要的作用, 迫使油膜克服径流朝风向的方向漂移。

图9 溢油事故后72h 油膜分布Fig. 9 Distribution of oil film thickness at 72 hours after the oil spill accident

从长江口水库取水口油膜厚度随时间变化过程看(图8), 对东风西沙水库取水口, 溢油事故发生后16h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为10μm。与主导风北风结果相比, 东南风使得油膜向西北漂移更远的距离。对太仓水库取水口, 溢油事故发生后4h 油膜到达, 8h 后离开, 持续影响4h, 期间最大厚度为1.8μm。对陈行水库取水口, 溢油事故发生后2h 油膜到达, 12h 后离开, 持续影响10h, 期间最大厚度为7μm。对青草沙水库取水口, 溢油事故发生后10h 油膜到达, 39h 后离开,持续影响29h, 期间最大厚度为2.3μm。

对青草沙饮用水水源保护区, 溢油事故发生后9h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为10.45μm。对陈行饮用水水源保护区,溢油事故发生后1h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为149.46μm。对长江太仓浏河饮用水水源保护区, 溢油事故发生后3h 油膜到达, 55.5h 后离开, 持续影响52.5h, 期间最大厚度为52.27μm。对东风西沙饮用水水源保护区, 溢油事故发生后17.5h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为37.89μm。

对长江刀鲚国家级水产种质资源核心保护区,溢油事故发生后2h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为84.47μm。对长江刀鲚国家级水产种质资源实验保护区, 溢油事故发生后0.5h 油膜到达, 影响时间超过溢油发生后的72h, 期间最大厚度为110.00μm。溢油事故对上海吴淞炮台湾国家湿地公园、上海崇明东滩鸟类国家级自然保护区、上海九段沙湿地国家级自然保护区、上海市长江口中华鲟省级自然保护区均无影响。这与油膜的平面分布图是一致的。

2.3 不利风向西北风4.8m·s–1 风速

溢油事故后72h 油膜分布与主导风北风作用下的分布相似, 这里不再给出。

从长江口水库取水口油膜厚度随时间变化过程看(图8), 西北风阻止了油膜向上游漂移, 溢油事故对东风西沙、太仓水库取水口无影响, 对陈行水库和青草沙水库影响微小, 最大油膜厚度小于1μm。与主导风北风结果相比, 西北风与北风作用下陈行水库取水口油膜厚度差别显著。

对青草沙饮用水水源保护区, 溢油事故发生后3h 油膜到达, 23h 后离开, 持续影响20h, 期间最大厚度为5.12μm。对陈行饮用水水源保护区, 溢油事故发生后10h 油膜到达, 12.5h 后离开, 持续影响2.5h, 期间最大厚度为4.90μm。对长江太仓浏河饮用水水源保护区和东风西沙饮用水水源保护区无影响。

对于其他环境敏感保护目标, 西北风减轻了溢油事故点上游和北港、崇明东滩外侧敏感目标的影响, 但加重了南港和南槽的环境敏感保护目标的影响。

2.4 不利风向西南风3.2m·s–1 风速

溢油事故发生后, 油膜受涨潮流的作用向上游运动, 同时受西南风的作用向东北漂移, 大部分油膜绕过青草沙水库的西南角, 进入北港。在溢油事故后72h(图10), 油膜受向海的径流和西南风的作用, 大部分油膜在北港向下游漂移扩散, 位于崇明岛南岸的北港中段中北侧, 小部分油膜在北槽上端的中北侧。可见, 西南风使油膜漂移至北港, 并和径流一起作用向海输运。

图10 溢油事故后72h 油膜分布Fig. 10 Distribution of oil film thickness at 72 hours after the oil spill accident

从长江口水库取水口油膜厚度随时间变化过程看(图8), 溢油事故对东风西沙水库和太仓水库取水口无影响。对陈行水库取水口, 溢油事故发生后2h油膜到达, 8h 后离开, 持续影响6h, 期间最大厚度为6.3μm。对青草沙水库取水口, 溢油事故发生后9h 油膜到达, 54h 后离开, 持续影响45h, 期间最大厚度为8.5μm。

对青草沙饮用水水源保护区, 溢油事故发生后8h 油膜到达, 67h 后离开, 持续影响59h, 期间最大厚度为17.79μm。对陈行饮用水水源保护区, 溢油事故发生后1h 油膜到达, 41h 后离开, 持续影响40h,期间最大厚度为122.60μm。对长江太仓浏河饮用水水源保护区, 溢油事故发生后3.5h 油膜到达, 6.5h后离开, 持续影响3h, 期间最大厚度为15.02μm。溢油事故对东风西沙饮用水水源保护区无影响。

对其他环境敏感保护目标的影响, 西南风减轻了溢油事故点上游和下游南港、南槽环境敏感目标的影响, 但加重了北港和崇明东滩外侧敏感目标的影响。

3 结论

本文应用河口海岸三维水动力模式ECOM-si,耦合溢油漂移扩散模块, 模拟和分析长江口宝钢码头船舶发生溢油事故后油膜的漂移扩散, 以及对环境敏感保护目标的影响。基于长江口崇明东滩气象站2007—2021 共14 年的实测风速风向资料, 给出影响溢油漂移的主导风和不利风向。在冬季多年平均1月径流量11700m³·s–1和主导风、3 个不利风作用下,数值模拟和详细给出了上海港宝钢原料区域码头溢油事故发生后油膜的平面分布、油膜到达和离开4个水库取水口和饮用水水源保护区以及其他环境敏感保护区的时间、持续影响时间和最大油膜厚度。

长江口宝钢码头溢油事故发生后, 油膜的运动主要受水流的平流、风的漂移和随机扩散作用。在1 月主导风北风5.6m·s–1风速作用下, 油膜沿长江口南岸向下游漂移扩散, 在溢油事故后第72h 大部分位于南港中下段和南槽中上段的南侧, 小部分进入北槽南侧。在油膜的运动过程中, 绝大部分位于长江口南岸, 在主导风北风的作用下向南漂移起着重要的作用。在不利风向东南风4.0m·s–1风速作用下,油膜西北方向漂移, 聚集于南支北岸, 同时受径流作用向下游输运, 在溢油事故后72h 油膜位于南支的上下段的中北侧和北港上段的中北侧。在不利风向西北风4.8m·s–1风速作用下, 西北风减轻了溢油事故点上游和北港、崇明东滩外侧敏感目标受到的影响, 加重了对南港和南槽的环境敏感保护目标的影响。在不利风向西南风3.2m·s–1风速作用下, 在溢油事故后72h 大部分油膜在北港向下游漂移扩散,位于崇明岛南岸的北港中段中北侧, 小部分油膜在北槽上端的中北侧。西南风减轻了溢油事故点上游和下游南港、南槽环境敏感目标的影响, 加重了北港和崇明东滩外侧敏感目标的影响。油膜随涨落潮流作上下游来回振荡过程中, 径流使油膜向海输运,风使油膜朝风向方向漂移。不同风向作用下油膜的分布和对环境敏感保护目标影响显著不同, 表明风在油膜漂移扩散中起着十分重要的作用。