复杂河网水系油粒子模型开发及溢油污染模拟

王 鹏，申 霞，王船海，华祖林，马乙心，刘晓东

（1.河海大学浅水湖泊综合治理与资源开发教育部重点实验室，江苏南京 210098；2.河海大学环境学院，江苏南京 210098；3.南京水利科学研究院，江苏南京 210029；4.河海大学水文水资源学院，江苏南京 210098）

1 研究背景

随着全球经济的发展，世界各国对能源的需求与日俱增，海上石油运输量迅猛增长，水体溢油事故风险大为提升。为了预测溢油事故发生后的运动轨迹和扩散范围，科学评估其对生态环境的影响程度，提高溢油事件处置的应急响应能力，通常采用数值模拟方法对溢油自身扩展运动、在水流和风场作用下的漂移和扩散以及各种风化过程进行预测。溢油模型的研究始于1960年代，主要经历三个阶段：首先是油膜扩展模型研究阶段，以Fay理论为代表，考虑重力、惯性力、表面张力和黏性力作用，通过不同阶段油膜扩展的动力机制建立油膜直径经验表达式，取得了开创性成果［1］；第二阶段为基于对流扩散方程的数值离散模拟溢油运动，其可能会引入与物理扩散无关的数值扩散，甚至完全掩盖溢油的实际物理扩散过程，计算结果欠佳；第三阶段为Johansen［2］和Elliott［3］等人提出的油粒子模型，该方法将水体中的溢油视为众多离散油粒子的组合，粒子的平流过程具有拉格朗日性质，而剪切流和紊流引起的扩散过程属于随机运动，可用随机走动法来模拟，即将紊流视为一种随机流场，每个油粒子在紊流场中的运动类似于流体分子的布朗运动，运动的随机性使得整个粒子“云团”在水体中扩散，由于该方法对溢油扩散过程的模拟更加接近真实情景，目前应用较为广泛［4-6］。

在国外，油粒子模型最先开发并应用于海洋、海湾等水体的二维、三维溢油模拟，而在溢油概率较低的河网水系未得到应用。然而，我国河网地区的航运业普遍非常发达，随着船舶流量和货运量的快速上升，船舶溢油事故风险骤增［7］。长江、淮河、珠江等流域下游广泛分布着大量的复杂河网水系，这些地区河道沟渠纵横交错，山丘区树状河网和平原区网状河网并存，水系结构高度分异。此外，流域内水闸、泵站等水利工程众多，水系连通性受人工高度控制，加上海洋潮汐的驱动作用，导致河道水流流向往复不定，水动力特征十分复杂。由于河网水流运动过程通常采用一维水动力模型模拟，仅能输出沿河道方向的断面平均流速，无法据此判断油粒子运动至河网节点后流动方向，导致现有油粒子模型不能直接在河网地区应用。因此，在复杂河网水系采用油粒子模型开展溢油污染模拟，关键在于解决油粒子漂移至河网节点后的迁移和分配问题。

本文以一维河网水动力模型为基础，针对油品在该类水体的扩展、输移和风化特征，开发适用于河网水系溢油污染模拟的油粒子模型。一方面，使用随机走动法模拟溢油的惯性扩展、黏性扩展和表面张力扩展等机械扩展过程，解决了油粒子模型不能模拟溢油自身扩展阶段的难题［8］。另一方面，通过定义河网节点出流河道流向因子，判断油粒子运动至节点后的流动去向，实现了油粒子模型与一维河网水动力模型的结合。

2 材料与方法

2.1 模型构建本研究在“数字流域系统”［9-10］基础上进行油粒子模型开发工作，该系统的模型库包括流域水文模型、污染负荷模型、水动力模型和水质模型。由于河网地区各类水体的水深均较浅，系统采用一、二维耦合的水动力模型对水流运动过程进行模拟［11］，为油粒子模型提供流速信息。

溢油进入水体后，其行为和归宿包括两个过程：空间上的传输以及同时发生的油品风化。空间传输由自身扩展和输移两个阶段构成：扩展指油膜由于内部张力而导致的面积增大；输移过程指在环境动力要素的作用下溢油的迁移运动，包括水平方向的漂移、扩散过程。油品的风化引起溢油组成性质的改变，包含蒸发、乳化、分散、溶解、光氧化、生物降解、吸附沉降、水体混合扩散以及生物体内的代谢作用等。

2.1.1 对流过程 对流过程考虑水流流速和风漂流流速对油粒子漂移的综合影响。其对流过程的基本方程如下［12］：

式中：Si为油粒子位置；Uc为水流表面流速，m/s；Uw为风漂流流速，m/s；α为水流运动对油膜漂移的影响因子，通常取1.1～1.2。

（1）一维河网模型。一维河网模型水流流速采用一维河网水动力模型计算得到，风引起的油膜漂移速度大小表示为：

式中：β为风漂流系数，通常取3%～ 4%；w10为水面上10 m处风速，m/s；φ为风向与水流流向的夹角。

河网是由纵横交错的河道构成的水系，从物理结构上看，河道交汇处形成了若干节点。根据河道水流流向与河网节点的关系，可以将与节点相连的河道分为入流河道和出流河道。其中入流河道指的是水流流入节点的河道，水流流出节点的河道称为出流河道。由于河网地区地势低平，河道比降小，加上受潮汐和水利工程调控的影响，造成河道水流流向不定，因此与节点相连的河道是入流河道，还是出流河道，需要根据河网水动力模型的计算结果进行动态判断。

对于出、入流河道交汇形成的节点，在河网数值模拟中通常将其作为零维模型对象处理，由于零维模型无法模拟水动力过程，导致油粒子随水流运动至河网节点后，立即随出流河道的水流继续运动，因此节点大小和属性不会对油粒子漂移和归宿产生影响。同时，由于出流河道数量可能不止一条，为了预测油粒子通过节点后的运动轨迹，除了利用一维河网水动力模型预测流速之外，还需要判断油粒子运动至河网节点后的去向。

由于河网水系的河流宽度普遍较窄，溢油可以在较短的时间内在宽度方向上扩散均匀，加之绝大部分溢油漂浮于水面，则油粒子在节点处的流向概率与出流河道的宽度及其流速成正比。此外，考虑到这类地区的河道水深通常相差不大，因此，模型假设油粒子流入某条河道的概率与该条河道的出流流量大小成正比，即以与节点相连的各出流河道流量为权重判断油粒子流向。具体计算方法如下：

首先，将与某一节点相连的所有出流河道的流量从小到大排列，定义各条出流河道的流向因子为dfi，并按下式计算：

式中：dfi为第i条出流河道的流向因子；qi为第i条出流河道的流量；n为出流河道数量。

其次，定义各条出流河道的流向因子判断区间。其中，第1 条出流河道的流向因子判断区间为［0，df1］，第i条出流河道的流向因子判断区间为

最后，对于运动至河网节点处的每一个油粒子，生成［0，1］的均匀随机数Rd，如果Rd位于某条出流河道的流向因子判断区间，则油粒子将随水流漂移至该条出流河道。

（2）平面二维模型。对于平面二维模型，在风力作用下，油粒子漂移方向与风向成0°～40°夹角，此时风引起的油膜漂移速度可表示为［12］：

式中：D为考虑风向偏角的转换矩阵；θ为地转科氏力引起的风向偏角。

2.1.2 机械扩展和紊动扩散 溢油机械扩展及剪切流和紊流引起的粒子紊动扩散过程，均采用随机走动方法模拟。通过在流场中追踪各质点的运动轨迹，得到每一时刻各个油粒子所处的空间位置，统计各时刻油粒子的位置可得到各时刻溢油的空间分布。

根据分子紊动扩散理论，在一维空间情况下，随机走动方差与扩散系数之间的关系可以表示为［12］：

式中：x为随机走动距离，m；K为扩散系数，m2/s；Rn为均值为0，标准差为1.0的正态分布随机数；Δt为时间步长，s。

因此，对于一维河网模型，油粒子每一个时间步长的随机走动速度可采用下式计算：

式中：Vt为紊流扩散的随机走动速度，m/s；De为油膜机械扩展系数，m2/s；DT为紊动扩散系数，m2/s。

同理，对于平面二维模型，每一个时间步长的随机走动速度可表示为［8］：

式中δ为［0，π］之间的均匀分布随机角。

机械扩展系数和紊动扩散系数的计算采用Sayed［8］提出的方法。对于机械扩展系数，将溢油的机械扩展过程分为惯性、黏性、表面张力三个阶段，首先计算油粒子实际厚度及油粒子厚度的临界值，再根据两者之间的大小关系确定扩展阶段，最后计算不同扩展阶段对应的机械扩展系数De。紊动扩散系数考虑摩阻速度和水深的影响。

2.1.3 风化过程 溢油事故发生后，除了伴随着扩展、对流、扩散等动力过程外，油品还经历如蒸发、乳化、分散、溶解、光氧化及生物降解等风化过程，使油膜质量、油膜物理化学性质等发生一系列变化，这些变化主要和油品自身性质以及海况条件如风、波浪、水流、气温以及生物活动等有关。本研究建立的溢油预测模型主要用于突发溢油污染应急预测，分析溢油发生后短期内的油膜迁移扩散规律，为短期模型，而光氧化、生物降解等过程相比其它风化过程是一个非常缓慢的过程。此外，分散会影响溢油的垂向分布特征，但是不会对溢油在水体中的整体质量变化产生影响。因此，本研究重点对溢油初期风化作用较突出的蒸发和乳化两个过程进行模拟。

（1）蒸发过程。蒸发是油品中的石油烃轻组分从液态变为气态的过程，是溢油初期阶段油品与大气物质交换的重要过程，是使油品残留量大幅减少的重要途径。蒸发过程与油品组分、油膜厚度以及环境状况等因素相关。对于原油泄露，模型采用目前广泛使用的Stiver等［13］提出的计算模式。对于柴油或汽油等燃料油泄露导致的蒸发损失，采用Fingas 模型［14］计算，该模型将石油及其产品的蒸发表示为时间、温度和180℃时质量蒸馏比的函数，认为蒸发过程满足平方根函数或对数函数形式。

（2）乳化过程。溢油乳化是另一项重要的风化过程，指油类吸收水而形成油包水乳化液的过程。乳化使油滴体积增加2～5倍，主要受风速、波浪、油膜厚度、环境温度、油风化程度等因素的影响。乳化作用通常在油膜拓展较大、厚度较薄时发生。风浪能量打碎油膜，水滴分散到油中，形成油包水的乳化液，呈黑褐色黏性泡沫乳油状漂浮于水面。乳化物含水率采用Mackay 等提出的乳化作用方程［15］，该方程也被其它大多数溢油模型所采用。

2.1.4 河岸吸附 河网地区的河道宽度普遍较小，因此溢油在自身扩展阶段就可能接触到水陆边界。当溢油接触到河岸时，有可能被河岸本身或滩涂、坑洼和植物吸附。研究表明［16］，河岸对溢油的吸附作用有三种情况：完全吸收、完全反射和部分吸收。然而，由于河岸在材质、形状、水力条件、植被等方面存在较大差异，对其吸收能力只能进行粗略分析。为了简化计算，本研究对于河岸吸附作用模拟只考虑完全吸收和完全反射两种情况，通过假设油粒子吸附概率，判断其吸附状态［17］。

2.2 输入与输出根据以上对溢油模型对流、扩散和风化过程的分析，模型输入数据包括以下几个方面：（1）溢油事件时空特征：溢油发生时间、溢油点坐标等；（2）溢油质量及油粒子数量；（3）油品初始理化性质，具体包括油品的初始温度、密度、蒸发参数；（4）气象数据，包括风速、风向等。

为了反映溢油对水环境的影响范围和程度，模型输出部分包括各计算单元的油粒子数量、质量、体积。其中油粒子体积考虑蒸发和乳化过程的影响，其剩余体积可表示为：

式中：V0为油粒子的初始体积，m3；Vi为蒸发和乳化后油粒子的剩余体积，m3；Fvi为第i个油粒子的蒸发率；Ywi为第i个油粒子乳化物的含水率。

初始密度由用户根据油品种类设定，考虑温度、蒸发和乳化作用的影响。

2.3 研究区选择太湖流域位于长江下游地区，流域面积36 900 km2，流域内河网密布，湖泊众多，水面率达17%，平原地区河道密度达3.2 km/km2，属于典型的平原河网水系。为了提升流域在水资源保护、水环境治理等方面的管理水平和决策能力，经过二十多年的研制和不断完善，以基于双对象共享结构的数字流域系统为支撑［9］，开发了“太湖流域水量水质决策支持系统”。系统模型库由水文模型、污染负荷模型、水动力模型、水质模型组成，可对流域产汇流过程、污染物在陆域和水体的迁移和转化过程进行模拟。系统分别采用零维、一维和二维模型对全流域的97 个小型湖泊、952 条河道及太湖湖体进行概化，并包含188 座闸泵工程。其中一维模型共概化河道断面4306 个，河道总长7958 km，二维模型生成网格2339个。太湖流域概化河网及太湖计算网格如图1所示。

图1 太湖流域概化河网及太湖计算网格

2.4 模型率定和验证以“太湖流域水资源综合规划数模研制项目”为依托［18］，采用沿长江13座水利工程、沿杭州湾3 座水利工程、望亭水利枢纽及太浦闸全年引排水量资料作为“太湖流域河网水动力模型”的边界条件，对2000年太湖及12个水位代表站、4个流量代表站的全年水位和流量过程进行了率定，同时利用环太湖流量巡测资料和流域供、用、耗、排水量调查资料对出入太湖水量和流域水量进行了平衡分析；此外，还利用1998、1999年水量实测和调查资料，对模型计算结果进行了验证。率定和验证结果均表明，太湖及大部分代表站的特征水位、全年水位和流量过程线预测值与实测值相差较小，各控制线水量平衡分析成果能够反映太湖流域的水流实际运动状况［19］。

由于溢油事故发生具有突发性和不可预见性，尤其对于水流流向多变的河网水系，很难在这类地区获得水位、流量、溢油浓度的高频率同步观测资料，因此，基于观测数据对溢油模型预测结果进行率定验证的难度较高。本文通过设定溢油情景的方法开展数值试验，重点对油粒子在河网水系迁移和分配过程的合理性进行分析。

2.5 溢油情景设置将太湖流域河网水系作为溢油污染数值试验场景，耦合油粒子模型与“太湖流域河网水动力模型”，通过设置假想溢油事件，开展全流域溢油漂移和归宿模拟。采用流向因子判断法模拟油粒子在河网水系的迁移和分配过程，分析不同典型时刻油粒子空间分布特征，探讨溢油预测结果的合理性及预测方法的科学性，评估自然风化过程和人为打捞措施对油粒子漂移和扩散的影响。

太湖流域自2002年开始实施“引江济太”工程，以提升流域水资源与水环境承载能力。太浦河作为“引江济太”工程的主要输水通道、上海和嘉兴等地的重要水源地及芜申线（三级航道）重要组成部分，同时具有泄洪、排涝、供水、航运等多种功能，一旦发生溢油事故，将对水源地水质及周边地区供水安全产生不利影响。因此，假定在太浦河太浦闸下游约5 km 处发生货船碰撞溢油事故，泄漏油品为柴油，泄漏量为2 t，泄漏持续时间10 min，油粒子个数5000 个。模型计算时间为2000年1月1日8∶00至2000年1月20日8∶00，假设溢油事件发生在1月1日10∶00，时间步长5 min。

3 结果与讨论

按上述设定的溢油情景开展河网突发溢油污染模拟。首先，以某一节点为例，分析油粒子运动至河网节点后的流动去向，验证出流河道油粒子数量与出流流量大小的关系。图2为2000年1月5日油粒子漂移至某河网节点处，并随水流继续向下游漂移的过程。如图2所示，油粒子在1月5日00∶00 运动至太浦河与4条支流交汇形成的节点，随后一部分油粒子沿太浦河继续向下游运动，另一部分油粒子沿支流2运动。此时，支流1、支流3和支流4均为入流河道，太浦河及支流2为出流河道，两条河道的油粒子数量分别为1683 个和114 个，流量分别为297.3 m3/s 和20.3 m3/s。由此可见，出流河道油粒子数量与其流量大小成正比，说明流向因子判断法能够以出流河道流量为权重，正确模拟油粒子在河网水系的迁移和分配过程。

图2 油粒子在河网节点的迁移和分配

对于复杂河网水系，由于其河流宽度普遍比较狭窄，加之油膜在溢油初始阶段的机械扩展速度较快，可以在溢油后的较短时间内沿河宽方向扩散均匀，因此，油粒子沿出流河道的漂移数量与出流河道宽度成正比，即出流河道宽度越大，向该河道漂移的油粒子数量越多。此外，出流河道流速大小也是影响油粒子漂移方向的关键因素，即出流河道流速越大，油粒子向该河道漂移的概率越高。同时，河网地区各河道水深通常相差不大。因此，综合考虑油粒子流向概率与河宽、流速及水深的关系，以出流河道流量为权重判断油粒子流向的方法是合理的。

目前，油粒子模型在河流溢油污染模拟中开展了一些应用研究［20］，取得了较好的预测效果，但是在油粒子漂移方向判断上，其与现有应用于海洋、海岸等开敞水域的油粒子模型并无本质区别，也是基于二维或三维水动力模型的流场预测结果，按照对流方程计算其漂移轨迹，仍然无法应用于众多河流交汇形成的复杂河网水系，使油粒子模型的应用场景受到一定限制。本研究通过定义出流河道流向因子及其判断区间，在流向概率与流量大小间建立联系，从实际预测效果来看，较好地解决了一维河网水动力模型预测结果无法决定油粒子漂移方向的问题。

为了进一步分析自然风化过程对溢油漂移和归宿预测的影响，将不同典型时刻油粒子的空间分布绘于图3。如图所示，在溢油发生的初始阶段，由于各支流流向大部分为流入太浦河，因此，从1月1日10∶00 到1月5日00∶00，油粒子总体随水流沿太浦河向下游移动，部分溢油短时间进入支流，最终又汇入太浦河，因此形成了两条溢油带，总长度约为900 m（图2（a）），在蒸发和河岸吸附等风化因素的作用下，油粒子总质量减小为854 g，相对于初始质量减少了57%。随后，油粒子开始向周边河网水系扩散，进入太浦河支流，影响范围逐渐扩大（图2（b））。1月8日11∶00 溢油进入拦路港，此时距事故发生169 h，油粒子总质量为528 g。1月9日13∶00 油粒子进入黄浦江，总质量进一步减至398 g。最后，受上游径流和下游潮汐的双重影响，溢油沿黄浦江往复震荡并逐渐向下游入江口漂移。

图3 典型时刻油粒子空间分布

除了自然风化过程，人为打捞措施也是影响油粒子漂移和扩散的重要因素，模型通过设置打捞范围和打捞率，模拟围油栏、吸油毡、消油剂等溢油回收和消除措施。假设在溢油事故发生24 h 后采取打捞措施，打捞率设置为50%。图4为溢油打捞情景下典型时刻油粒子空间分布图，由图可见，油粒子进入拦路港和黄浦江的时间与未采取打捞措施的方案几乎没有差别，进一步说明对流过程是控制溢油在河网水系输移的主导因素。从油粒子空间分布来看，打捞措施对减少局部水体的油粒子数量有较大影响，并且该影响随着溢油输移范围的扩大而增强。

图4 溢油打捞情景下典型时刻油粒子空间分布

4 结论

近年来，油粒子模型在溢油输运和归宿模拟方面获得了巨大成功，已被广泛应用于海洋、海湾等开敞水域的二维、三维溢油数值模拟。然而，由于河网地区的水流运动模拟通常采用一维水动力模型，导致现有的油粒子模型不能与其耦合，无法在河网水系开展溢油污染模拟。本研究通过定义河道流向因子及其判断区间，确定油粒子运动至河网节点处的流动去向，开发了一种能够与一维河网水动力模型耦合的油粒子模型。通过在太湖流域开展溢油数值试验，分析了自然风化过程和人为打捞措施对油粒子漂移和扩散的影响。结果表明，流向因子判断法能够正确模拟油粒子在河网水系的迁移和分配过程，基于该方法开发的油粒子模型适用于水流条件多变的复杂河网水系，可对溢油在该类地区的扩散范围及其水质影响进行科学评估。目前，该研究成果已成功集成于《太湖流域水环境综合治理总体方案》重点项目“太湖流域水资源监控与保护预警系统”，具备了在该流域开展溢油污染模拟的应用条件，为提升复杂河网水系突发水污染预警与决策的科学水平提供了有力的技术支撑。