象山港位于浙江省北部沿海, 由象山港峡湾、牛鼻水道和佛渡水道三大部分构成(图 1), 海域面积为563 km², 平均水深为10 m, 是一个纵长约为70 km的典型的狭长型半封闭海湾^[1]。自20世纪90年代初, 该海域大力发展临港工农业, 推进网箱养殖等多种渔业生产方式, 区域海洋经济水平迅猛提高, 现已发展成为重要的海水增养殖基地和渔业资源基地^[2-4]。随着临港工农业、海洋渔业等的快速发展, 大量溶解无机氮(Dissolved Inorganic Nitrogen, DIN)经河流、水闸等入海口汇入港内^[5-7], 同时由于象山港是一个半封闭海湾, 港内(尤其是狭湾处)水交换能力较弱, 水交换时间较长^[1], 使得DIN难于向外海输运并滞留湾内, 致使象山港海域富营养化状态日益严峻^{[4, 8-10]}。Zhu等^[11]指出DIN是象山港海域主要的污染物且自1982年至2011年DIN含量呈现逐渐增加的趋势。根据《2017年象山港海洋环境公报》数据, 象山港DIN含量在0.16~1.19 mg/L, 平均含量为0.65 mg/L, 劣于第四类海水水质标准, 赤潮等有害藻华现象频发, 对海洋生态、经济、社会效益产生了严重影响^[12-14], 因而对DIN实行总量控制是很有必要的。

计算、优化分配DIN的环境容量是进行DIN总量控制的基础。环境容量是指目标海域在满足某一海水水质标准的要求下, 一定时间范围内污染物的最大允许排海通量^[15-17]。针对象山港DIN环境容量, 前人多基于浓度场分担率法进行计算^[16-17], 该方法虽简便易行、计算量小, 但其无法对各入海口的排海通量进行优化分配。排海通量最优化法针对该方法的不足, 通过优化分配各入海口的排海通量, 使各水质控制点在满足一定海水水质标准的前提下, 其允许排海通量总和达到最大进而得到目标海域的环境容量。目前该方法已应用于渤海、胶州湾、莱州湾等海域环境容量的计算^{[15, 18-20]}。本文将建立象山港三维水质模型, 对该海域DIN浓度进行数值模拟, 基于排海通量最优化法优化分配各入海口的排海通量并计算该海域DIN的环境容量。

1 研究方法 1.1 模型配置

本文通过耦合ROMS(Regional Ocean Model System)及NEMURO(North Pacific Ecosystem Model for Understanding Regional Oceanography)模式, 建立了象山港三维水质模型。其中ROMS模式基于三维非线性斜压原始方程开发, 在近海水动力环境模拟中是一种常用模型; NEMURO模型是一个低营养级水质模型, 最初被应用于北太平洋中部, 后被应用于包括东中国海在内的多个海区^[21-23]。模型模拟海区经纬度范围为121.41°E—122.13°E, 29.25°N—29.88°N, 水平方向采用正交经纬网格, 网格数为360×360, 分辨率为(1/500)°, 垂直方向采用地形跟随的S坐标, 共分10层。地形数据由分辨率为600 m的海图数据插值得到^[24]。模型综合考虑了风场、潮流、密度流及陆源污染物入海口的共同作用, 其中风场、海表热通量、淡水通量等数据来自COADS(Comprehensive Ocean Atmosphere Data Set)数据集; 潮汐强迫考虑了8个主要分潮M₂, S₂, K₁, O₁, N₂, K₂, P₁, Q₁的影响; 开边界处水位及流速数据来自SODA(Simple Ocean Data Assimilation)数据集; 陆源污染物入海口考虑了凫溪、颜公河、郭巨大碶等10个入海口(图 1)的影响, 各入海口径流量及DIN数据由宁波市环境保护科学研究设计院提供。凫溪(S6)为港内年径流量最大的河流, 图 2给出了其径流量及DIN排海浓度的季节变化。同时为简化模型, 本文将DIN视作保守物质, 水质模型中部分参数借鉴王玉衡^[25]及Liang等^[26]的工作。模型对象山港DIN时空分布进行了气候态模拟且运转一年后稳定, 下文工作主要基于第二年模拟结果展开。

1.2 排海通量最优化法

排海通量最优化法的目标是求解一类最优化问题, 其要求各污染源的允许排海浓度之和达到极大值, 所以目标函数为:

$ \max \sum\limits_{j = 1}^n {{R_j}, } $

(1)

同时各水质控制点的污染物浓度符合目标海域海水水质标准要求, 即:

$ C_i^0 \le {C_i}, $

(2)

而水质控制点污染物浓度大小取决于其与各入海口污染物排海浓度之间的水质响应函数, 就保守物质而言, 水质响应函数一般呈线性关系^[15], 所以其约束条件为:

$ \left( {\begin{array}{*{20}{c}} {{a_{11}}}\\ {{a_{21}}}\\ \vdots \\ {{a_{i1}}} \end{array}\begin{array}{*{20}{c}} {{a_{12}}}\\ {{a_{22}}}\\ \vdots \\ {{a_{i2}}} \end{array}\begin{array}{*{20}{c}} \cdots \\ \cdots \\ \vdots \\ \cdots \end{array}\begin{array}{*{20}{c}} {{a_{1j}}}\\ {{a_{2j}}}\\ \vdots \\ {{a_{ij}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{R_1}}\\ {{R_2}}\\ \vdots \\ {{R_j}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{\beta _1}}\\ {{\beta _2}}\\ \vdots \\ {{\beta _i}} \end{array}} \right) \le \left( {\begin{array}{*{20}{c}} {{C_1}}\\ {{C_2}}\\ \vdots \\ {{C_i}} \end{array}} \right), $

(3)

此外, 可对各入海口的污染物最低排海浓度做出约束:

$ {R_j} \ge {R_{\min }} \ge 0. $

(4)

其中i为水质控制点个数, j为入海口个数, R_j为第j个入海口的污染物排海浓度, R_min为各入海口的污染物最低排海浓度,$C_i^0$为各水质控制点处的污染物浓度, C_i为各水质控制点处由目标海域海水水质标准所要求的污染物浓度, α_ij为第i个水质控制点处污染物浓度与第j个入海口的污染物排海浓度之间水质响应函数的系数, β_i第i个水质控制点处污染物浓度与各入海口的污染物排海浓度之间水质响应函数的截距之和。求解上述目标函数(式1)与约束条件(式2—4)所界定的最优化问题即可得到目标海域污染物的环境容量^[15]。

2 模型验证

首先, 利用湾顶和湾口处的水位周日观测数据进行潮位验证(图 3a, b), 其中内湾B2站(121.56°E, 29.50°N)水深11.5 m, 观测时间为2012年12月19日至20日; 湾口D2站(122.07°E, 29.71°N)水深3.0 m, 观测时间为2012年4月9日至11日^[27]。两站相对误差分别为15.08%、16.51%, 变化趋势基本一致; 其次, 参考Jiang等^[2]于2010年1、4、7、11月观测所得表层水温数据, 模型能够重现象山港水温的季节变化(图 3c); 第三, 选取宁波市环境保护科学研究设计院提供的DIN浓度现场观测数据与模拟DIN浓度进行比较(图 3d), 两者相对误差(表 1)最大为–19.15%, 最小为–1.82%, 且DIN空间分布均呈现自湾顶至湾口降低的变化趋势。最后, 对水位及表层u、v数据进行潮汐、潮流调和分析(图 4), 象山港潮汐类型属于非正规半日潮, M₂分潮占主导地位, 其振幅与迟角自湾口向湾顶逐渐增大。潮流呈往复流特征, M₂分潮流椭圆长轴走向基本与岸线平行。潮流流速值最大出现在湾口, 往湾顶流速逐渐减小, 与韩松林等^[28]刻画的潮汐潮流特征较为一致。总体来看, 该模型模拟结果与实测数据基本吻合, 能够反映该海域水动力场及DIN浓度的变化规律。

3 计算结果 3.1 水质响应函数的建立

参考虞兰兰等^[20]选取水质控制点的原则, 设置各入海口(S1, S2, …, S10)口外400 m处作为水质控制点(P1, P2, …, P10)。在保持各入海口DIN排海浓度不变的基础上, 利用水质模型依次计算其中一个入海口的DIN排海浓度增大为原来的2、4、6、8、10倍时水质控制点处DIN浓度的变化, 经函数拟合即可建立起该水质控制点处DIN浓度与各入海口的DIN排海浓度之间的水质响应函数。各水质响应函数均呈现明显的线性关系, 主要是由于在模拟过程中, DIN被视为保守物质^{[15, 19-20]}。受篇幅限制, 本文只列出凫溪(S6)口外的水质控制点P6与郭巨大碶(S1)、联胜新碶(S2)、狮子口闸(S3)、下陈江(S4)之间的水质响应函数(图 5)。

3.2 入海口DIN排海通量的线性规划

结合各水质控制点水质要求及其与各入海口之间的水质响应函数, 目标函数(式1)与约束条件(式2—4)所确定的最优化问题有解的充分条件为各入海口DIN最低排海浓度R_min介于0~0.17 mg/L。基于此将R_min自0开始, 以0.01 mg/L递增至0.17 mg/L设置18组线性规划实验, 实验发现不同R_min的取值对应需要削减的入海口并不相同, 由此可将上述18组线性规划实验结果归纳为a(0≤R_min≤0.10 mg/L)、b(0.10 < R_min≤0.12 mg/L)、c(0.12 < R_min≤0.14 mg/L)、d(0.14 < R_min≤0.17 mg/L)四类情境:情境a中, 即各入海口DIN最低排海浓度R_min介于0~0.10 mg/L时, 仅凫溪(S6)与颜公河(S7)需要削减DIN排海通量, 该条件下, DIN环境容量最大, 约为274.17 t/a; 情境b中, 除了上述两个入海口, 西周排污口(S8)也需要削减, 对应的DIN环境容量为249.02 t/a; 情境c中需要削减的入海口更多, 相应的DIN环境容量减小为229.41 t/a; 情境d中, 除郭巨大碶外, 各入海口均需要不同程度地削减DIN排海通量, 此时DIN环境容量最小, 约为127.16 t/a, 其与张永胜^[17]基于浓度场分担率法计算的DIN环境容量结果143.21 t/a较为接近。综上, 象山港DIN环境容量为127.16~274.17 t/a, 且其随着R_min的增大(由情境a至d)逐渐减小(表 2)。

结合表 2可以发现随着R_min的逐渐增大, 各入海口DIN排海通量需要依次削减, 因而可以按削减优先程度对各入海口做如下分级: Ⅰ级:凫溪(S6)、颜公河(S7); Ⅱ级:西周排污口(S8); Ⅲ级:下陈江(S4)、朱娘桥河(S5); Ⅳ级:联胜新碶(S2)、狮子口闸(S3)、雅林溪(S9)、贤庠河(S10)。水体交换时间指某一海域内保守物质浓度稀释至初始浓度值的一半时所需要的时间, 是描述海湾物质输运的重要指标之一^[29]。图 6为各入海口的水体半交换时间及DIN年排海通量分布情况, 由图 6可知Ⅰ级入海口所处海域水体半交换时间最长且DIN年排海通量最大, 表明该海域水体交换能力较弱, DIN难于向外海输运并滞留于该海域, 因而Ⅰ级入海口的DIN排海通量削减优先程度最高; 虽然Ⅲ级入海口(下陈江与朱娘桥河)所处海域水体半交换时间长于Ⅱ级(西周排污口), 但由于西周排污口的DIN年排海通量比III级入海口高了近2.3倍, 因此西周排污口削减优先程度高于两者; Ⅳ级入海口靠近外海, 水体半交换时间较短、水体交换能力较强, 其削减优先程度最低。

4 结论

本文基于ROMS及NEMURO数值模型, 综合考虑凫溪、颜公河、郭巨大碶等代表性入海口的影响, 建立了象山港三维水质模型, 经检验, 模拟结果与观测结果基本一致, 可用于象山港海域DIN环境容量的计算。在各水质控制点满足海域功能区划要求的前提下, 通过调整R_min的取值设置了a、b、c、d四类情景并计算了对应情境下的DIN环境容量。结果表明:

(1) 不同情境下, 象山港海域DIN环境容量为127.16~274.17 t/a, 且随着R_min的增大, DIN环境容量逐渐减小。

(2) 随着R_min的逐渐增大, 各入海口DIN排海通量依次需要削减, 按削减优先程度可将各入海口划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ四级, 即Ⅰ级凫溪、颜公河, Ⅱ级西周排污口, Ⅲ级下陈江、朱娘桥河, Ⅳ级联胜新碶、狮子口闸、雅林溪、贤庠河。削减优先程度与各入海口所处海域的水体半交换时间及其DIN年排海通量有关:水体半交换时间越长、DIN年排海通量越大, 削减优先程度越高。