崖城油气田附近海域溢油漂移轨迹数值模拟预测

王璟，娄安刚，曹振东

（中国海洋大学环境工程与科学学院，山东青岛266100）

崖城油气田附近海域溢油漂移轨迹数值模拟预测

王璟，娄安刚，曹振东

（中国海洋大学环境工程与科学学院，山东青岛266100）

通过采用不规则的三角网格和有限体积法的FVCOM模式，建立三维潮流数值模型。利用大海域计算得到的调和常数值作为开边界的输入值，模拟出崖城附近海域的潮流和潮位变化情况。在潮流、潮位验证正确的前提下，利用欧拉—拉格朗日追踪方法，建立了溢油轨迹预测模型，进行崖城油气田附近海域溢油中心轨迹的预测，同时预测了溢油漂移的平均速率和油膜抵达敏感区的时间，为油气田实施应急措施提供技术支持。

FVCOM模型；崖城；溢油；数值模拟

1 引言

崖城位于海南岛最南端，距离三亚市西40 km处，其附近海域拥有大量的油气田，在油气田开发和运营的过程中可能会对周围的海域和一些敏感资源造成一定的污染[1]，因此崖城附近海域水动力场的改变直接影响到这些工程造成的污染物输送。本文以崖城13-4平台附近发生突发性溢油为例，通过建立三维潮流模型模拟出崖城附近海域的潮流潮位情况，并以此为基础分别在平均风和最强风下预测出溢油的漂移轨迹、平均速率、抵岸时间等。

在之前对南海的研究中，多以北部湾[2-3]为主要研究对象，而对海南岛南部海域的研究较少。本文重点针对海南岛南部海域，对其潮流潮位随时间变化情况加以研究。

2 潮流数值模型

2.1 流体动力学基本方程组

运用FVCOM三维模型进行数值模拟，采用sigma坐标，三维非线性潮波运动方程为：

式中：g为重力加速度；f=2Ωsinφ为柯氏参数；Ω为地转角速度；φ为地理纬度；ζ为自静止水面算起的水位高度；h为自静止水面算起的水深；H=h+ζ；u、v、w为流速，p为压强，ρ为海水密度，τ为切应力。

2.2 边界条件

海面边界条件：

固体侧边界条件：

在开边界上水位边界条件：

这里H和g分别是调和常数的振幅和迟角，下标C代表着M2、S2、K1、O1四个分潮中的一个。ω代表着分潮频率，f为交点因子，u为交点订正角，V0是天文潮的初位相。

2.3 计算海域和网格设置

由于计算域开边界上调和常数实际中难以得到，本文首先建立了大海域海流模型（见图1），由大海域计算得出的调和常数值[4]，作为计算海域开边界的输入值。

图1 大海域网格及地形图

本文所建立的海域数学计算域范围及网格见图2，从图中可以看出，三角形网格较好的反映出计算域内的海岸线和地形特征。其中坐标范围为16.5°—19°N，108°—110.5°E，水深及岸线根据海图确定。崖城13-4气田周围海域网格加密，最小网格边长800 m。垂向δ分20层。

图2 计算海域网格及地形图

水界强迫水位：利用大海域的计算结果作为工程海域水界的输入资料。水界可提供4个分潮(M2、S2、K1、O1)的调和常数。由于是不正规全日潮海域[3]，水界可输入下列调和方程：

式中，右方第1项表示平均日分潮，第2项表示平均半日分潮。

3 模拟结果验证分析

3.1 潮流与潮位的验证

在计算海域中，选取A、B、C三点（坐标位置见表1）分别作为潮流验证点和潮位验证点。在这些点分别将数值计算的结果与实测的资料进行了验证，其验证结果见图3。

表1 验证点坐标位置

图3 流向验证曲线

从图3c—3h可以看出，流速的大小以及方向，转流发生时刻的计算值与实测值基本一致；图3a和3b显示A、B两站位的潮位振幅和位相计算值亦与实测值基本一致。

潮位和潮流的验证结果表明建立的潮流模型是可行的，适合本海区。

3.2 评价海域潮流计算结果

图4、5给出涨、落潮中间时的表层计算潮流场的潮流场图。由图中可以看出崖城13-4气田附近潮流分布表现为往复流：涨潮为西北向流，落潮为东南向流，落潮潮流流速大于涨潮潮流流速。通过计算得知，崖城13-4气田附近区域涨潮最大流速为48 cm/s，平均流速为31cm/s，落潮最大流速为60 cm/s，平均流速为38 cm/s。海南岛东南部海域流速远低于其他海域流速，这是由于东南部有一由100 m变深至1000 m的海底陡坡，与文献[5]的结论一致。

图4 涨潮中间时计算潮流场

图5 落潮中间时计算潮流场

图6和图7为崖城13-4气田附近区域涨潮时和落潮时垂向最大流速分布图，由图可以看出，最大流速随着深度的增加而减小，与文献[6]的结论一致。

图6 涨潮中间时最大流速垂向分布

图7 落潮中间时最大流速垂向分布

4 溢油预测模型

4.1 油膜漂移轨迹

在环境动力模型提供的环境动力参数的基础上，采用欧拉-拉格朗日追踪方法，进行油膜中心轨迹的预测[7]。油膜中心漂移速度，取决于海面风速与表层流，是空间和时间的函数，其值用油膜中心点所在网格点上的速度内插而得。空间每个网格节点上的x、y方向上的速度在某时刻为：

式中，Vrx、Vry为网格点上表层潮流流速加上风海流的x、y方向分量，表层潮流流速由环境动力学模型求出，风海流的模拟根据埃克曼公式计算：

式中，w为风速，Ф为纬度。Vwind为网格点上的风速,α为风因子，θ0为风向，θ为油粒子受风影响的漂移偏角。θ的取值与风速的大小有关，公式为：

油粒子漂移轨迹计算公式为：

式中：S0为初始时刻，S为油膜中心点所在位置，Vl(x(t),y(t),t)为拉格朗日追踪速度，

由于空间和时间不同，流况不同，有时风速、风向也不同，所以在不同地点、不同时刻发生溢油后所追踪到的油膜中心运移轨迹就不同。

4.2 油膜扩展输移预测

剪流和湍流引起的扩散过程属于随机运动，可用随机走动法实现模拟。由于每个粒子的随机运动而导致整个粒子云团在水体中的扩散过程。对于水体表面随机扩散过程可用下式描述[8]：

式中：ra为α=(x,y)方向上的湍动扩散距离；R为[-1，1]间均匀分布随机数。kα为α方向上的湍流扩散系数，Δt为时间步长。

第i个粒子在Δt时段内的位移可表示为：

式中：rx，ry为在x，y方向上的随机移动距离；ui，vi为拉格朗日速度在x，y方向上的分量。

4.3 油的化学过程

溢油在输移和扩展的过程中，也同时发生各种生物和化学变化过程，使得海上溢油的量不断减少。

4.3.1 溢油的挥发

溢油挥发过程受油性质、油厚度、风及油组分控制。采用Stiver和Mackay提出的一个暴露模式来计算油的挥发：

式中，B为系数，常取10.3；TG为挥发曲线梯度；T为油的表面温度，通常与大气温度相近；T0为初始时油挥发温度；Pa为大气压；V为油分子体积；R为大气常数；θ为挥发系数，常取2.5×10-3，Uw为风速；T0、TG的数值常参考如下常数：

由YC13-1气田混合凝析油的试验得知API度(15.56℃)为43.70。

4.3.2 溢油的乳化

乳化过程受风速、波浪、油的厚度、环境温度、油风化程度等因素的影响，一般用含水率来表示乳化程度。

式中：

YWi为第i个油粒子含水率，Uw为风速，WAx为油的含蜡量，Asph为油的沥青质量含量，η0为油的无水动力粘性系数，YWsat为稳定含水量，K1、K2为常数，分别为5.0×10-7和1.2×10-5。ηi为乳化后油的运动粘性系数。

式中，ηoil为乳化前油的运动粘性系数。

4.4 溢油量及溢出方式

4.4.1 溢油点

崖城13-4气田主体区，坐标17°37'35.697"N，109°09'44.123"E。

4.4.2 溢油方式及溢油量

点源一次性连续排放2 h，溢油量为50 m3，平均溢油速率为25 m3/h。

4.5 风场

根据提供的风速风向分布图以及风向玫瑰图，本文选取表2中所示年平均最大风频风场（风速：m/s），同时选取表3中所示5个不利极值风，考虑溢油对近海海域可能造成的最大影响。

4.6 预测结果

（1）油膜漂移轨迹

图8、9为崖城13-4气田主体区发生溢油后，分别在平均风和极值风（见表2、3所示风场）下所模拟出油膜的漂移轨迹示意图。

图8 平均风下油膜漂移轨迹示意图

图9 极值风下油膜漂移轨迹示意图

表2 年平均风场

表3 年极值风场

（2）油膜漂移的平均速率、抵岸时间及抵达功能区的时间

通过计算可知，无论在什么时候发生溢油，油膜经过一定时间会溢出计算域或者抵达岸边。如果溢油量偏少并且油膜的漂移时间较长，油膜会在溢出计算域前或抵岸前消失。表4、5给出平均风和极值风况与流场耦合的条件下，在崖城13-4气田主体区发生溢油后，油膜漂移的距离，平均速率、抵岸时间和抵达重点海域功能区的时间等。由表5可以看出，发生溢油之后，在出现极值不利风向下，油膜会在较短时间内抵达需要重点保护的海南岛南部重点海域功能区（最短为偏西南风向下，21.4 h即可到达敏感海域，对其造成影响）。因此需要各方在溢油事故发生后及时采取有效措施对溢油进行拦截。

表4 平均风下溢油不同条件下油膜漂移时间等

表5 极值风下溢油不同条件下油膜漂移时间等

4 结语

本文采用FVCOM三维潮流模式模拟海南岛南部海域的潮流场，看出崖城油气田附近海域的潮流表现为往复流，涨潮为西北向流，落潮为东南向流，并且落潮潮流大于涨潮潮流。在潮流模拟正确的情况下进行了崖城附近海域溢油轨迹的模拟，并预测了油膜抵达油气田附近重点海域功能区的时间。在出现极值不利风向下，油膜会在较短时间内抵达需要重点保护的海南岛南部重点海域功能区（最短为偏西南风向下，21.4 h即可到达敏感海域，对其造成影响），油气田附近一旦发生溢油，应立即采取行之有效的应急措施，保护溢油可能到达的重点保护海域，防止其遭受污染。

[1]万修全,鲍献文,吴德兴,等.胶州湾及其邻近海域潮流和污染物扩散的数值模拟[J].海洋科学,2003,27(5)：31-36.

[2]徐振华,雷方辉,娄安刚,等.北部湾潮汐潮流的数值模拟[J].海洋科学,2010,34(2)：10-14.

[3]赵昌,吕新刚,乔方利,等.北部湾潮波数值研究[J].海洋学报, 2010,32(4)：1-11.

[4]曹德明,方国洪.北部湾潮汐和潮流的数值模拟[J].海洋与湖沼, 1990,21(2)：105-113.

[5]孙洪亮,黄卫民.北部湾潮汐潮流的三维数值模拟[J].海洋学报, 2001,23(2)：1-8.

[6]叶安乐.分潮流最大流速发生时刻随深度变化的特征[J].海洋通报,1984,3(2)：1-8.

[7]娄安刚,王学昌,孙长青,等.胶州湾海面溢油轨迹的数值模拟[J].黄渤海海洋,2001,19(1)：1-8.

[8]娄安刚,吴德星,王学昌,等.三维海洋溢油预测模型的建立[J].青岛海洋大学学报：自然科学版,2001,31(4)：473-479.

The numerical simulation and forecast of spilled oil trajectory in Yacheng oil and gas field sea area

WANG Jing,LOUAn-gang,CAO Zhen-dong
(College of Environmental Science&Engineering,Ocean University of China,Qingdao,266100 China)

Based on the unstructured grid finite-volume coastal ocean model(FVCOM),the complete three-dimensional fine-solution equation group of shallow-sea tide is used to compute the tidal current and sea-level variations in Yacheng sea area.The open boundary conditions are provided by the open oceanic model.The simulation of the tidal current and sea-level are verified with the observation data.By Euler-Lagarangian tracing method,a trajectory forecasting model for spilled oil is established to predict the trajectory of the oil-film centre,the averaging moving rate and the reaching time of the oil to the sensitive area in Yacheng oil and gas field sea area. It offers some effective technologies for the emergency response to the oil and gas field.

FVCOM model;Yacheng;oil spill;numerical simulation

book=19,ebook=19

X55

：A

：1003-0239（2012）03-0073-08

2011-05-09

王璟（1987-），男，硕士研究生，主要从事海洋环境动力和数值模拟研究。E-mail：qiaokeli16520245@163.com