平潭竹屿湾水交换能力和溢油风险数值模拟

张振伟，刘必劲*，曾银东，张世民

(1. 厦门理工学院土木工程与建筑学院，福建 厦门 361024; 2. 福建省海洋预报台，福建 福州 350503；3. 国家海洋局厦门海洋预报台，福建 厦门 361008)

平潭竹屿湾(图1)位于平潭综合实验区中部，北靠霞屿村鹅头山(临近城市核心区)，东接平潭旧城区，南与敖东镇及南寨山接壤，西南接竹屿出海口。国家相关部门支持在平潭竹屿湾开展生态整治和修复工作，工作内容主要包括海湾深度清淤疏浚、海绵污染防控系统建设、竹屿湾水环境监控能力建设等。生态整治和修复工作实施后，可以有效提高平潭防灾减灾能力、改善竹屿湾水环境等。本次数模主要针对竹屿湾水交换能力和溢油风险影响范围进行研究，为竹屿湾的水环境监控能力建设及管理提供支撑。

有限体积海岸海洋模型(Finite Volume Coast Ocean Model, FVCOM)模型是应用有限体积法对守恒形式的三维水动力方程进行离散而构建的海洋动力计算模型[1-2]。模型基于非结构三角形网格进行离散，这种网格可以较好的拟合不规则水域边界。许多学者开展了海域水体交换和溢油等方面的研究。赵亮等(2002)应用水动力模型(Estuary, Coast Ocean Model, ECOM)基于水质点运移的物理交换，针对胶州湾的水交换能力开展了计算研究[3]。李真(2009)应用FVCOM模型对围垦引起的罗源湾的水动力环境进行了模拟研究[4]。王兴刚等(2015)开展了连云港港主体港区水交换能力的研究[5]。迟万清等(2018)研究了不同风向作用下胶州湾水体交换规律[6]。王雪等(2017)应用FVCOM模型研究了胶州湾各个子区域之间的水交换情况[7]。刘晓东等(2017)应用FVCOM模型开展了厦门围头湾海域溢油风险计算分析[8]。

1 材料与方法

1.1 FVCOM模型

为了提高浅水区域的垂向分辨率，FVCOM模型垂向坐标应用了σ坐标变换。FVCOM模型中σ坐标有多种形式可以选择，本次模拟中我们采用了垂向均匀分布的σ坐标。FVCOM模型的动量方程和连续方程表达为：

(1)

(2)

(3)

水平耗散项表达为：

(4)

(5)

(6)

式(1)至(6)中：x、y、σ为三维笛卡尔直角坐标系的正东(m)、正北(m)和垂直方向坐标；u、v、ω为x、y、σ方向上的速度分量(m/s)；ζ为自由水面(m)；D为总水深(m)；H为静水水深(m)；q为湍动能(m/s)；l为湍流宏观尺度(m)；Fx、Fy分别为动量x方向、动量y方向的水平扩散项(m/s2)；ρ为水体密度(kg/m3)；f为科氏力系数(s-1)； g为重力加速度(m/s2)；Km为垂向涡粘系数(m2/s)；Kh为垂向热扩散系数(m2/s)；Am为水平涡粘系数(m2/s)，Ah为水平热扩散系数(m2/s)；它们可以由修正的MY-2.5湍流闭合子模型计算[9-10]。

图1 平潭竹屿湾及其邻近海域水深地形和计算网格Fig. 1 Bathymetry of Zhuyu Bay and its adjacent waters and grids of the model

1.2 污染扩散模型

FVCOM模型中污染扩散方程表达为：

(7)

式(7)中：C为物质浓度，KH为垂向涡粘系数；Fc为水平扩散项；C0为点源浓度。

(8)

式(8)中：ts、te分别为物质释放的开始和结束时间；i为三角网格节点编号，N为物质释放节点的总数。

水交换能力反映了水体的自净能力，与水环境、水生态密切相关。评价水体交换能力的指标有很多，如平均滞留时间[11]、半交换时间[12]、冲洗时间等。平均滞留时间是指水质点离开研究区所需的平均时间。半交换时间是指保守物质浓度降至初始浓度一半时所需时间。冲洗时间是指研究区内保守物质浓度降为零时所需时间。假定t0时研究区域内某保守物质的浓度为C0，而后某时刻t(t>t0)存留于研究区域内的物质浓度变为C(t)，保守物质平均滞留时间(θ)定义为：

(9)

式(9)中：当C(t)=50%C0时，Δt=t-t0即为水体的半交换时间；当C(t)=0时，Δt=t-t0即为水体的冲洗时间。

1.3 粒子追踪模型

拉格朗日粒子追踪模型，可以模拟粒子在风、波浪、海流等作用下的漂移运动轨迹，假设粒子不会发生碰撞，也不会发生混合，则粒子扩散方程写为：

(10)

u=uc+κ·uw+ur，v=vc+κ·vw+vr

(11)

式(10)、(11)中：原坐标为(x0，y0)的粒子经时间Δt=t-t0后，漂移到坐标(x,y)；u和v分别是粒子运动的东、北分量，它由流速uc、风速uw、粒子随机运动速度ur组成，κ为风对粒子拖曳系数，取0.025。通过跟踪各粒子坐标(x,y)的各位置，确定运移范围，统计分析粒子扫海面积。

2 结果与讨论

2.1 潮汐潮流计算结果

本研究的计算区域范围为24.82°—25.90°N，119.10°—120.23°E。为保证模拟精度，海域水深根据中华人民共和国海事局出版的最新海图及竹屿湾附近海域实测水深插值得到。模型的水深地形和网格如图1所示。在开边界处模型计算网格分辨率约为4 km，在近海岸处网格进行了加密，竹屿湾海域网格分辨率约为30 m左右。模型垂直方向采用均匀分布的符号坐标，为与实测层数一致垂向分为7个σ层，摩擦系数选用ORIG选项进行计算。

模型开边界的潮位时间序列从OTIS全球潮汐模型 (OSU Tidal Inversion Software)数据库中提取，在生成潮位时间序列时选择了M2、S2、K2、N2、K1、O1、P1、Q1八个主要分潮。

模型模拟了时间自9月1日—14日的潮汐潮流运动情况，该时间范围包含了大潮和小潮的观测期限，小潮观测时间(UTC)9月4日03时—5日03时。大潮观测时间(UTC)自9月12日03时—13日03时，水文调查站位见图2。图3给出了平潭海洋站大潮测流期间潮位模拟与实测值的对比，由图可知潮位时间序列计算结果与实测结果符合较好。

图2 平潭竹屿湾及其邻近海域水文调查站点位Fig. 2 Locations of tidal current observation stations in Zhuyu Bay and its adjacent waters, Pingtan

图3 平潭竹屿湾2018年9月潮位观测结果与模拟结果对比Fig. 3 Comparison of observed and modeled water elevations of Zhuyu Bay in September 2018

图4 大潮和小潮流速、流向验证结果Fig. 4 Modeled current validation against observations during spring and neap tides

图4给出了大潮流速、流向和小潮流速、流向验证结果，由图中结果可以看出，计算流速、流向总体上与实测值吻合较好。海坛海峡内潮波属于前进波，潮位越高流速越大。高潮位时海峡内海流自北向南流，流速大；落潮半潮位时海流在分流尾屿区域形成南北分流，分流尾屿北侧海流向北流，而南侧海流则向南流，海峡中部流速较小；低潮位时海峡内海流自南向北流，水道上流速较大。因潮差大，海峡两岸露滩面积大。自低潮再经过约3 h，至涨潮半潮位时，海峡北东口海流自北向南流与海峡北上涨潮海流在分流尾屿区域汇合，海峡中部流速较小[13]。潮位潮流验证结果表明本研究建立的数值模型，能够很好的模拟研究区内潮位、潮流运动特征，可以应用模型开展区域水环境特征的计算研究。

图5给出了大潮高潮竹屿湾内潮流场的矢量图。由图可以看出，竹屿湾内的潮流流速在大潮高潮落潮半潮位和高潮涨潮半潮位时的流速最大，流速超过1.0 m/s，说明潮流动力很强。在大潮高潮涨平和落平时刻流速小，最小流速量值在0.3 m/s左右。因竹屿湾呈狭长状态，分析竹屿湾中特征点的流速流向过程线可以看出潮流呈往复流特征，潮汐呈驻波特征。

图5 平潭竹屿湾及其邻近海域大潮水深平均潮流场Fig. 5 Tidal current velocity vector and elevation in Zhuyu Bay and its adjacent waters

2.2 水体交换计算结果

利用建立的竹屿湾水动力模型，构建竹屿湾对流扩散模型，进行竹屿湾水交换能力指标的计算。计算时以竹屿湾口为界，在湾内释放可溶于水的保守物质，湾外物质浓度设为0.0 mg/L，湾内物质浓度设为1.0 mg/L，释放时间选为小潮高潮时刻(UTC时间9月3日08时)。计算水半交换时间是对于计算域内每个网格，当网格节点的浓度第一次小于0.5 mg/L时，记下该时刻，这一时刻与保守物质释放时刻的差即为水半交换时间。对于平均滞留时间，需要按照式(9)进行积分计算，对于每个网格节点，记录其保守物质第一次变为0.0 mg/L的时刻，积分的时间为保守物质释放时刻至保守物质变为0.0 mg/L的时刻，这样可以得到每个网格节点的平均滞留时刻。对于冲洗时间，每个网格节点也是不一样的，这里没有针对每一个节点开展计算，而是取竹屿湾交换能力最差的位置处的冲洗时间作为整个海湾的冲洗时间。图6、7给出了竹屿湾水体半交换时间和平均滞留时间的分布情况，由图可知，竹屿湾大部分水域水体半交换时间小于1.0 d，平均滞留时间3.0 d左右。竹屿湾全部水体的冲洗时间为15.0 d。总体上，竹屿湾水体交换能力较强。

Liu 等(2004)针对不同释放时间对平均滞留时间造成的误差做了分析[14]，其分析结果证明平均滞留时间与初始浓度无关，但保守物质的释放时刻不同，对平均滞留时间的计算结果会造成影响，其影响范围在5%～10%之间。

2.3 溢油扩散计算结果

竹屿湾溢油事故风险主要来自于海域的运输船只，假定运输船只发生溢油事故，3 h内漏油5 t，模拟泄露油品为燃料油。溢油点选在竹屿湾湾口航道区域。采用瞬时排放的方式，计算时不考虑粘结，也不考虑挥发，模拟溢油点位于竹屿湾口水道，坐标为(25.515°N，119.697°E)，油粒子位置每隔10 min输出一次。

溢油数模旨在了解一旦发生溢油事故，溢油的扩散分布趋势。本次溢油模拟主要模拟静风和不利风向(NNE风向，风速9.2 m/s)条件下，高平潮、低平潮、涨急、落急等典型潮时发生溢油时，溢油泄漏48 h内各典型时刻的影响范围及到达附近环境敏感点的时间。海坛海峡及附近区域敏感目标包括海坛岛北部保留区、福清湾农渔业区、山岐澳中华鲎保护区、塘屿列岛海洋保护区等。

图6 竹屿湾水体半交换时间分布Fig. 6 Distribution of half-life time in Zhuyu Bay

图7 竹屿湾平均滞留时间分布Fig. 7 Distribution of average residence time in Zhuyu Bay

表1给出了静风和不利风时刻，油粒子扫海面积的统计结果，油粒子在某一单元内出现，则该单元的面积作为油粒子扫海面积的贡献值。图8给出了不利风条件下，油粒子48 h扩散范围图。

表1 平潭竹屿湾溢油扫海面积Tab. 1 Area of oil film after oil spill over 48 h in Zhuyu Bay

图8 不利风向下竹屿湾及其邻近海域4个不同潮时刻溢油48 h轨迹分布Fig. 8 Distribution of oil particles trajectory after oil spill over 48h at 4 tide moments under the averse wind conditions in Zhuyu Bay and its adjacent waters

计算分析表明静风条件下，油粒子总体上与潮流运动趋势一致。涨急、高平潮和低潮情况下，油粒子主要呈现出在海坛海峡内的往复运动，总体上油粒子在靠近海坛海峡平潭一侧运动。落急时刻，油粒子一旦从竹屿湾口水道运动至海坛海峡水道，油粒子会迅速向海坛海峡南口区域运动，其扩散范围达到了平潭坛南湾、平潭塘屿岛、平潭草屿岛、高山湾等区域，其面积可达40.800 km2。

由图8可知，不利风对油粒子扩散影响显著。不利风条件下，油粒子在NNE风向的作用下，显现出向海坛海峡福清一侧运动的趋势。油粒子的扩散范围可以到达平潭坛南湾、平潭塘屿岛、平潭草屿岛、高山湾等区域。在涨急、高平潮和落急时刻，油粒子的扫海面积均超过了20 km2，落急时刻扫海面积最大，达到了43.100 km2。低平潮时刻释放的油粒子扩散范围最小，其扫海面积为6.880 km2。因此，不利风对油粒子的漂移扩散影响显著，风的作用扩大了油粒子的影响范围。

3 结论

本研究应用FVCOM海洋动力模型开展了平潭竹屿湾水体交换和溢油扩散计算分析。通过在小潮高潮时刻释放示踪剂，对竹屿湾水体交换能力开展了研究，计算结果表明竹屿湾大部分水域的水半交换时间小于1.0 d，最大3.5 d；平均滞留时间大部分区域小于3.0 d，最大4.0 d；水体冲洗时间为15.0 d。

溢油扩散计算结果表明，静风条件下落急时刻发生溢油时48 h油粒子扫海面积最大，达到了40.800 km2。油粒子的扩散范围到达了平潭坛南湾、塘屿岛、草屿岛及高山湾等区域，对周围海域影响大。其他时刻发生溢油时油粒子主要在海坛海峡平潭一侧往复活动。不利风(NNE向)作用下，溢油粒子扫海面积明显比静风条件大，不同时刻发生溢油时油粒子均可以扩散到平潭坛南湾、塘屿岛、草屿岛及高山湾附近海域，低平潮扫海面积最小6.880 km2，落急时刻发生溢时油粒子扫海面积最大，为43.100 km2。其他时刻发生溢油时，油粒子扫海面积均大于24.200 km2。