附录

附录

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

本文所建立的厦门湾水动力模型模拟区域的总海域面积约为1 074 km2。海岸线数据采用美国海洋和大气管理局(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中进行前处理。

附图1 厦门湾海域计算网格

Supplementary Fig.1 The computational grid of the Xiamen Bay

东部和南部外海开边界条件设置为潮位强迫条件, 开边界网格节点的潮位数据通过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 m3/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 m2/s。垂直涡流黏度选择-模型计算, 最小值和最大值分别为0.000 018和0.4m2/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, 则悬浮塑料颗粒浓度变化可由以下常微分方程描述:

dMPSS/d= prss–sessv+ressv, (B-1)

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

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

prss = papro/d

sessv = vsm×MPSS×86 400 (B-2)

ressv = ressa/d

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

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

ressa = If bssmx×bresu==1

Then min(resrat, (Max(0, XSED–1))/d×86 400)(B-3)

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

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

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

附图2 厦门湾海洋微塑料垃圾入海点源

Supplementary Fig.2 The sources of microplastic around the Xiamen Bay