海洋是自然以及人类最宝贵的财富, 是生命初始的摇篮, 伴随着人类对海洋的开发利用, 海洋环境问题随之产生。其中, 海洋塑料污染是各类海洋环境问题中最突出的问题之一(王金平等, 2021)。据估计全球约80%的海洋垃圾由塑料构成, 塑料垃圾在进入海洋后, 在光照、波浪冲刷、生物作用等影响下破碎成包括微塑料在内的塑料颗粒(Barnes et al, 2009)。这些塑料颗粒会在复杂的海洋动力(风、洋流及海浪等)作用下进行输运迁移, 并广泛分布于海洋水体、沉积物及生物体中, 危害海洋生物的健康, 也能通过食物链富集传输危害人类自身(陈刘等, 2020), 对整个海洋生态系统造成极大威胁, 严重影响海洋经济的发展。然而, 作为新型海洋污染物, 海洋塑料(包括粒径小于5 mm的微塑料)入海通量计算过程中涉及的源汇过程、迁移途径、输运过程等尚未完全清楚, 需进一步研究完善(李道季, 2020)。因此, 研究塑料垃圾在海洋中的运动特性和输移机理, 对确定海洋塑料垃圾的最终分布和归趋有着至关重要的作用, 也是本文的研究重点。

目前, 研究海洋塑料垃圾的输运迁移规律, 主要面临两个关键性的技术挑战: 一是复杂的海洋动力环境。塑料进入海洋后会在风浪流的共同作用下发生分解、絮凝、漂移、沉降和再悬浮等运动, 复杂多变的海洋动力因素给研究海洋塑料颗粒的输运迁移规律造成了不便。二是海洋塑料自身的理化性质及其演变过程对塑料的输运过程有重要影响。在同样的水动力环境下, 不同分布形态(漂浮、悬浮、沉积)的塑料颗粒因不同粒径大小、形状、密度等产生的输运特性差异明显(Besseling et al, 2017)。

耦合海洋水动力模拟和粒子追踪模型的海洋污染物扩散模型是研究表层漂浮塑料颗粒在表层洋流系统、潮流和风等作用下远距离输送的主要研究方法之一。其中, 粒子追踪模型将水体中的塑料颗粒运动看成风和海流作用下发生的拉格朗日漂移, 通过求解拉格朗日方程得到塑料颗粒的运动轨迹和源汇区域。Van Sebille等(2012, 2015)和Lebreton等(2013)基于洋流动力模型的预测结果表明, 在主要海洋盆地的亚热带纬度上形成了五个海洋塑料垃圾积聚区, 并与航次实地调查结果吻合(Law et al, 2010, 2014)。Chubarenko等(2016)计算发现, 对海水表层的漂浮球状微塑料颗粒, 风阻作用导致其传输速度为海流输送速率的3倍。Zhang等(2020)研究发现, 漂浮和悬浮微塑料在长江口及其临近海域会呈现非常不同的输运规律, 与潮流、风阻和微塑料释放点(入海)位置等强相关。在该研究中, 微塑料的物理性质主要根据实地采样数据进行设置, 因此, 未考虑海洋塑料自身的理化性质多样性对其输运迁移规律的影响。Neumann等(2014)的研究表明, 海表风应力会增大北海南部海域漂浮污染物的粒子移动速度和传输范围。针对2011年日本海啸, Maximenko等(2018)通过数值模拟和原型观测发现, 迎风面大的塑料垃圾输运速度比较快, 而迎风面小、相对较重的塑料垃圾则主要在盆地中部汇集。Ter Halle等(2016)指出立方体微塑料在漂浮过程中受风力的影响将发生持续翻转, 而扁平状的微塑料与水面接触面积大, 较难发生翻转。由此可见, 现有海洋塑料颗粒扩散模型假定颗粒大部分漂浮在海洋表面, 主要研究其在风浪流联合作用下的水平输移。

Reisser等(2015)利用多层拖网对北大西洋表层以下5 m的水体进行调查, 每隔0.5 m进行现场采样。经过对样品长度、类型以及实验分析发现, 海水中塑料颗粒的浓度随深度呈指数下降, 且其浓度在5 m左右水层接近于零。Eriksen等(2014)提出塑料碎片最终将沉降到海底。Mountford等(2019)使用欧拉方法建立了三维海洋塑料碎片扩散模型, 模拟结果表明大量的塑料碎片会在海洋中进行垂直输运, 因此分布于不同水层中。

水平和垂直输运主要是针对漂浮和悬浮在海洋中的塑料颗粒, 而对于埋藏在沉积物中的塑料颗粒而言, 由于其与风浪流等水动力因素响应关系较弱, 目前研究主要集中于其空间分布和来源的分析。Barrett等(2020)分析了澳大利亚湾深海沉积物中的塑料含量, 发现与其他深海沉积物相比, 塑料颗粒数量明显更高, 并且沉积物中塑料碎片的数量随着海洋表面塑性颗粒数量的增加以及海底倾斜角的增加而增加。Dodson等(2020)研究了美国弗吉尼亚州和北卡罗来纳州沿海环境中塑料颗粒的分布情况, 发现即使在海洋表面未受到污染的河口和海滩沉积物中同样存在塑料颗粒。Zuo等(2020)以华南珠江口的三个红树林湿地为研究对象, 研究了潮间带红树林生态系统的海洋塑料污染, 发现塑料颗粒的丰度与华南珠江口的人口密度和国内生产总值显著正相关。由此可见, 沉积物中塑料污染水平与其邻近的城市和国家的发展紧密相关, 探究沉积物塑料通量与临近地区社会经济发展程度的关系是当前研究的一个重点方向。

因此, 本文以福建省厦门湾为研究对象, 基于MIKE3的生物仿真模拟技术(Agent Based Modelling, ABM), 建立厦门湾的三维水动力数值模型, 分别研究漂浮、悬浮和沉积塑料颗粒在复杂海洋动力因素下的迁移特性和输运机理, 并结合水动力条件对厦门湾海域塑料颗粒的分布趋势开展分析和预测。该研究对促进我国沿海海洋生态环境的保护有一定现实意义, 同时也为塑料颗粒在海洋中的迁移扩散研究和预测提供理论依据。

福建省厦门湾位于海陆交界的海岸带, 是地球表面最为活跃的自然区域, 也是资源与环境最为优越、具有海陆过渡特点的独立环境体系。并且, 厦门湾海域不仅接收来自厦门海岸的大部分生活污水和工业废水, 还接收来自九龙江流域的农业径流和顺江而下的大量塑料垃圾(Peng et al, 2013), 是当前研究海洋塑料垃圾输运迁移特性的热点区域之一。

1 研究方法

目前, 针对海洋塑料颗粒迁移问题的建模方法主要有拉格朗日和欧拉方法。MIKE3的生物仿真模拟技术(Agent Based Modelling, ABM模块)耦合了粒子追踪模型(PTM, 拉格朗日方法)和生态实验室模型(ECO Lab, 欧拉方法), 可与MIKE 3的水动力模型(HD)进行集成计算, 模拟研究生物个体对水环境参数做出的响应。本文将水动力HD模块与基于拉格朗日算法的ABM模块集成, 用于模拟厦门湾漂浮塑料垃圾在风、径流、潮流等作用下的输移特性; 将其与基于欧拉法的水质生态模块(ECO Lab)耦合, 用于模拟厦门湾悬浮塑料颗粒在海洋中的物理输运过程。针对悬浮塑料颗粒, 本文模型除了能考虑其在水动力作用下的水平扩散输移, 还可考虑其在海洋表面、海底的沉降、沉积等复杂物理过程。

值得一提的是, 本文的主要研究思路为根据塑料颗粒的性质对其进行分类(漂浮或悬浮颗粒), 并根据分类选择不同的建模方法, 从而对塑料颗粒的时空分布特性形成较系统深入的认识。即, 本文塑料颗粒在水体中的漂浮状态是预先知道/假定的, 而它们(漂浮和悬浮塑料颗粒)之间的相互转换或者耦合关系不是本文的研究重点。

1.1 厦门湾水动力模拟

MIKE3的水动力模拟是基于三维不可压缩的雷诺平均纳维-斯托克斯方程(RANS方程)的解, 服从布辛涅司克(Boussinesq)和静水压的假设, 满足连续性方程和动量方程, 而温度和盐度等的迁移则遵循一般的迁移扩散方程。求解域的离散化采用有限体积法, 通过将连续体细分成不重叠的单元/元素来离散化空间域。在三维情况下, 主要使用分层网格划分技术, 即水平域使用非结构化网格, 而在垂直域则使用结构化网格。时间积分是半隐式的, 即水平项被隐式处理, 垂直项被隐式或部分显式和部分隐式处理。具体可参考MIKE3的用户使用手册。

本文的研究区域选取厦门湾及其周围海域, 根据牛敬硕(2019)及谢森扬(2016)对厦门湾水动力特性的研究, 选取T1潮位测站及V1潮流测站的实测数据与模型模拟结果进行对比; 测站的平面分布如图 1所示。模型计算时每隔一小时输出一次结果, 模拟值和实测值的对比结果如图 2、3所示, 其中实线为模型计算值, 散点为实测值。具体的模型设置等参见附录A。

从图 2中可以看出, 潮位整体拟合情况良好, 实测值和模拟值的潮汐相位基本一致, 高低潮发生时刻也基本一致, 最大振幅在大潮时约为2.6 m, 在小潮时约为1.3 m。在一天之内实测值和模拟值均出现两次涨潮和落潮的过程, 体现了所在海域的规则半日潮特征。同时, 潮流测站的模拟值与实测值亦较为吻合(图 3); 流速大小基本一致, 表层流速最大可达0.8 m/s, 潮流变向时刻基本一致。从图中还可以看出潮流测站点所在海域呈现明显的往复流特征, 这在一定程度上也促进了塑料颗粒在近岸和外海之间的传输, 且随着水深的增加最大流速值逐渐减少, 底层最大流速大约为顶层的3/5。由此可见, 模型的模拟结果可以很好地体现厦门湾的潮流特征。

1.2 厦门湾漂浮塑料垃圾的输移模拟

基于验证较好的水动力数值模型, 在厦门湾九龙江口释放漂浮塑料垃圾(拉格朗日粒子), 通过MIKE3的粒子追踪模型捕捉其在厦门湾的迁移轨迹, 并将其与实测值进行对比, 验证所建立模型的准确性。厦门湾塑料漂移轨迹的实测值由2020年8月1日《九龙江口-厦门湾海洋垃圾漂移轨迹及分布预测预报单》提供(厦门市生态环境局, 2022)。厦门大学自2017年以来开展了九龙江口-厦门湾的海上垃圾监测和预报工作, 该报告单中提供了每日6时至12时从九龙江口漂浮到厦门湾的海上垃圾轨迹信息。

表面风和流是漂浮粒子物理迁移的主要驱动力, 相关信息(如流场)由上文所建立的水动力模型提供。其中, 风场除了会对水流产生影响(即风生流), 漂浮在自由水面的(较大)粒子, 如上述实测轨迹中的较大型海洋塑料垃圾, 还会受到风的直接拖曳作用, 即风阻力。Critchell等(2016)提出可以使用风飘移系数来计算表面风速中有多少被转移到粒子运动速度中, 即风漂移系数的取值与粒子大小、迎风面等因素密切相关。同时, 可通过设置风飘移角来考虑风阻力的实际作用方向, 即考虑科式力的影响。通过校准, 本文设置风漂移系数为0.01, 风飘移角度为28°。

此外, 粒子在水体中的沉降速度也是影响其运动特征的重要物理量, 其取值同样与粒子的大小、形状、密度等物理性质密切相关。由于验证模型的研究对象为始终漂浮于自由水面的塑料垃圾, 因此, 此处假定粒子的密度始终小于海水密度, 设定其沉降速度为−1 m/s, 忽略由于塑料密度与海水密度不同引起的垂直运动。这也符合本文关于漂浮塑料颗粒的定义, 即始终漂浮于自由水面的塑料颗粒, 而塑料颗粒的垂直运动特性则主要通过下节研究悬浮塑料垃圾的运移来明确。

图 4给出了海上塑料垃圾漂浮轨迹的实测值与模拟值的对比: 图 4a中绿色的线是当日海上塑料垃圾的实际漂移轨迹, 而模拟轨迹(图 4b)中绿色轨迹线是在相同时间段和风况下, 本文模型在九龙江口释放的漂浮塑料粒子的预测运动轨迹(九龙江口污染点源的具体位置参见附录C)。通过对比漂浮塑料颗粒群的实测和模型的迁移路径和长度可以看出, 实测值和模拟值较为接近, 本文漂浮塑料颗粒模型的精度较好。

1.3 厦门湾悬浮塑料垃圾的输移模拟

Koutnik等(2021)将海洋中悬浮塑料颗粒的物理运输过程简化为大气沉降、对流扩散、沉积及再悬浮等四个过程。悬浮塑料颗粒在水体中的对流扩散可通过建立物质输运扩散模型(如MIKE3的对流-扩散模块)进行模拟。然而, 为了同时考虑其他过程和相互作用, 本文将水动力模块和Eco Lab生态模块相耦合, 模拟悬浮塑料粒子浓度的时空分布, 并结合实测资料对所建模型进行验证。

在Eco Lab中, 状态变量(如悬浮塑料粒子的浓度)的输运方程为

其中, C为Eco Lab水质状态变量的浓度。u, v, w分别为流速在x, y, z方向的分量, D_x, D_y, D_z分别为x, y, z方向的扩散系数, S_c为源和汇。P_c为Eco Lab过程, 包括上文所述的塑料粒子沉降和再悬浮等。多个过程之间以线性或者非线性的方式进行耦合,

其中, P_i为与微塑料迁移运动相关的物理过程, n为特定状态变量所包括的过程数, 本文n = 3。描述这些过程引起的塑料粒子浓度变化率的数学表达参见附录B。

图 5给出了模型计算结果与厦门湾特定站点处塑料采样分析得出的浓度的垂直和水平分布对比(Tang et al, 2018), 测站具体信息参见附录C。浓度垂直分布为模拟时刻开始第16 d的计算结果, 而水平分布对比中每个站点的水平浓度为模型在开始模拟后20 d内结果的平均值。从图 5可以看出, 从海面往下塑料粒子浓度近似为指数式下降, 符合实地采样规律, 且塑料粒子的水平分布也基本符合实测情况。2、3号测站点(厦门西海屿处)模拟值偏大的可能原因是其所处位置周围塑料的污染源点源较多, 且其所处位置也是潮流的主要通道, 接收了来自九龙江河口和同安湾的塑料垃圾。

图 6给出了厦门湾区域在模拟期间悬浮塑料粒子的最终分布特征, 从图中的结果可以看出, 厦门湾悬浮塑料粒子的分布范围介于0~2 200 particles/m², 且具有明显的空间差异性, 这也符合图 5b的采样规律。研究发现, 悬浮塑料颗粒在整个计算区域内均有分布, 并由西北向东南扩散, 在厦门岛沿岸和九龙江径流入海处出现高浓度区域, 且在青屿水道处也出现了高浓度带的分布现象。这是由于悬浮塑料颗粒在厦门湾内迁移扩散的规律主要与不同位置的流速和地形有关。根据定义的悬浮塑料颗粒运动过程的数学表达, 在相近的地形中, 其浓度分布与流速呈负相关的关系(假定流速方向相同): 在流速大的地方悬浮塑料颗粒迁移运动剧烈, 其重悬浮的速度大, 该位置的悬浮塑料颗粒浓度就会较低, 而相对流速小的位置就会呈现相反的规律。纵观整个厦门湾内漂浮塑料颗粒的分布, 发现其与各位置处的地形和距离污染源的远近也有关系: 以九龙江河口区与厦门岛西海域的交界处为例, 该位置周围有三个污染源, 所以该位置虽然流速高, 但悬浮塑料颗粒的浓度也高, 在一些内陷的豁口处也会存在高浓度微塑料聚集区。

由以上验证及分析可知, 本文所建模型可以较好地反映所关心海域(即厦门湾)的水动力场特征, 包括潮流、潮位等, 并能够较好的预测漂浮和悬浮塑料垃圾在风、潮流等作用下的输移轨迹和空间分布规律, 可用于下文中拓展研究各种类型塑料颗粒在不同风场、流场下的物理迁移规律。

2 结果与讨论

本章节利用上文中验证较好的数值模型, 选取三种不同物理性质的海洋塑料颗粒(即形状、大小及沉降速度等不同)研究其在不同海洋动力(如不同风况)作用下的运动输移特性。值得一提的是, 针对漂浮海洋塑料垃圾, 我们仍假定其始终漂浮在海洋表面, 主要关心其在风、流耦合作用下的水平迁移。而针对悬浮塑料垃圾, 则认为其在水体中同时经历了水平迁移扩散、垂直输运、沉降等复杂物理过程。严格来讲, 这两种塑料垃圾及其在水体中的行为方式并不是完全独立的, 譬如某些漂浮塑料垃圾的密度小于海水密度, 会在浮力和波浪破碎等造成的垂直湍流等的共同作用下向下沉降, 或悬浮在水体中, 或沉积于海床。根据实际采样情况, 通过编写经验公式对悬浮塑料粒子浓度、沉降通量等进行修正以考虑这些耦合作用, 是本文下一步需要深入研究的问题。

2.1 风况对漂浮海洋塑料垃圾水平飘移的影响

本研究中厦门湾的海洋风场数据采用中国海洋信息网的中国台站观测数据。整理风场数据得到, 南北方向是2 m/s及以上风速有效方向的向海区间, 且在风速为4.8 m/s以上时, 东南向的持续时间最长。而东北向是厦门的常风向, 东北向的风大约占到该月风向总时间的1/3, 最大风速约为10 m/s。最不利风向为东南风, 因为东南的向岸风会将外海的塑料垃圾吹向陆地, 造成厦门湾湾内塑料垃圾的增长。

根据以上分析, 本节在每个污染物点源释放漂浮塑料颗粒(点源具体位置见附录C), 释放速度为每时间步长1 000个, 通过改变不同的风参数, 探究最大风、最不利风及常风作用下, 厦门湾塑料垃圾漂移情况。

计算时, 最大风况的风向与实测风向保持一致, 风速大小始终为10 m/s。如图 7所示, 排污口释放的漂浮塑料粒子最终在厦门湾沿岸聚集的面积有明显增加, 同一时刻悬浮塑料颗粒局部浓度也比下文考虑其他风况作用时要高。这是因为风速的增大加速了塑料颗粒水平方向的迁移运动, 更容易堆积到岸线。

最不利风况的风速与实测风速保持一致, 风向始终为东南风, 计算结果如图 8所示。研究表明最终漂浮塑料颗粒在厦门湾沿岸聚集的面积有明显增加, 特别是在最不利风与潮流作用方向相同时, 厦门湾西北沿岸会快速聚集大量塑料颗粒。

常风向的风速与实测风速保持一致, 风向始终为东北风, 此时, 涨潮和落潮时潮流的方向与风向均存在一定夹角, 塑料颗粒会在风和潮流共同作用下迁移。如图 9所示, 与实测风向对比, 常风向作用下西南区域岸线会聚集大量的塑料颗粒, 使九龙江河口塑料污染的影响更大。

2.2 不同类别悬浮塑料垃圾的迁移特性

接下来考虑具有不同形状及沉降速度的悬浮塑料颗粒的垂直扩散特性(表 1)。排污口1号设置为排放源(图 10箭头所指处, 排污口的具体位置见附录C), 模拟结果如图 10所示。可以看出, 不同粒径和形状的塑料颗粒在海洋中的运动迁移特性有明显不同, 沉降速度大的塑料颗粒会相对快速地沉降在距离颗粒释放位置附近(即悬浮颗粒的浓度相对较小), 而沉降速度小的颗粒会随着水流漂流得更远, 沉降位置距离释放点较远。此时, 本文研究的对象主要为最长边小于50 mm的较小塑料颗粒, 而对于粒径较大的大型海洋塑料垃圾, 风的直接拖曳作用较明显, 可能对输运迁移特性造成影响, 是本文下一步需要深入系统研究的课题。

2.3 海洋塑料垃圾沉积通量的估算

本文1.3节及2.2节中建立的Eco Lab模型将塑料颗粒的沉降过程简化为当水柱中塑料颗粒浓度达到或者超过临界值时, 塑料颗粒开始向海底沉降。浓度临界值根据经验或者根据实测值校准得到。该方法的优点是简单且易实现, 但未考虑近岸海域塑料沉降到海底期间的复杂物理过程。寻找水柱表层塑料颗粒丰度/浓度(数据较易在实际采样中获得)、水柱中塑料颗粒浓度、当地沉积塑料颗粒通量间的对应关系, 并依此提出合理的参数化方案, 引入到粒子追踪或者Eco Lab模型的模拟计算中, 是提高数值模拟精度的一个有效途径。因此, 本节拟根据实际采样情况及上文模拟计算得出的总体分布规律, 编写塑料颗粒沉积通量的经验公式, 并以此探究沉积塑料通量与临近地区社会经济发展程度的相关性。同时, 该研究可为下一步建立改进的塑料颗粒迁移扩散模型提供理论依据。

参考Kukulka等(2012), 建立海洋塑料浓度的垂向分布模型:

式中, z为以海面为基准竖直向上的笛卡尔坐标, N为海水中z位置处的塑料丰度, N₀为海洋表层的塑料丰度。ω为塑料的终端上升速度, 对北大西洋塑料样本的初步试验发现, 其在静水中的终端上升速度为0.005~0.035 m/s, 而对扁平塑料碎片其取值为(0.014± 0.007) m/s (Morét-Ferguson et al, 2010)。本文参考Kukulka等(2012), ω = 0.01m/s。

A₀为塑料颗粒的垂向扩散系数, 由下式计算:

式中, k为冯卡曼系数, 取为0.4。海表面风速u^*= 0.0012W₁₀, 其中W₁₀为距海面10 m处的风速, H_s为海面的特征波高。根据厦门湾常年风况和波浪情况(林毅辉等, 2009; 朱君等, 2017), 假定海面条件为微风的情况下W₁₀ = 3.6 m/s, H_s = 0.2 m。

将塑料丰度在水柱中由海底至海水表面进行积分, 得到整个水柱中塑料通量T的大小,

根据Mani等(2016)的研究, 近岸海洋中塑料垃圾的来源主要与其所在城市生产的塑料量有关, 并且大约有80%的垃圾来自于陆源输入。因此, 在不考虑研究区域所在地的垃圾治理情况下, 假设

式中α为经验系数, 本文中塑料污染源的个数、空间分布等复杂关系主要通过调整经验系数α来考虑。e(t)为当地塑料产量随时间变化公式, 根据1949~2017年厦门塑料年产量的数据, 拟合e(t)的表达式为e(t) = e^{0.042(t–1949)} (龙邹霞, 2019)。

海洋塑料会在洋流等作用下转移到深海, 但是近岸区域塑料颗粒沉入沉积物中的机理尚不明确, 近年来对塑料沉积物分层通量的研究表明, 塑料垂向通量的分布总体上是上层大于下层(Zheng et al, 2020), 但受沉积物与水体之间不稳定的动态沉积过程和人类活动的相关影响, 塑料沉积物垂向分布并无明显变化。此外, 研究表明海洋表层塑料浓度高的位置与其对应海底沉积塑料浓度也偏高(Barrett et al, 2020)。因此, 假设沉积物通量总值S随时间变化具有如下形式,

根据龙邹霞(2019)在厦门湾鸡屿处得到的沉积物中塑料颗粒通量的数据对经验公式(7)进行验证, 如图 11所示。从图中发现, 除了在20世纪70和80年代外, 厦门湾沉积物中塑料含量与经验公式吻合良好, 近乎呈线性增长。20世纪70年代初期由于厦门塑料产业的发展, 塑料制品快速增长, 塑料沉积速率同时增大。而20世纪80年代则由于马普公约的签订, 塑料污染治理初见成效, 塑料速率增长变缓。人为因素的干扰对经验公式的编写和普适性造成了困难, 需要进一步的研究。

3 结论

本文使用数值模拟的方法, 研究了厦门湾海洋塑料垃圾迁移运动特性并对沉积物中塑料通量进行了预测。在对海洋塑料垂向分类的基础上, 基于MIKE模型, 使用拉格朗日方法和欧拉方法分别描述水体中漂浮和悬浮塑料的迁移运动, 并提出了计算沉积塑料通量的经验公式。首先, 建立了针对厦门湾的水动力模型。其次, 在正确验证水动力模型的基础上, 考虑漂浮、悬浮、沉积塑料的迁移特征分别建立了相应的数学模型。最后, 考虑了海洋塑料的分布特征和迁移规律与水动力因素之间的响应关系, 并对厦门湾沉积塑料通量进行了预测, 本文的主要结论如下:

(1) 厦门湾主要潮流通道为厦门岛东侧水道, 大小金门水道和青屿水道, 且潮流通道处的潮流呈现明显的往复流特性。涨潮时, 潮流方向大多向北和西北。落潮时, 潮流方向大多向南和东南。

(2) 漂浮塑料颗粒的迁移运动受风的影响很大, 在风的影响下进行水平方向扩散。在最大风、常风向、不利风的作用下, 大量漂浮塑料颗粒向岸边聚集, 导致局部浓度剧增。

(3) 悬浮塑料颗粒的垂向分布呈指数式下降的规律, 不同理化性质的悬浮塑料颗粒, 其水平扩散规律不同。沉降速度大的塑料颗粒水平方向传输距离更近, 而随着沉降速度减小, 悬浮塑料颗粒能够传送到更远的地方。悬浮塑料颗粒的分布还与流速有关, 一般来说流速较大的位置悬浮塑料颗粒的浓度较低, 且在该地区不容易发生悬浮塑料颗粒的沉积现象。

(4) 建立了厦门湾沉积塑料通量随时间变化的经验公式, 发现沉积塑料含量随时间呈乎呈线性增长。预测值与实测值总体吻合良好, 但经验公式未考虑人为因素的干扰, 需要进一步的研究。

附录 A 基于MIKE3的厦门湾水动力模型设置

本文所建立的厦门湾水动力模型模拟区域的总海域面积约为1 074 km²。海岸线数据采用美国海洋和大气管理局(NOAA)的GSHHG数据库, 数据集的坐标为WGS-84。水深数据来自于美国地球物理中心(NGDC) 2008年发布的ETOPO1数据库, 数据集的坐标同样为WGS-84。为保证模型的精度, 厦门湾近岸区域和九龙江河口处水深和地形采用小比例尺的海图数据进行修正, 海图选取2017年中华人民海事局发布的围头港-厦门港的6517100号海图, 比例尺大小为1:75 000。

由于地表水模拟软件(SMS)中网格编辑器生成的网格和MIKE相比更加规整, 对于海岸等不规则边界的适应性更强, 可以提高模型运算速度, 且其自带的网格质量验证功能可以很好地提升网格质量。因此, 将海岸线和水深数据在SMS软件中预处理, 生成MIKE的网格输入文件, 调整质量后的网格如附图 1所示。模型水平方向使用渐变的非结构化三角网格, 在近岸和九龙江河口处进行网格加密, 共生成10 644个节点和19 975个网格, 垂向上使用σ坐标等比例划分成5层, 网格的水深值根据插入水深散点的距离进行线性插值, 最后将制作好的地形文件导入MIKE中进行前处理。

东部和南部外海开边界条件设置为潮位强迫条件, 开边界网格节点的潮位数据通过13个主要分潮(M2、S2、N2、K2、K1、O1、Q1、P1、MF、MM、M4、MS4、MN4)的调和常数计算得到。而开边界调和常数的值则取自俄勒冈州立大学(OTIS)潮汐资料中的TPXO7.2数据集。TPXO模型基于二维正压流体方程, 运用广义反演法进行实测数据同化, 并采用最小二乘法进行数据拟合。该模型同化数据包括卫星测高数据(T/P、Topex Tandem、ERS、GFO)和实测验潮数据, 模型预测范围覆盖全球海平面, 模型的结果已经广泛应用于我国的沿海地区(范文蓝等, 2018; Wang, 2019)。本文采用TMD (Tide Modal Driver)工具箱计算得到开边界网格节点每隔一小时的预报数据, 其他时间的潮位数据通过MIKE模型线性插值确定。同时考虑到九龙江河口径流会携带大量来自上游的微塑料垃圾, 根据骆智斌等(2008)对九龙江河口的水动力潮流潮汐模型研究结果, 在九龙江边界加入多年平均径流量371 m³/s作为流量控制边界, 陆地边界设置为闭边界, 即将法向流速设置为0。

风场数据则来自于国家海洋科学数据中心提供的中国台站观测数据集, 将厦门测点站每隔1 h的实测数据应用于整个模型区域, 规定风摩擦系数为0.001 255。为避免模型产生大的波动, 设置软启动间隔为7 200 s, 在此期间内风速从零增加到指定值。在计算中使用可变的时间步长间隔, 使收敛条件判断数(Courant Friedrich Levy, CFL)在所有计算节点中小于临界CFL值。规定临界CFL值为0.8, 最小时间步长为0.01 s, 最大时间步长为30 s。

模型通过设置干湿临界水深和淹没临界水深处理移动边界的问题, 每个单元根据水深值和单元表面是否淹没分成干燥、部分干燥和湿三种不同的状态, 其中干燥单元将会在计算中被移除, 部分干燥单元忽略其动量通量, 只计算质量通量。湿单元则计算动量通量和质量通量。本文规定干水深临界值为0.005 m, 湿水深临界值为0.1 m, 淹没水深临界值为0.05 m。

在近海潮波模拟中, 测站点多分布在浅水区, 水深一般都不大, 在这种情况下, 底部的摩擦效应较深水而言要重要得多, 本文通过粗糙高度设置底床糙率, 不同位置的粗糙高度取值不同, 范围为0.01~ 0.05 m。另外, 模型还考虑了涡流黏度和地球自转产生的科式力的影响。水平涡流黏度通过Smagorinsky公式指定, 规定Smagorinsky系数为0.28, 涡流黏度的最小值和最大值分别为0.000 018和10 000 000 000 m²/s。垂直涡流黏度选择k-ε模型计算, 最小值和最大值分别为0.000 018和0.4 m²/s。

为了使模型更加平滑地启动, 本文将水动力模型设置为冷启动, 将模型初始运行时刻外海开边界潮位的平均值0.62 m作为初始水位, 初始流速设置为0。

B Eco Lab生态模型过程变量数学表达式编写

本文将海洋中悬浮塑料颗粒的物理运输过程简化为以下四个过程:

(1) 塑料颗粒的大气沉降: 风会携带塑料颗粒进入大气中, 大气到海洋表面的塑料颗粒沉降十分常见; (2) 对流扩散: 通过平流和扩散在水柱中运输塑料颗粒; (3) 沉积: 密度大于海水的塑料颗粒会从水柱中沉积到底床表面; (4) 再悬浮: 当海水流速大于海底沉积塑料颗粒的再悬浮临界速度时, 塑料颗粒会重新悬浮到水柱中。

将悬浮塑料颗粒(1)、(3)、(4)的三个基本运输过程在MIKE生物仿真模块中通过开放式的数学公式定义, 并耦合对流扩散模型则可对悬浮塑料颗粒的运动过程进行较为完整的模拟。

在ECO Lab创建悬浮塑料颗粒迁移模板, 分别定义状态变量(state variables, 描述生物个体包括微塑料颗粒对水环境参数做出的响应的物理量, 如悬浮微塑料的浓度等)、常数(constants, 不随时间变化的物理量, 如悬浮微塑料的沉降速度等)、过程(processes, 与微塑料迁移运动相关的物理过程, 如前文所述的大气沉降、对流扩散等)、作用力(forcing, 外力作用, 如水动力模块输出的流场、风场信息)、公式表达(expression)、附加输出(derived outputs)等部分。其中, 塑料颗粒模板定义的状态变量为悬浮塑料颗粒的浓度值和沉积塑料颗粒的面密度。记悬浮塑料颗粒的浓度值为MPSS, 沉积物塑料颗粒的面密度为MPSED, 则悬浮塑料颗粒浓度变化可由以下常微分方程描述:

其中, MPSS为悬浮塑料颗粒的浓度, 单位为mg/L; dt为时间步长, 单位为d; prss、sessv、ressv分别为大气沉降过程塑料颗粒输入浓度变化率, 悬浮塑料颗粒的沉积率和再悬浮率, 单位均为g/(m³·d)。

式(B-1)等号右边分别为塑料颗粒的大气沉降过程、悬浮塑料颗粒的沉积过程和沉积物中塑料颗粒的再悬浮过程的总和, 这三个过程具体的公式表达如下:

式中, papro为大气沉降到海水表面的塑料颗粒沉降速率, 参考Eco Lab关于污染物扩散模型的设置, 假定其主要与空气中塑料颗粒的浓度有关, 其与风速风向的相关性则需要开展进一步研究; sessv为单位体积沉降速率; dz为垂直方向上划分的层高; vsm为悬浮塑料颗粒的沉降速度, 与塑料颗粒的相对粒径, 形状大小(球状或碎片等)等有关; 本文将沉降过程简化为当悬浮物浓度达到一定值时, 塑料粒子开始向海底沉降, 其浓度临界值可以根据经验或者实测值决定。

ressa为单位面积塑料颗粒悬浮率, 可由下式逻辑表达式求得,

该逻辑表达式表示了沉积物中塑料颗粒再悬浮的条件既要满足水柱中悬浮塑料颗粒达到饱和浓度, 同时也要满足流速大于再悬浮的临界流速。其中, resrat为最大再悬浮率; bssmx表示水中悬浮是否未饱和, 其中1为饱和, 0为未饱和; bresu表示是否悬浮, 1为悬浮, 0为非悬浮。塑料颗粒的重悬浮率与水流大小及沉积塑料颗粒的沉降速度有关, 相同材质的塑料, 粒径大的更容易悬浮起来, 因为粒径大的塑料粒子受到的浮力更大, 更容易在水流的作用下产生扰动。注: 本文将再悬浮过程进行简化, 只考虑了交换过程中再悬浮量的大小, 而未考虑水流的流向等的影响。

C 塑料垃圾的来源及入海点

对于沿岸城市向海洋排放塑料垃圾的研究, 本文综合分析后采用陆源输入、海滨旅游、排污口排污和船舶垃圾四个参数作为海洋微塑料垃圾的来源, 将厦门湾主要的排污口处、人为活动区域密集的厦门岛西海域和九龙江径流处设置为微塑料排放的点源。同时在厦门北航道和厦门岛近岸分别设置移动点源和面源以考虑滨海旅游和船舶垃圾入海的效应。其中各点源的位置如附图 2所示。