绿潮是海洋大型藻类爆发性生长聚集形成的藻华现象。近年来, 绿潮在全球沿岸海域爆发并造成严重危害的现象变得越来越频繁, 发生地理范围也日趋增大, 已经成为一种世界性的海洋灾害^[1], 其监测与防灾减灾一直是世界范围内的难题。从2007年黄海海域首次观测到漂浮浒苔开始, 浒苔绿潮每年都会周期性的在黄海海域出现, 影响着江苏、山东沿海地区。国家海洋局发布的《中国海洋灾害公报》显示, 2008年5—8月, 黄海海域暴发的浒苔绿潮灾害造成直接经济损失达13.22亿元。

江苏海域是绿潮灾害的多发海域, 2007年至今, 江苏近岸海域连续多年发生绿潮灾害, 并具有分布面积广、覆盖范围大、持续时间长的特点。绿潮发生对江苏近岸海域生态系统产生直接影响, 破坏海洋环境并间接对渔业生产造成影响^[2]。有研究表明, 黄海绿潮爆发的源头就在江苏近岸海域, 紫菜养殖筏架是绿潮暴发的重要源头, 紫菜养殖筏架回收过程中, 数千吨绿潮藻被刮落入海是绿潮暴发的关键因素^[3]。

如今, 在绿潮灾害成为一种周期性现象的同时, 更为重要的问题是:如何预防这一灾害的发生, 在灾害发生后如何快速地作出反应。卫星遥感在监测绿潮方面具有重要的作用, 是获取绿潮信息的有效手段^[4-9], 但遥感观测也有局限性, 当被监测海面上空云量较多时, 可能会出现连续几天不能获得绿潮信息的问题, 而且遥感不能对绿潮的路径做出定量的预测。数值模拟可以对海洋现象进行连续模拟, 在有良好资料的基础上, 利用数值模拟来预测绿潮路径可以补充遥感手段的不足^[10-11]。

本文通过收集江苏海域水深地形资料、潮汐潮流资料、风场资料及卫星遥感绿潮监测资料等, 建立了江苏重点海域绿潮漂移扩散模型, 可以为该海域绿潮灾害的预警预报提供有益参考。

1 模式简介

MIKE3^[12]是丹麦水利研究所(Danish Hydraulic Instititute, DHI)开发的专业三维自由水面流动模拟系统, 适用于湖泊、河口、海湾和海岸地区的水动力及其相关现象的三维仿真模拟, 可进行水利港口工程设计及规划、复杂条件下的水流计算、泥沙沉积与传输、水质模拟预报和环境治理规划等多方面研究应用^[13-16]。

本文利用MIKE3 FLOW MODEL FM模型, 建立三维潮流数学模型, 主要对区域的流场进行模拟。模型基本方程如下:

(1) 连续方程:

$ \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} + \frac{{\partial w}}{{\partial z}} = S $

(1)

其中, u、v、w:水流速度在x、y、z方向的分量; x、y、z:笛卡尔右手坐标系下的坐标; S:源汇项。

(2) 动量方程:

$ \frac{{\partial u}}{{\partial t}} + \frac{{\partial {u^2}}}{{\partial x}} + \frac{{\partial vu}}{{\partial y}} + \frac{{\partial wu}}{{\partial z}} = fv - g\frac{{\partial \eta }}{{\partial x}} - \frac{1}{{{\rho _0}}}\frac{{\partial {p_a}}}{{\partial x}} - \frac{g}{{{\rho _0}}}\int\limits_z^\eta {\frac{{\partial \rho }}{{\partial x}}} dz \\\;- \frac{1}{{{\rho _0}h}}\left( {\frac{{\partial {s_{xx}}}}{{\partial x}} + \frac{{\partial {x_{xy}}}}{{\partial y}}} \right) + {F_u} + \frac{\partial }{{\partial z}}\left( {{v_t}\frac{{\partial u}}{{\partial z}}} \right) + {u_s}S $

(2)

$ \frac{{\partial v}}{{\partial t}} + \frac{{\partial {v^2}}}{{\partial x}} + \frac{{\partial uv}}{{\partial y}} + \frac{{\partial wv}}{{\partial z}} = - fu - g\frac{{\partial \eta }}{{\partial y}} - \frac{1}{{{\rho _0}}}\frac{{\partial {p_a}}}{{\partial y}} - \frac{g}{{{\rho _0}}}\int\limits_z^\eta {\frac{{\partial \rho }}{{\partial y}}} dz \\\; - \frac{1}{{{\rho _0}h}}\left( {\frac{{\partial {s_{yx}}}}{{\partial x}} + \frac{{\partial {x_{yy}}}}{{\partial y}}} \right) + {F_v} + \frac{\partial }{{\partial z}}\left( {{v_t}\frac{{\partial v}}{{\partial z}}} \right) + {v_s}S $

(3)

其中, t:时间; η:表面水位; d:静水深; h:总水深, 满足关系式h=η+d; f:科氏力参数, 满足关系式f= 2ΩsinΦ, Ω为地球自转角速度, Φ是维度; g:重力加速度; ρ、ρ₀:水体密度和海水密度; S_xx、S_xy、S_yx、S_yy:辐射应力分量; v_t:垂向紊流涡粘系数; P_a:大气压力; u_s、v_s:源汇项水流速度; F_u、F_v:斜压状态下的水平压力项。

2 水动力模型配置及验证 2.1 模型范围和网格

数值模型网格南起浙江象山, 北至山东威海, 基本覆盖了研究所关注的区域。模型水平方向采用三角形网格, 并在苏北浅滩地形变化梯度比较大的区域适度加密网格(图 1)。整个计算区域共有89 155个三角形单元, 45 193个节点, 网格最小空间步长约为200 m。模型垂向采用σ坐标, 共均匀分为6个σ层。

计算域内地形方面, 苏北浅滩以及近海的水深采用中国人民解放军海军司令部航海保证部出版的2013年版海图上的水深数据, 外海的水深数据采用美国国家地球物理中心(National Geophysical Data Center, NGDC)提供的ETOP01全球地貌数据, 空间分辨率为1′×1′。水深基准面统一为平均海平面, 计算区域水深如图 1所示。

2.2 定解条件 2.2.1 初始条件

由于潮波运动是一种摩阻运动, 故采用冷启动, 即潮位为零, 流速为零, 由此产生的误差在计算过程中会自行消除, 在需要模拟的过程提前5 d开始计算, 以使模型稳定。模型为正压模型, 不考虑温度和盐度影响。

2.2.2 边界条件

闭边界采用不可入条件, 法向流速为零, 即$ \vec V \cdot \vec n = 0 $。

外海开边界共设置40个边界点, 给定潮位, 由M₂、S₂、N₂、K₂、K₁、O₁、P₁、Q₁八个主要分潮的调和常数通过计算合成。分潮调和常数采用日本国家天文台NAO99b潮波预报模型结果内插, 并经调整后用于本模型。

$ \xi = {A_0} + \sum\limits_{i = 1}^8 {{H_i}} {F_i}\cos \left[ {{\sigma _{it}} - {{\left( {{v_0} + u} \right)}_i} + {g_i}} \right] $

(4)

式中, A₀为平均海平面, F_i、(v₀+u)_i为天文要素, H_i、g_i为某分潮的调和常数, 即振幅与迟角。

模型考虑了风场对海表流速的影响, 海表面风场采用美国环境预报中心和国家大气研究中心联合推出的NCEP/NCAR Reanalysis 1再分析数据。

2.3 模型参数 2.3.1 时间步长

模型主时间步长取30 s, 根据CFL数来调整实际计算步长以保证模型稳定, 最小时间步长取0.01 s, CFL条件数取为0.8。

2.3.2 干湿水深

模型采用干湿判别法处理动边界, 即当网格水深小于0.1 m时, 设动量为零, 只考虑物质通量; 当网格水深小于0.005 m时, 该网格不参与计算。

2.3.3 涡粘系数

水平涡粘系数采用Smagorinsky公式, Smagorinsky系数C_s为0.28;垂向涡粘系数采用k-ε方程求解。

2.3.4 径流

模型考虑长江径流影响, 长江径流量取大通水文站流量数据。

2.4 模型验证

利用收集的近年来吕四、如东、大丰、滨海、连云港海域10个站位的实测海流资料以及吕四、洋口、连云港3个站位实测潮位资料对模型结果进行验证, 模型验证点位基本覆盖了所关心的区域, 验证资料站位分布见图 2。

2.4.1 潮位验证

潮位验证结果如图 3~图 5, 由图可见, 江苏沿海三个站点的潮位计算值与实测潮位较吻合, 其中计算值和实测值的误差小于10%, 计算潮位与实测潮位相位基本一致。

对模拟潮位进行调和分析, 得到M₂、S₂两个主要分潮的同潮图(图 6), 其中虚线为等振幅线, 实线为等迟角线。与海洋图集上的分潮同潮图相比, 两个分潮的振幅和迟角拟合效果都较好, 分潮无潮点位置基本一致。

2.4.2 海流验证

吕四、如东、大丰、滨海和连云港10个站位的表层海流验证结果如图 7~图 11, 由图可见, 模型计算流速与实测流速值吻合较好, 涨落急流速和流速变化过程较为一致, 流向的模拟结果与实测值也具有良好的一致性。

通过以上模型的率定与验证表明:水动力模型采用的物理参数和计算参数基本合理, 计算方法可靠, 能够较好地模拟研究潮波运动和水动力场特征。水动力模型可以为绿潮漂移扩散计算提供基本的水动力场输入条件。

3 绿潮漂移模型建立

在不考虑绿潮藻自身生态过程的情况下, 绿潮藻在海水中的移动可以看作是质点在表层海流和海面风共同作用下的漂移、输运过程, 可以采用拉格朗日粒子追踪方法对绿潮藻的漂移路径进行预测。

3.1 计算方法

绿潮漂移计算采用欧拉-拉格郎日非耦合法, 即利用水动力场的计算结果, 结合公式(5)和公式(6)计算绿藻质点漂移情况。

$ X = {X_0} + \left( {U + \alpha {W_{10}}{\rm{cos}}A} \right)\Delta t + R{\rm{cos}}B $

(5)

$ Y = {Y_0} + \left( {V + \alpha {W_{10}}{\rm{sin}}A} \right)\Delta t + R{\rm{sin}}B $

(6)

其中, X₀、Y₀、X、Y分别代表某绿潮质点的前后时刻的坐标; U、V分别为X、Y方向的流速分量; W₁₀为海面上的风速; A为风向; α为风拖曳系数, α的取值范围一般在0.01~0.05^[17-18], 本文结合对海上漂流浮标漂移轨迹模拟的多次数值试验, 确定α取值为0.04; R为随机扩散项, R=rE, r为0~1之间的随机数, E为扩散距离, $ E = \sqrt {12D\Delta t} $, D为水平方向扩散系数, 其量值等于水平涡粘系数, 由水动力模型计算获取; B为随机扩散方向, B=2πr; Δt为时间步长, 与水动力模型同步计算。

其方法步骤: 1、确定初始坐标位置; 2、根据流场计算结果, 用插值公式求该点的流速值, 代入上式求下一时刻的坐标位置。以此方式循环递推计算绿潮质点漂移轨迹。

3.2 绿潮漂移扩散模拟

根据国家海洋局东海监测中心发布的《东海监测中心卫星遥感监视监测简报》中绿潮遥感监测信息, 对比验证漂移模型后报结果。由于绿潮(浒苔)的繁殖扩展速度非常之快, 有时相差一天的卫星图片表现出来的绿潮分布形状大为不同, 很难辨别出前一天的绿潮在今天的位置。本文经过认真筛选后, 选取2014年5月21日—23日一次较为完整的绿潮过程, 可以明显地辨别出同一片绿潮藻的漂移扩散状况, 作为绿藻漂移实测位置, 与模拟结果作对比检验。

根据卫星遥感监测结果, 提取绿潮漂移信息, 本次绿潮过程为:

(1) 5月21日13时MODIS影像发现, 在盐城至连云港市海域分布有浒苔绿潮, 分布范围: 120°32′~ 121°11′E、34°12′~35°09′N, 覆盖面积69.4 km², 分布面积1 547 km²。

(2) 5月22日13时MODIS影像发现, 在盐城至日照市海域分布有浒苔绿潮, 分布范围: 120°27′~ 121°27′E、34°12′~35°29′N, 覆盖面积137.4 km², 分布面积3 928 km²。

(3) 5月23日12时45分MODIS影像发现, 在盐城至日照市海域分布有浒苔绿潮, 分布范围: 120°12′~ 121°35′E、34°16′~35°36′N, 覆盖面积181.6 km², 分布面积5427 km²。

通过遥感监测结果可以看出, 本次过程, 浒苔绿潮向东北方向漂移扩散, 从盐城-连云港海域逐渐漂向日照海域; 绿潮面积呈扩大趋势, 5月21日分布面积为1 547 km², 5月22日分布面积为3 928 km², 5月23日分布面积为5 427 km²。

利用绿潮漂移扩散模型对绿潮漂移扩散情况进行模拟:绿潮藻的运移主要受表层海流和海表面风的共同影响, 因此, 设置绿潮质点于海表面。将绿潮质点释放的初始时刻设置为5月21日13时, 初始位置设置为卫星遥感所监测到的绿潮分布位置(图 12), 在绿潮初始位置范围内一次性均匀释放1 600个粒子。输入水动力模型模拟海流和风场条件, 启动模型进行积分计算; 预测24 h和48 h后绿潮的分布位置。

24 h预测结果如图 13所示, 绿潮向东北方向漂移扩散, 分布范围: 120°30′~121°19′E、34°21′~35°29′N, 分布面积约为3 675 km²。48 h预测结果如图 14所示, 绿潮继续向东北方向漂移扩散, 分布范围: 120°30′~ 121°30′E、34°24′~35°36′N, 分布面积4 713 km²。

将绿潮24 h和48 h后的模拟结果与卫星遥感结果进行比较(表 1), 结果表明绿潮数值模拟结果与遥感监测结果较为一致: 24 h绿潮模拟分布面积为3 675 km², 模拟分布区域与遥感监测结果位置重合率为63%; 48 h绿潮模拟分布面积为4 713 km², 模拟分布区域与遥感结果位置重合率为61%。这表明模型能够较准确预测绿潮斑块的漂移轨迹。

表 2同时给出了风拖曳系数α取不同值时的模拟结果, 可以看出, α的取值对绿潮模拟分布面积和位置均有较大影响; 当α取值过小时, 绿潮粒子扩散不充分, 导致模拟绿潮分布面积偏小, 且分布范围也与卫星遥感结果存在较大偏差; α取值过大, 会使模拟绿潮分布位置在风的作用下出现较大偏移, 降低模拟的精度。

4 结论

本文通过收集江苏海域水深地形资料、潮汐潮流资料、风场资料及卫星遥感绿潮监视资料等, 建立了江苏重点海域水动力数值模型并进行了模型验证。验证结果表明, 水动力模型能够较好地模拟本海域潮波运动和水动力场特征, 可以为绿潮漂移扩散计算提供基本的水动力场输入条件。在水动力模型的基础上建立了绿潮漂移扩散模型, 利用该模型对绿潮漂移扩散过程进行了模拟; 通过与卫星遥感监测结果相对比, 漂移扩散模型对绿潮24 h和48 h的漂移路径和扩散范围的模拟精度较高, 能为江苏海域绿潮的预警预报提供一定参考。