海洋牧场是指在一定海域内, 采用规模化渔业设施和系统化管理体制, 利用自然的海洋生态环境, 将经济海洋生物聚集起来, 对鱼、虾、贝、藻等海洋资源进行有计划和有目的的海上放养^[1-4]。网箱养殖是海洋牧场规模化养殖的一种重要手段, 通过在自然海区中设置网箱, 依靠自动流水和人工投饵的生产方式进行高密度海水产养殖^[5]。网箱养殖机动、灵活、简便、高产、水域适应性广^[5-6], 但在半封闭的内湾中网箱养殖分布密集, 形成大范围内超负荷的养殖状况, 可能会造成周围海域海水水质的恶化^[7-10]。研究表明, 网箱养殖对水质的影响主要表现为营养负荷的增加^[8-10], 养殖过程中产生的残饵、粪便含有氮、磷等营养物质 ^[11-14]。超出近岸海域环境承载力和自净能力的氮、磷等营养物质的排放会引发水体的富营养化及其他形式的海洋生态环境破坏和异常^[13-15]。

本文通过研究“耕海1号”海洋牧场网箱养殖所在海域的环境条件, 在水动力环境潮流场数值模拟的基础上, 构建了典型海洋牧场网箱养殖污染预警系统, 计算了海洋牧场所在海域受纳水体对养殖生产排泄物的响应关系, 从而评估海洋牧场网箱养殖活动对周边海域的影响。

1 研究对象及控制方程

“耕海1号”位于烟台市四十里湾海域(121°27′ 52.577″~121°28′01.265″E, 37°31′02.859″~37°31′17.662″N), 位置如图 1所示(图件底图为卫星遥感图, 来源于山东省地理信息公共服务平台)。本文研究对象是其大型智能化网箱养殖区, 该区域计划养殖鲷鱼、马面鲀及其他鱼类, 养殖种类及相关参数见表 1。项目于2020年开始运营, 截至目前, 斑石鲷和真鲷年均产量约269.5 t, 马面鲀约24.7 t, 黑鲪、黄带拟鲹、半滑舌鳎等其他鱼类约101.5 t。养殖过程中均采用人工颗粒饲料, 饵料系数约1∶1.2。

1.1 污染物扩散模型

采用平面二维数值模型研究工程海域的潮流场运动及海域污染物扩散影响, 模型采用非结构三角网格剖分计算域; 采用标准Galerkin有限元法进行水平空间离散, 在时间上, 采用显式迎风差分格式离散动量方程与输运方程。

1) 模型控制方程

质量守恒方程:

$ \frac{{\partial \zeta }}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {hu} \right) + \frac{\partial }{{\partial y}}\left( {hv} \right) = 0 . $

(1)

动量方程:

$ \begin{aligned} \frac{\partial u}{\partial t} & +u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}-\frac{\partial}{\partial x}\left(\varepsilon_x \frac{\partial u}{\partial x}\right)-\frac{\partial}{\partial y}\left(\varepsilon_x \frac{\partial u}{\partial y}\right) \\ & -f v+\frac{g u \sqrt{u^2+v^2}}{c_z^2 H}=-g \frac{\partial \zeta}{\partial x} \end{aligned} $

(2)

$ \begin{aligned} \frac{\partial v}{\partial t} & +u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}-\frac{\partial}{\partial x}\left(\varepsilon_x \frac{\partial u}{\partial x}\right)-\frac{\partial}{\partial y}\left(\varepsilon_y \frac{\partial u}{\partial y}\right) \\ & +f u+\frac{g v \sqrt{u^2+v^2}}{c_z^2 H}=-g \frac{\partial \zeta}{\partial y} \end{aligned} $

(3)

式中, $\zeta $为水位; h为静水深; H为总水深, H=h+$ \zeta $; u为x向垂向平均流速; v为y方向垂向平均流速; g为重力加速度; f为科氏力参数(f = 2ω sin $\varphi $, $\varphi $为计算海域所处地理纬度); ${c_z}$为谢才系数, ${c_z} = \frac{1}{n}{H^{\frac{1}{6}}}$, n为曼宁系数; ${\varepsilon _x}$、$ {\varepsilon _y} $为x、y方向水平涡动黏滞系数。

2) 定解条件

初始条件:

$ \left\{\begin{array}{c} \left.\zeta(x, y, t)\right|_{t=t_0}=\zeta\left(x, y, t_0\right)=0 \\ \left.u(x, y, t)\right|_{t=t_0}=\left.v(x, y, t)\right|_{t=t_0}=0 \end{array} .\right. $

(4)

边界条件:

开边界: 全球模型调和求得开边界的M₂、S₂、K₁、O₁、M₄和MS₄ 6个分潮调和常数值输入计算。

$ \zeta = \mathop \sum \limits_{i = 1}^N \left\{ {{f_i}{H_i}\cos \left[ {{\sigma _i}t + \left( {{V_{0i}} + {V_i}} \right) - {G_i}} \right]} \right\} , $

(5)

式中, f_i、σ_i是第i个分潮(这里共取6个分潮: M₂、S₂、K₁、O₁、M₄和MS₄)的交点因子和角速度; H_i和G_i是调和常数, 分别为分潮的振幅和迟角; V_0i+V_i是分潮的幅角。

1.2 预测模型及计算方法

描述该海域的水体运动和污染物的扩散选用二维浅水方程组和物质平衡方程。本次模拟选用MIKE21中transport模块进行计算^[16-19]。污染物在海水中的迁移扩散过程, 由二维对流、扩散方程表示:

$ \begin{aligned} \frac{\partial C_i}{\partial t}+U \frac{\partial C_i}{\partial x}+V \frac{\partial C_i}{\partial y}= & \frac{\partial}{\partial x}\left(E_x \frac{\partial C_i}{\partial x}\right) \\ & +\frac{\partial}{\partial y}\left(E_x \frac{\partial C_i}{\partial y}\right)+K_i C_i+S_i \end{aligned} $

(6)

式中, C_i为污染物浓度; u、v为x、y方向上的流速分量; E_x、E_y为x、y向上的扩散系数; K_i为污染物降解系数; S_i为污染物源汇项。

方程包括3大项: 物理输移扩散项、生化项及源汇项。定解条件如下:

$ A: 在闭边界上没有物质通量, 即\frac{{\partial C}}{{\partial t}} = 0. $

(7)

$ B: 在外海开边界上流出满足\frac{{\partial C}}{{\partial t}} + {V_n}\frac{{\partial C}}{{\partial n}} = 0. $

(8)

$ C: 在外海开边界上流入满足C(x, y, t)=0 .$

(9)

1.3 水质检测方法

收集的水样经0.45 μm滤膜过滤后进行水质检测分析。检测项目包括溶解氧、化学需氧量(chemical oxygen demand, COD)、无机氮、活性磷酸盐, 分析方法按照《海洋监测规范第4部分: 海水分析》(GB17378.4—2007)^[20]和《海洋调查规范第4部分: 海水化学要素调查》(GB/T 12763.4— 2007)^[21]执行。

2 水动力模型设置及验证 2.1 模型设置 2.1.1 计算域设置

本项目所建立的海域数学模型计算域范围见图 2(a), 计算域坐标范围为37°04′14.22″~ 40°58′08.25″N, 117°29′33.27″~122°41′36.62″E。

模拟采用三角网格, 用动边界的方法对干、湿网格进行处理。整个模拟区域内由14 639个节点和23 426个三角单元组成, 最小空间步长约为20 m。为了能清楚了解本工程附近海域的潮流状况, 将本工程附近海域进行局部多次加密, 加密区域见图 2(b)。

2.1.2 水深和岸界

水深和岸界根据中国人民解放军海军航海保证部制作的11500号、11700、11800号、11300号、11900号、12100号海图以及工程周边实测水深和岸线确定。

2.1.3 模型开边界条件

全球潮汐模型调和常数插值计算得到开边界的M₂、S₂、K₁、O₁、M₄和MS₄ 6个分潮调和常数值输入计算。

$ \zeta = \mathop \sum \limits_{i = 1}^N \left\{ {{f_i}{H_i}\cos \left[ {{\sigma _i}t + \left( {{V_{0i}} + {V_i}} \right) - {G_i}} \right]} \right\}, $

(10)

式中, f_i、σ_i是第i个分潮(这里共取6个分潮: M₂、S₂、K₁、O₁、M₄和MS₄)的交点因子和角速度; H_i和G_i是调和常数, 分别为分潮的振幅和迟角; V_0i+V_i是分潮的幅角。

闭边界: 以大海域和工程周边岸线作为闭边界。

2.1.4 计算时间步长和底床糙率

模型计算时间步长根据CFL条件进行动态调整, 确保模型计算稳定进行, 最小时间步长1 s。底床糙率通过曼宁系数进行控制, 曼尼系数n取32~ 55 m^1/3/s。

2.1.5 水平涡动黏滞系数

采用考虑亚尺度网格效应的Smagorinsky公式计算水平涡黏系数^[22], 表达式如下:

$ A = c_s^2{i^2}\sqrt {2{S_{ij}}{S_{ij}}} , $

(11)

式中, c_s为常数。

2.2 模型验证

模拟区内潮位、潮流验证点分布见图 1。

2.2.1 潮位验证

利用烟台港潮位站历史观测资料经调和分析后, 选用M₂、S₂、K₁、O₁、M₄和MS₄ 6个分潮的调和常数预报出大潮期的潮位与计算结果进行验证。模拟区内潮位验证点见图 2和表 2, 潮位验证曲线见图 3。

2.2.2 潮流验证

计算海域现有4个实测流观测资料(表 2)。图 4分别是L01—L03 (2015-05-21—2015-05-22)和C01 (2019-05-22—2019-05-23)流速、流向模拟值与实测值比较图。由图中可以看出, 流场的数值模拟结果流速和流向, 都与实测资料变化基本一致, 流速相对误差基本在20%以内, 吻合较好。

以上潮位和潮流验证结果表明, 相应验证点上潮位和潮流模拟结果与实测潮位和潮流资料基本吻合, 能够较好地反映工程周边海域潮流状况。

2.3 潮流响应特征

海流数据来源于中国海洋大学姜胜辉团队于2019年8月3日9: 00至8月4日10: 00大潮期间在“耕海1号”附近海域进行的海流同步观测资料^[23], 海流观测调查资料显示: 涨潮流流向主要集中出现在偏东南向, 落潮流流向主要集中出现在西北向。各层潮流类型判别数均小于2, 其潮流性质主要为半日潮流, 潮流运动形式以往复流为主。大潮期表层平均最大流速为16.6 cm/s, 流向为326.8°; 中层平均最大流速为15.4 cm/s, 流向为341.6°; 底层平均最大流速为11.5 cm/s, 流向为346.0°。

涨潮时, 湾内平均流速在17 cm/s左右, 整体流向自东北向西南方向, 受崆峒岛的影响, 工程海域流向为自北向南, 流速在10~15 cm/s, 如图 5(a)所示; 落急时, 湾内平均流速在22 cm/s左右, 流向沿西北方向向外流出湾口, 工程附近海域流向为自东南流向西北, 流速在15~20 cm/s, 如图 5(b)所示。大潮涨急、落急时工程建设前后流速对比如图 6所示, 工程建设后工程附近大部分区域表现为流速减小, 最大减小量为8 cm/s, 工程两侧小部分区域表现为流速增大, 最大增加量为6 cm/s。受桩基阻隔影响, 项目中间围栏养殖区域流速也相对减小, 减小程度约为3 cm/s。

3 养殖污染物预测分析与实际环境影响评价

网箱养殖完全依赖人工投饲, 投喂的饵料大部分被鱼类摄食, 还有一部分则散失在水体中, 导致水体中的有机物增多, 氮、磷浓度增加。另外, 鱼类的排泄物也会对水体造成影响。海洋牧场养殖污染物入海后要经过一系列物理和化学变化, 不同海域水体自净能力和污染物容纳量各不相同, 主要取决于流体动力过程和生物化学过程^{[19, 24-25]}。海域的物理自净能力有两种机制, 一是潮流的动力搬迁作用, 二是水体的稀释扩散作用^[25-27]。描述该海域的水体运动和污染物的扩散选用二维浅水方程组和物质平衡方程, 模拟选用MIKE21中transport模块进行计算。模拟计算中所用的海流数据来源于中国海洋大学姜胜辉团队于2019年8月3日9: 00至8月4日10: 00大潮期间在“耕海1号”附近海域进行的海流同步观测资料^[23], 观测站位距离本项目约6 km。

3.1 养殖污染物环境浓度预测分析 3.1.1 排污点位置

养殖排污范围为图 7(a)蓝色区域, 本次在网箱养殖外沿选取的排污代表点进行计算, 排污点位置见图 7(a)。

3.1.2 海洋牧场产生的污染源强分析

网箱养殖排污系数参考《第一次全国污染源普查水产养殖业污染源产排污系数手册》^[28]。根据排污系数及各鱼类设计年产量计算得出: 鲷鱼(包括斑石鲷和真鲷)污染物的年排放量为COD 21.70 t, 总氮21.61 t, 总磷3.62 t; 马面鲀污染物的年排放量为COD 3.80 t, 总氮3.78 t, 总磷0.29 t; 其他鱼类(半滑舌鳎、黑鲪和黄带拟鲹)污染物的年排放量为COD 11.39 t, 总氮11.34 t, 总磷0.88 t。项目网箱养殖污染物的年排放总量为: COD 36.73 t, 总氮36.89 t, 总磷4.80 t。具体各类养殖生物的排污源强计算结果见表 3。考虑最大排污情况, 网箱养殖外沿每个代表点的排放源强为: COD 1.42 g/s, 总氮1.42 g/s, 总磷0.19 g/s。

3.1.3 海洋牧场周围海域分析

污染物背景浓度参考中国海洋大学姜胜辉团队于2019年8月3日9: 00至8月4日10: 00大潮期间在“耕海1号”附近海域进行的海流同步观测资料^[23], 这里选取调查资料中距离本项目最近的27号站位的实测值作为海水水质的背景值, 各参数背景值见表 4。

3.1.4 养殖污染物环境浓度预测分析

养殖活动主要产生的污染物为COD、总氮、总磷, 本次按照扩散方程对COD、无机氮、活性磷酸盐的浓度进行预测。

根据计算, 叠加背景值后, 全潮CODr、无机氮、活性磷酸盐的浓度包络线范围见图 7(b)—(d), 各级预测浓度的扩散范围统计见表 5。网箱附近COD最大预测浓度约1.033 mg/L, 小于2 mg/L; 无机氮最大预测浓度约0.095 mg/L, 小于0.2 mg/L; 活性磷酸盐最大预测浓度约0.009 1 mg/L, 小于0.02 mg/L。COD、无机氮、活性磷酸盐的预测浓度水平均满足《中华人民共和国海水水质标准(GB3097—1997)》中一类水质标准。

3.2 海洋牧场运行后实际产生的环境影响

根据前文描述, 网箱养殖年产总量395.67 t, 达到了设计年产量的80%左右。通过对项目周围海域7个站位的海水水质进行检测分析(分析结果见表 6), 发现自2020至2023年, 项目周围海域海水COD质量浓度为0.312~1.37 mg/L, 磷酸盐质量浓度最大值达到0.004 73 mg/L, 亚硝酸盐氮最大值达到0.002 75 mg/L, 硝酸盐氮最大值达到0.127 mg/L, 氨氮在0.001 47~0.027 5 mg/L。COD值小于2 mg/L, 无机氮小于0.2 mg/L, 磷酸盐小于0.015 mg/L, 均满足《中华人民共和国海水水质标准(GB3097—1997)》中一类水质标准。项目运行后, 水质检测结果与预测结果一致, 这表明选用MIKE21中transport模块进行污染物浓度预测的合理性和科学性, 同时也进一步表明“耕海1号”海洋牧场对周围海域水环境的影响在可控范围内。相较于烟台北部海域和胶州湾的最大流速, 四十里湾的潮流流速缓慢^[29], 污染物在水体中的扩散迁移范围固然有限, 但缓慢的潮流有可能会导致养殖污染物沉降在周边区域的沉积物中, 因此后续研究中网箱养殖对养殖系统及周边海域沉积环境的潜在影响同样也应当受到关注。

4 结论

构建了“耕海1号”海洋牧场网箱养殖污染预警系统, 通过数值模拟智能评估了“耕海1号”海洋牧场对所在海域海水水质影响, 计算了海洋牧场所在海域受纳水体对养殖活动的响应关系, 结果如下。

(1)“耕海1号”海洋牧场养殖过程中会对周边海洋环境产生一定的影响, 主要表现为项目运营期养殖鱼类排泄物会导致邻近区域COD、无机氮、活性磷酸盐等污染物浓度有所升高, 但增量有限, 符合《海水水质标准(GB3097—1997)》中一类水质标准要求。

(2)“耕海1号”海洋牧场养殖过程中COD、无机氮、活性磷酸盐等污染物影响范围很小, 扩散最远仅180 m。

(3) 海洋牧场运行4 a后, 通过对周围海域的海水水质进行检测分析, 发现实际水质与预测结果一致, 海洋牧场对海水环境的影响在可控范围内。

致谢: 文中涉及的生态环境调查数据、海流数据等由中国海洋大学姜胜辉团队提供, 在此一并致谢。