长江宜昌段溢油事故影响数学模型分析研究

胡　旭，辛小康，尹　炜，陈和春，王继保，李光浩

（1.三峡大学　水利与环境学院，湖北　宜昌　443002；2.长江水资源保护科学研究所，武汉　430051）

长江是我国主要的河流运输通道，客货运输量大。交通运输部发布《2016年交通运输行业发展统计公报》显示[1]，内河运输完成货物量达35.72亿吨，占水路总运输货物量的55.97%；内河港口货物吞吐量完成47.76亿吨，增长3.1%。随着国家依托黄金水道推动长江经济带发展战略的实施，内陆自贸实验区相继开放，货运量在今后仍将继续增加[2]。在同条件及预警技术管理水平不变前提下，这无疑增大了溢油事故风险发生的概率。众所周知，溢油事故产生的石油类污染物会给整个河流水体造成巨大水环境污染，严重破坏水环境平衡，危害水生生物。相关研究表明，中华鲟在含油浓度为10 mg/L时，死亡时间仅为2天[3]。

为了更好地应对内河溢油事故风险，有必要对内河溢油事故进行模拟研究。目前，工程界和学术界对溢油事故的模拟研究多采用数值模拟方法[4-5]。现有研究对象多数为入海口及海洋地区。对于环境差异较大、水文条件特殊、河道情况复杂多变的内河，溢油事故研究较少[6-9]。本文采用MIKE21软件，建立平面二维溢油模型，模拟葛洲坝至下游枝江江段丰水期和枯水期突发溢油事故，并分析不同溢油量条件下油膜的扩展、输移范围及其对下游造成的影响。

1　模型计算原理

1.1　水动力模块

MIKE21是由丹麦水力研究院研究开发的应用于近海海域、海湾、港区及河道等水环境模拟的系统。水动力模块（HD）是其核心，可以模拟由于各种力的作用而产生的水位及水流变化。其控制方程是基于不可压缩流体和Reynolds值均布的Navier-Stokes方程，并遵从Boussinesq假设和静水压力假设。二维非恒定浅水方程组如下：

式中，t为时间，η为水面高程，h为总水深，g为重力加速度，ρ为水的密度，ρ0为（淡）水的参考密度，f=2Ωsinϕ为科氏力系数（Ω为旋转角速率，ϕ为地理纬度），Pa为大气压强，Sij为辐射应力张量，S和(us,vs)分别为点源的排放量和速度，uˉ和vˉ分别为沿x，y方向上的流速在深度上的平均值；(τsx,τsy)和(τbx,τby)为水面风应力张量和河床床面应力张量，Tij为侧向应力。

该模块利用隐式交替方向ADI（Alternating di⁃rection implicit）技术对上述连续方程和动量方程进行离散，所得矩阵用追赶法求解，各项微分和主要系数均采用中心差分格式，Taylor级数展开的截断误差可达二至三阶精度。

1.2　溢油模块

MIKE21的Oil Spill模块（OS模块）是在HD模块基础之上建立的溢油模型，用于模拟溢油扩散运动规律和归宿。溢油模块采用“油粒子”对溢油量进行概化，采用拉格朗日方法描述油粒子的输移轨迹，包括扩展、漂移、扩散等过程[9]，其输移方程见公式（4）-（7）。同时，在输移过程中油粒子也会发生如蒸发、乳化、溶解等风化过程，该软件通过计算油粒子质量损失来体现。

（1）扩展：油膜扩展运动采用修正的Fay重力-粘力公式计算。公式如下：

式中：Ao为油膜面积；t为时间；Kα为扩散系数；Vo为油膜体积；Ro为油膜半径，h0为初始油膜厚度，取10 cm。

（2）漂移：影响油粒子漂移速度的主要因素是水流和风作用力。漂移速度为:

式中：UW为水面上10 m处的风速，Us为表面流速，cw为漂移系数，一般取值0.02～0.03。

（3）扩散：由于单个粒子不能被分成几片，因此扩散的过程被解释为在随机方向上的运动。对于二维的情况，可以将随机走动的距离形式表示为一个时间步长α方向上的扩散距离。计算公式如下：

式中：Sα为在α方向上的一个时间步长内可能扩散走动的距离，Dα为α方向上的扩散系数，R为-1到1的随机数。Δtp为计算时间步长。

2　研究区域概况

本文模拟范围为长江葛洲坝至下游枝江江段。上游起始断面为葛洲坝大坝，下游终止断面为枝江昌门溪，范围全长80 km。根据《全国重要江河湖泊水功能区划（2011-2030年）》，该江段包含长江宜昌饮用水源、工业用水区和长江宜昌中华鲟保护区两个水功能区。根据水功能区管理要求，长江宜昌饮用水源、工业用水区水质目标为Ⅱ类，水质目标要求较高。模拟范围如图1。

图1　模拟范围内各水功能区、监测断面、概化溢油事故点相对位置示意图Fig.1　Relative position of each water function area,monitoring section,and generalized point of oil spill accident

3　模型参数设定及计算方案

3.1　网格划分

本文采用2015年实测水下地形数据。葛洲坝至枝江区间的长江干流，河道形状相对规则，变化梯度不大。为了保证更高计算精度和收敛性，模型网格划分时采用四边形网格，网格尺寸为200 m×40 m；支流清江入汇口段，地形复杂，区域较小，采用三角形网格划分，三角形网格控制面积小于1 000 m2。计算区域共剖分计算网格单元16 093个，并根据该实测水下地形进行插值，如图2和图3所示。

图2　葛洲坝-枝江区间水质数学模型计算网格Fig.2　Computing grid for water quality mathematical model between Gezhouba and Zhijiang section

图3　葛洲坝-枝江区间计算网格地形插值Fig.3　Terrain interpolation for computing grid between Gezhouba and Zhijiang section

3.2　水动力模块边界条件及参数设定

葛洲坝是三峡的反调节水库，为日调节水库，设计水文条件可参照三峡出库流量。三峡水库2003年才开始初期蓄水，2010年才开始175 m试验性蓄水，水文情势尚不稳定。于是本文便采用三峡水库近10年最枯、最丰月的月均流量作为设计水文条件。本模型枯水期上游边界设计流量为三峡2月份平均出库流量5 280 m3/s，参照枝江水位站水位—流量关系曲线，得出模型下游边界设计水位为33.42 m。丰水期上游边界设计流量为三峡8月份平均出库流量21 357 m3/s，对应下游边界设计水位为36.15 m。参照多年平均水温，枯水期设定水温10℃，汛期设定水温20℃。

设定湿水深度为0.100 m，淹没深度为0.050 m，干水深度为0.005 m，根据辛小康[10]等对长江宜昌江段水动力模型的率定结果，曼宁系数为0.031，涡粘性系数的Smagorinsky系数为0.28，扩散系数D设定为1 m2/s。

3.3　溢油模块参数设定

本文首先计算各个油粒子的位置变化和组成变化，然后统计各个网格上的油粒子个数和组分含量，模拟出油膜的浓度时空分布和组分变化。MIKE21中溢油模型将单个油粒子质量定义为挥发性轻组分的质量、不挥发性重组分的质量、蜡状物的质量以及沥青的质量。本文只考虑轻组分和重组分，两者比例为3∶2，即1个1 000 kg油粒子，包含轻组分600 kg，重组分400 kg。由于石油类污染物易发生蒸发、溶解、乳化、沉淀、生物降解、光氧化等特点，本文参考姜卫星[11]2007年对狭长型感潮河流溢油模型的研究，采用纵向扩散系数为0.25 m2/s、横向为0.10 m2/s。风速和风向，可根据气象资料获得。

枯水期葛洲坝下游江段石油类背景浓度采用宜昌黄陵庙断面2010-2014年2月份实测平均值0.03 mg/L。丰水期采用该断面2010-2014年8月份实测平均值0.01 mg/L。

3.4　计算方案

本次模拟设计溢油点位置为该区域航运交通要塞三江引航道和大江引航道的汇合点，坐标经纬度为（E111°16′4.8″，N30°42′10.8″），垂向水深为零。参考我国《船舶油污染事故等级》（JT/458-2001）[12]，考虑最不利情形，即油船造成的水域油污染事故。根据油船大事故和重大事故的划分依据，拟取设计总溢油量分别为5 t和10 t，并分别模拟枯水期和丰水期两种不同水文条件下的溢油事故情形，共4种模拟工况，提取每种工况1 h，3 h，10 h后的溢油模拟计算结果进行对比分析。

4　计算结果与分析

4.1　枯水期

枯水期瞬时溢油为10 t时，事故发生后1 h，3 h和10 h以后油膜位置和面积如图4所示。可以看出，当溢油10 t后1 h，油膜尚未完全展开，油膜厚度≥0.3µm的范围为600 m×300 m；油膜厚度大部分为100µm～0.1 mm的范围为400 m×200 m。溢油后3 h，油膜厚度≥0.3µm的范围为3600 m×300 m；油膜厚度为100µm～0.1 mm的范围为3 200 m×280 m。溢油10 h以后，油类物质经过扩散，污染带范围进一步增加，但由于乳化、蒸发、风化、光解等综合作用，油膜厚度有所减小，油膜厚度≥0.3µm的范围可达到14 000 m×400 m；油膜厚度为100µm～0.1 mm的范围为600 m×200 m。

图4　枯水期溢油10 t后1 h，3 h，10 h油膜厚度分布图Fig.4　Oil film thickness distribution at 1 h,3 h and 10 h after 10 t spillage during low water period

枯水期瞬时溢油5 t后1 h，3 h和10 h油膜演进情况见图5。可以看出，当溢油后1 h，油膜厚度≥0.3µm的范围为400 m×240 m；油膜厚度大部分为100µm～0.1 mm，其范围为300 m×200 m。溢油后3 h，油膜厚度≥0.3µm的范围为2 200 m×300 m；油膜厚度为100µm～0.1 mm的范围为2 000 m×260m。溢油10 h以后，油类物质经过扩散，并有少量附着在河岸，污染带范围进一步增加，油膜厚度有所减小，油膜厚度≥0.3µm的范围可达到10 000 m×400 m，而油膜厚度在100µm～0.1 mm几乎消失。

图5　枯水期溢油5 t后1 h，3 h，10 h油膜厚度分布图Fig.5　Oil film thickness distribution at 1 h,3 h and 10 h after 5 t spillage during low water period

4.2　丰水期

丰水期瞬时溢油10 t时，事故发生后1 h，3 h和10 h以后油膜位置和面积如图6所示。可以看出，当溢油10 t后1 h，油膜厚度为100µm～1 mm，范围为800 m×100 m。溢油后3 h，油膜面积达到4 400 m×200 m，油膜厚度超过100µm。溢油10 h以后，油类物质经过扩散，污染带范围进一步增加，油膜范围达到10 000 m×400 m，同样由于乳化、蒸发、风化、光解等综合作用，油膜厚度减小，基本小于100µm。

图6　丰水期溢油10 t后1 h，3 h，10 h油膜厚度分布图Fig.6　Oil film thickness distribution at 1 h,3 h and 10 h after 10 t spillage during high water period

丰水期瞬时溢油5 t时，事故发生后1 h，3 h和10 h以后油膜位置和面积如图7所示。可以看出，当溢油5 t后1 h，油膜尚未完全展开，油膜厚度为100µm～1 mm，范围800 m×100 m。溢油后3 h，油膜面积达2 600 m×200 m，油膜厚度超过100µm。溢油10 h以后，油类物质经过扩散，飘移，污染带范围进一步增加，油膜范围达到8 000 m×300 m，而油膜厚度减小，基本小于100µm。

图7　丰水期溢油5 t后1 h，3 h，10 h油膜厚度分布图Fig.7　Oil film thickness distribution at1 h,3 h and 10 h after 5 t spillage during high water period

4.3　结果分析

根据模拟结果可知，对于所有工况，漏油5 t以上均会产生明显的油膜厚度≥100µm的区域，且持续时间都超过3 h。溢油1 h后油膜处于或刚完成溢油扩散的初始阶段，油膜展开范围不大，厚度均大于100µm。溢油3 h后，油膜范围继续扩大，油膜影响范围均大于2 000 m×260 m，最大可达3 600 m×300 m。溢油10 h以后，油类物质经过扩散，污染带范围进一步增加，无论是枯水期还是丰水期，两种溢油量油膜影响范围均大于8 000 m×300 m，最大可达10 000 m×400 m，但由于物理、化学等综合作用，油膜厚度减小较明显，普遍小于100µm。在同等水力条件下，溢油量越多油膜相对更厚，但油膜扩散范围的差异较小，这可能与河道较狭窄有关。

对比枯水期和丰水期可知，枯水期由于水量偏小，流速小，油污范围较小，持续的时间较长，而丰水期水量大，流速大，油膜扫过的范围相对更大更远，油膜厚度减小速率更大，持续的时间也相对较短。枯水期油膜扫过的范围可达到猇亭一带，丰水期可达到枝城一带。无论是枯水期还是丰水期，当瞬时溢油量为5 t及以上时，均会对中华鲟自然保护区产生明显影响。

5　结论

我国内河流域环境普遍复杂多变，与海洋有较大差别，对于内河溢油事故模拟，需要考虑的因素众多。本文根据实测地形资料借用MIKE21软件对宜昌江段进行溢油事故模拟预测，考虑了不同水文条件，不同溢油量的情形，而模拟结果也表明这对油膜厚度、飘移范围、持续时间的差别有较为显著的影响，这对内河溢油应急对策具有指导意义，也为内河溢油事故对水环境影响评估提供了很好的参考价值。

[1]中华人民共和国交通运输部.交通运输行业发展统计公报[R].北京：综合规划司，2016.

[2]袁群.浅谈长江航运油污染现状、原因及其对策[J].水运管理，2004（5）：22-23.

[3]倪朝辉，翟良安.石油对鱼类等水生生物的毒性[J].淡水渔业，1997，27（6）：38-40.

[4]陈家兴，杜娟，付金锋.基于GIS的三峡库区溢油模型研究[J].中国水运，2010，10（9）：96-98.

[5]高龙骧.内河船舶溢油扩散的建模与仿真[D].上海：复旦大学，2013.

[6]赵琰鑫，王永桂，张万顺，等.河道溢油污染事故二维数值模型研究[J].人民长江，2012，43（15）：81-84.

[7]GARCIA M R，TOVAR F H.Computer modeling of oil spill trajectories with a high accuracy method[J].Spill Science&Technology Bulletin，1999，5（5/6）：323-330.

[8]巫丽俊.感潮江段事故溢油二维数值模拟[D].南京：河海大学，2006.

[9]张帆，黄立文，邓建，等.重庆主城区江段溢油模型及数值试验研究[J].武汉理工大学学报（交通科学与工程版），2011，35（1）：87-90.

[10]辛小康，叶闽，尹炜.长江宜昌江段水污染事故的水库调度措施研究[J].水电能源科学，2011，29（06）：46-48，95.

[11]姜卫星.黄浦江溢油事故的数值模拟研究[D].上海：同济大学，2007.

[12]中华人民共和国交通部.船舶油污染事故等级：JT/458-2001[S].北京：中国标准出版社，2001.