渤海氮磷营养盐和叶绿素浓度时空分布数值模拟与富营养化评估*

赵 宇 李克强① 孙 珊 陈 衎 谭光深 张 娟 王修林

(1. 中国海洋大学海洋化学理论与工程技术教育部重点实验室 中国海洋大学化学化工学院 山东青岛 266100; 2. 山东省海洋资源与环境研究院 山东烟台 264006)

当前, 富营养化是海湾主要的生态环境问题。全球经济中心大多坐落在环湾地区, 随着陆源污染物入海排放等人类活动压力逐渐增大, 海湾富营养化评估和污染防治成为全球陆海环境综合治理的重点。自1960 年代末特别是1972 年联合国第一次人类环境会议以来, 美国、日本、欧盟诸国相继开始近海富营养化综合治理, 其中, 厘清水体富营养化状态与变化趋势是综合治理的关键, 获取营养盐和叶绿素等富营养化评估要素长期的演变进程高质量数据是前提和难点。在海洋调查监测数据部分年代缺失条件下,数值模型模拟计算是重构营养盐和叶绿素等富营养化评估要素数据集的有效手段, 当前逐步发展建立了EFDC、WASP、FVCOM、ECOHAM、ERSEM、COHERENS 等多种数值模型计算方法(任湘湘, 2011;樊嘉蓉, 2022)。按照DPSIR 理论等, 海湾营养盐、有机质、叶绿素等富营养化评估要素, 主要是环湾地区各种人类经济社会活动压力共同作用的结果。这不仅表现在陆源污染物入海排放等改变氮磷营养盐等化学变量的浓度、结构和组成, 也表现在围填海、海岸工程等通过改变流速、流向等水动力场物理变量, 间接影响营养盐、有机质、叶绿素等生地化变量的时空分布的过程。因此, 海湾富营养化评估要素的数值模拟需要物理、化学和生物多种变量的模拟计算, 其中生物生态过程是长时间的演变过程, 需要进行多变量长时序模拟计算。

海湾富营养化评估要素(营养盐、有机质、叶绿素等)的模拟计算, 相关研究始于20 世纪40 年代(房月英, 2008)。早期研究聚焦富营养盐、叶绿素等单变量季节或时空分布模拟, 至20 世纪80 年代(Skogenet al, 2000; Molletal, 2003), 海洋三维水动力/生物地球化学模型模拟研究进一步完善和发展, 主要包括海湾富营养化要素控制过程及作用机制研究、模型逻辑架构构建及修正、参数率定、模拟计算等。海湾富营养化要素时空变化是水动力输运、生物地球化学迁移转化、浮游植物生长等过程综合作用的结果。对于海湾水动力输运过程研究, 通常采用现场观测结合数值模拟的方法研究流场中潮流和余流变化规律(李高阳, 2015)。对于营养盐生物地球化学迁移转化和浮游植物生长过程, 通常采用陆海同步调查、围隔生态系/培养瓶现场培养实验和数值模型模拟方法研究(Chenetal, 2022)。对于模型逻辑架构, 海湾水动力模型主要包括ROMS、FVCOM、POM、HAMSOM 等, 目前广泛采用ROMS 和FVCOM; 生物地球化学和生态动力学模型主要包括WASP、NPZD、ERSEM 等(刘学海, 2009)。对于模型参数率定, 通常采用伴随矩阵法(吕咸青等, 2002)、模拟退火法(张艳杰等, 2006)、蒙特卡罗法(朱陆陆, 2014)等。对于模拟计算技术, 通常采用粒子追踪法(翁怡婵, 2009)、差分法(Kristeket al, 2010)和分步模拟计算方法(齐宁等, 2020)等。

最近, Gao 等(2020)发展完善了海洋围隔生态系/实验瓶现场培养实验方法, 这不仅揭示了莱州湾浮游植物优势种与不同形态氮迁移转化过程之间的关系, 而且从动力学上量化了莱州湾几种典型浮游植物优势种演替的营养盐控制过程和作用机制。同时,Li 等(2017)根据渤海三维水动力/生物地球化学耦合模型, 应用陆海水环境同步监测调查数据, 将渤海溶解无机氮(dissolved inorganic nitrogen, DIN)等浓度时空分布模拟计算结果的准确性平均提高约25%。这样,围绕渤海富营养化评估要素的数值模拟计算, 针对渤海三维水动力/生物地球化学耦合模型在逻辑架构、过程耦合、参数地域化等存在的问题, 本文在改进渤海三维水动力/生物地球化学过程耦合模型架构基础上, 结合陆海同步调查监测数据校正、参数率定、多变量模拟计算等方法, 建立完善了渤海富营养化评估要素长期演变数据重构技术方法, 据此利用修订的营养状态质量指数(Nutritional Quality Index,NQI), 科学评估了渤海富营养化状态和变化趋势,为进一步提高中国推行的渤海重点海域“深入打好污染防治攻坚战”行动计划的实施提供坚实的科技支撑。

1 研究方法

海湾多变量模拟计算准确性不仅取决于模型的逻辑架构和监测数据的质量, 也取决于模型参数的灵敏度。这样, 海湾富营养化评估要素模拟计算主要包括三维水动力/生物地球化学耦合数值模型构建、参数率定、多变量模拟计算与准确性评价等步骤。

1.1 数值模型构建

王强(2004)和 Dai 等(2015)在 Hamburg Shelf Ocean Model (HAMSOM)(Backhausetal, 1983)基础上建立了渤海三维水动力/生物地球化学耦合模型,本文在此基础上, 对相关生物地球化学模块进行修正。模型中入海口相关数据参见戴爱泉(2015)研究结果。

根据溶解有机氮(dissolved organic nitrogen, DON)的现场培养实验结果(戴爱泉, 2015; Lietal, 2017),不同来源的DON 结构组成和降解速率存在显著差异。本文按照污水处理厂和河流不同来源DON 的降解速率差异性, 将生物地球化学过程模型中DON 状态变量分为污水处理厂输入生物易降解DON(BDON1)、河流输入生物易降解 DON(BDON2)和生物难降解DON(RDON), 方程修正如下:

式中, DON1 为陆源污水处理厂输入的生物易降解DON 浓度; DON2 为河流输入的生物易降解DON 浓度; DON3 为各排污口排入的难降解 DON 浓度;DONm为海源DON 浓度;KrDON为难降解DON 组分占比; KB1、KB2、KB3分别为DON1、DON2、DON3的降解速率常数;KtDON为溶解有机氮微生物降解温度效应系数; temp 为环境温度。

结合生物地球化学模块其他状态变量迁移转化过程, 建立的模型主要框架及过程如图1 所示。模型主要包括浮游植物、浮游动物、溶解无机态营养盐、溶解有机态营养盐、生物碎屑5 个子模块。共包含DIN、陆源生物易降解DON (BDON1, BDON2)、海源DON (BDONm)、生物难降解DON(RDON)、无机磷(DIP)、有机磷(DOP)、浮游植物(PHYTO)、浮游动物(ZOO)、碎屑(DPT)10 个状态变量。物质流统一用氮元素(单位: mmol/m3)表示, 浮游植物生物量可表示为活性颗粒态氮(PPON), chla浓度通过浮游植物叶绿素氮比值关系进行换算(表1)。主要的生物地球化学过程有:浮游植物生长、死亡和营养盐吸收分泌等过程; 浮游动物生长、死亡和营养盐释放等过程; 无机态和有机态营养盐吸收和转化过程; 碎屑矿化和沉降过程等。相关生物地球化学过程动力学方程详见附表Ⅰ。

表1 模型动力学过程所用参数意义及取值Tab.1 Significance and value of parameters used in model dynamics process

图1 渤海三维水动力/生物地球化学耦合模型架构示意图Fig.1 Schematic architecture of three-dimensional hydrodynamic/biogeochemical multi-process coupled model for the Bohai Sea

1.2 模型参数率定

1.2.1 参数率定 参数设置在主要参考文献值基础上, 采用渤海莱州湾海域现场培养实验结果(李克强等, 2007; Lietal, 2009, 2017), 利用DIN、DIP、DON、chla浓度季节变化平面分布模拟结果与实测结果对比进行校正, 利用年际变化模拟结果与实测结果对比进行验证, 模型参数结果如表1 所示。

1.2.2 数据来源

模型校验数据主要取自中国近岸海域环境质量公报(报告书)、中国海洋环境质量(状况)公报、北海区海洋环境公报、中国海洋生态环境状况公报和文献结果, 详见表2。其中, 2019 年氮磷营养盐和chla浓度来源于国家自然科学基金委渤黄海共享航次(航次编号: NORC2018-01; NORC2019-01)调查结果。

表2 模型模拟计算结果验证数据来源Tab.2 Data source for validation on the modelling results

1.3 检验评价

本文选用余弦相似性系数(similarity index, SI;方程 4)和相对标准偏差(relative standard deviation,RSD; 方程5)分别评价模拟值和实测值平面分布趋势相似性以及模拟值偏离实测值的程度(戴爱泉,2015), 从而综合定量评价模型模拟结果的准确性:

其中,Oi和Pi分别表示第i个水质点上的实测值和模拟值;n表示监测站位数量。

针对渤海DIN、DIP、DON 与chla浓度分布场,同时使用Kappa 系数评价模型模拟计算准确性。采用Cohen 提出的Kappa 系数评价分类标准(Cicchettietal,1990; Feinsteinetal, 1990)(表3)。

表3 Kappa 系数评价分类标准Tab.3 Classification criteria for Kappa coefficient evaluation

1.4 富营养化评价方法

渤海海水中溶解有机氮(DON)的浓度近年来持续升高, 某些海域已成为主要的氮形态(李志林等,2015)。在富营养化评价方法中, 用DON 代替化学需氧量(COD)能更好地评估富营养化状态。本文应用营养状态质量指数(NQI)法(陈于望, 1987), 采用DON代替COD 予以修订, 计算公式如下:

式中,C为下标各要素DON、DIN、DIP 与chla的浓度值, 上标*表示各要素DON、DIN、DIP 与chla的营养状态评价标准。其中, chla的营养状态评价标准采用赤潮发生叶绿素浓度阈值, 其他评价要素参照海水水质标准(GB 3097-1997), DON 采用DIN 标准值,各要素营养状态评价标准分别取值为这样, 通过各要素浓度与标准的比值可消掉量纲。NQI＞3 为富营养状态, 3≥NQI≥2 为中等营养状态, NQI＜2 为贫营养状态, NQI 指数越高富营养化程度则越高。

2 结果与讨论

2.1 渤海氮、磷和chl a 浓度时空变化模拟计算结果

2.1.1 DIN、DIP 和DON 平面分布 应用构建完善的三维水动力/生物地球化学耦合数值模型, 模拟计算了2019 年DIN、DIP、DON 春夏秋冬四个季节的浓度平面分布(图2)。经模型参数校正, 再与调查站位上的实测结果对比, DIN、DIP、DON 浓度模拟结果在数值大小及变化趋势上均具有较高的一致性,三者SI 平均值为0.77±0.08, RSD 平均值为26%±5%,Kappa 系数平均值为0.57±0.08。

图2 2019 年全年渤海DIN、DIP 与DON 浓度平面分布模拟计算结果与陆海同步监测调查结果比较图Fig.2 Comparison between the simulated and observed (in land sea synchronous monitoring survey) results of DIN, DIP, and DON in concentration plane distribution in the Bohai Sea in 2019

结果表明, 模拟计算值再现了渤海的DIN、DIP与DON 浓度平面分布春夏秋冬的季节变化特征。在春、夏两季, DIN 浓度因浮游植物生长被消耗利用, 渤海湾与辽东湾南部浓度都比较低, 莱州湾西部受黄河径流输入的影响, 成为夏季渤海DIN 浓度最高的海区(图2b, 2d); 秋季渤海DIN 浓度高值区不断扩大,渤海中部浓度相对较低, 但也高于夏季水平; 冬季渤海DIN 浓度达到最高峰, 浓度高值区较冬季继续扩大, 渤海中部平均浓度达到秋季的两倍水平, 仅有渤海海峡东部浓度较低(图2h)。由于浮游植物在春、夏季生长的消耗, 海水中的DIP 浓度较低(图2j, 2l); 秋季渤海湾南岸、莱州湾以及辽东湾北岸DIP 浓度较高,渤海中部浓度较低, 夏秋季是我国北方降雨期, 三大湾沿岸DIP 浓度上升与陆地径流输入有关, 从而形成了由沿岸高值区向渤海中部逐渐递减的趋势(图2l,2n)(王燕等, 2021); 冬季渤海中DIP 含量上升, 主要因为冬季时浅海风力增大, 导致跃层失去了稳定性,从而促使海水底层颗粒的有机磷分解所产生的DIP随着垂直混合被带到了表层(图2p)(黄爽, 2012)。在春、秋和冬季, DON 高值区主要集中于莱州湾西部、渤海湾、秦皇岛近岸、辽东湾北部及金州湾沿岸地区(图2r, 2v, 2x); 而夏季DON 受高值区出现在渤海中部(图2t), 模型没有很好的模拟再现这一特征, 这与DON 等有机质受浮游植物和微生物调控有关, 由于浮游植物夏季较高的生物量, DON 的源汇过程较为复杂, 模型模拟很难实现精准刻画, 但在调查站位上与调查监测结果比较, 31%的RSD 值还是符合海上调查监测平均误差范围(约30%)的。

在季节变化上, 2019 年渤海夏季DIN、DIP 浓度显著低于其他季节, 而秋、冬季DON 浓度高于其他季节(图2), 这与王婷(2009)在2002～2004 年的研究结果在变化规律上基本一致, 构建的模型很好的体现了黄河调水调沙引起的营养盐输入通量变化对DIN 和DIP 浓度平面分布产生的影响(王婷, 2009; Dingetal, 2020)。

2.1.2 chla浓度平面分布 应用模型模拟计算了2019 年chla春夏秋冬四个季节的浓度平面分布(图3)。在调查站位上与调查监测结果对比, chla浓度模拟结果在数值大小及变化趋势上均具有较高的一致性, SI平均值为0.77±0.04, RSD 平均值为23%±2%, Kappa系数平均值为0.63±0.13。

图3 2019 年渤海chl a 浓度平面分布模拟计算结果与陆海同步监测调查结果比较图Fig.3 Comparison between the simulated and observed (in land sea synchronous monitoring survey) results of chl a concentration plane distribution in the Bohai Sea in 2019

渤海的chla浓度平面分布模拟计算结果体现了浮游植物春夏秋冬的季节变化特征。春季chla浓度的高值区出现在黄河口附近海域、莱州湾南部、渤海湾与辽东湾北部、秦皇岛外海与金州湾, 渤海中央以及辽东湾南部海区为低值区(图3b); 夏季chla浓度高值区主要集中于黄河与辽河入海口附近海域, 中央海区浓度较春季增加一倍(图3d); 秋季由于光照、温度等因素影响, chla浓度急剧下降, 仅在辽东湾沿岸出现高值区(图3f); 冬季整个渤海的表层chla浓度继续下降, 只在渤海湾沿岸与秦皇岛沿岸地区浓度较高(图3h)。

在季节变化上, 2019 年渤海夏季chla浓度显著高于其他季节, chla与DIN 和DIP 浓度的季节性分布趋势呈负相关, 与DON 浓度季节性分布趋势呈正相关(图3)。主要的浓度高值区分布在黄河口临近海域,构建的模型很好的体现了黄河调水调沙引起的营养盐输入通量变化对浮游植物生物量平面分布产生的影响(Dingetal, 2020)。

综上, DIN、DIP、DON 与chla浓度模拟值与实测值的相对标准偏差、相似性系数和Kappa 系数平均值分别为24%、0.77 和0.60, 说明通过该模型能够较真实的反映渤海海域DIN、DIP、DON 与chla浓度分布情况, 可以用于渤海DIN、DIP、DON 与chla浓度长期变化模拟计算和数据重构。

2.1.3 年际变化 以1980 年至2020 年部分年代调查监测的DIN、DIP、DON 和chla浓度年际变化对模型进行验证。图4 中模拟值为模型对全渤海逐月模拟结果的平均值, 监测值为对应某一月份的所有监测值的平均值。模拟结果与调查监测结果比较表明,模拟结果在数值大小与变化趋势上均能够很好的再现渤海近40 年营养盐和chla的年际变化, DIN、DIP、DON 和chla年际变化模拟结果与调查结果的RSD平均值小于20%, SI 平均值高于0.80(图4)。

模拟结果表明, 20 世纪80 年代至90 年代中期DIN 浓度显著增加, 1995 年冬季DIN 浓度达到最大值(24.4 μmol/L), 在1996 年至2004 年间, DIN 浓度逐渐下降, 在2005 年至2013 年期间, DIN 浓度再次呈现上升趋势, 近几年DIN 浓度又呈下降趋势(图4a), 构建的模型能很好地再现这两次DIN 浓度峰值变化(王修林等, 2006; Wangetal, 2019)。分析表明, 农业化肥使用、污水排放等人类活动导致总氮河流输入和大气沉降通量增加, 是造成渤海DIN 浓度升高的主要原因, 自2002 年起水利部黄河委员会于夏季实施的调水调沙工程, 是导致渤海夏秋季营养盐输入和浓度增加的原因(Dingetal, 2020)。DIP 浓度在20 世纪80年代至90 年代初逐渐下降, 在1992 年夏季达到最小值(0.1 μmol/L), 然后增加至90 年代后期, 在1998 年冬季DIP 浓度达到最大值(1.1 μmol/L), 之后呈下降趋势(图4b), 模拟结果在年际变化和季节变化上均与文献结果基本一致(Ningetal, 2010; Xinetal, 2019)。20 世纪90 年代后期DIP 浓度的下降, 与含磷洗涤剂的禁用有关。在20 世纪80 年代到90 年代初DON 浓度总体处于缓慢上升趋势, 直至1994 年浓度有所下降;21 世纪10 年代末有所上升, 在2014 年又降低, 之后呈上升趋势(图4c), 这与陆源长期输入变化与海源浮游植物、微生物长期演变过程有关(唐永等, 2017)。

Chla浓度自20 世纪80 年代至21 世纪初呈波动缓慢升高趋势, 在年际变化上与营养盐浓度, 特别是与DIN 浓度变化呈现很好的相关性, 这与孟庆辉等(2022)的研究结果一致。Chla浓度季节变化主要是夏季峰(图4d), 与DIN 浓度呈反向变化关系, 这与夏季藻华对营养物质的消耗和冬季浮游植物死亡营养盐恢复的生态学规律一致; 与DON 浓度呈复杂源汇变化关系, 说明DON 浓度主要受陆源和海源的影响,也是支持浮游植物生长的重要物质(唐永等, 2017)。

2.2 模拟结果准确性评价

为验证模型改进后对模拟计算结果准确性的提高程度, 将改进前后的模拟计算量化验证结果进行比较, 结果如图5 所示。氮、磷营养盐与chla浓度的相似性系数平均值分别提高了0.07 和0.13, Kappa系数分别提高了0.09 和0.10, 说明改进后的模型对氮、磷营养盐与chla浓度的变化趋势的模拟要显著比模型改进前与实测结果更加吻合。二者的相对标准偏差也均有明显降低, 氮、磷营养盐与chla浓度分别降低了约11%和6%, 说明在数值大小上, 模型改进后的结果与实测结果也更为接近。综上所述, 通过参数优化和优质数据选取对模型进行改进, 提高了模型模拟计算的准确性。

图5 模型改进前后渤海氮、磷营养盐与chl a 浓度模拟计算评价结果的Taylor 图Fig.5 The Taylor diagram of the simulation results of N-P nutrients and chl a concentrations in Bohai Sea before and after model improvement

2.3 渤海富营养化评估

2.3.1 渤海富营养化状态年际变化 渤海富营养化状态指数(NQI)年际变化计算结果表明, 自1980 年至2020 年, NQI 整体呈先上升后下降的倒“U”型变化特征(图6)。其中, 1980 年至1995 年, NQI 上升达到极大值2.27, 1995 年至2013 年NQI 缓慢下降, 2013 年至2020 年迅速下降; 从数值上看, 1993 年至2017 年, NQI基本都大于2, 富营养状态处于中等水平, 当前处于贫富营养状态。利用渤海复合富营养化指数(CEI)(Linetal,2020)年际变化计算结果进行后验, 结果表明, NQI 与CEI 在变化趋势上相当一致(SI=0.92) (图6)。

图6 1980 至2020 年渤海复合富营养化指数(CEI)与富营养化状态(NQI)年际变化图Fig.6 Interannual changes of Bohai composite eutrophication index (CEI) and eutrophication status (NQI) from 1980 to 2020

2.3.2 渤海富营养化状态季节变化 在渤海全海域富营养化状态季节变化规律上, NQI 呈倒“U”型变化特征, 夏季富营养化程度最高; NQI 在5～10 月份超过2.0, 处于中等富营养状态。利用渤海复合富营养化指数(CEI)(Linetal, 2020)季节变化结果对比后验,NQI 与CEI 在变化规律上也较为一致(SI=0.74)(图7)。

图7 渤海复合富营养化指数(CEI)与富营养化状态(NQI)季节变化图Fig.7 Seasonal change in composite eutrophication index (CEI)and eutrophication status (NQI) in Bohai Sea

2.3.3 渤海富营养化状态平面分布 渤海富营养化状态指数(NQI)平面分布模拟计算结果如图8 所示。渤海富营养化状态平面分布总体呈现近岸较高、中部较低的趋势, 这与Lin 等(2020)的研究结果一致。近岸较高的富营养化状态主要与人类活动有关, 且海湾近岸地区水体交换能力差, 二者综合导致近岸海域富营养化程度显著高于中部海域。在季节变化上,富营养化状态与chla浓度模拟结果(图3)变化规律基本一致, 即随chla浓度春夏季升高, 富营养化水平升高, 随chla浓度秋冬季下降, 富营养化水平下降,同时还受海洋中的环境因子(光照、温度、营养盐、混合层深度等)季节变化调控。

图8 2019 年渤海富营养化状态(NQI)平面分布图Fig.8 Horizontal distribution of eutrophication status (NQI) in 2019 in the Bohai Sea

3 结论

本文通过改进生物地球化学模型架构, 建立完善了用于重构富营养化评价要素长期演变数据的渤海三维水动力/生物地球化学耦合数值模型。利用构建的模型成功模拟了渤海近四十年来DIN、DIP、DON与chla浓度时空分布[(相对标准偏差平均值)为24%,(余弦相似性系数平均值)为0.77,(Kappa系数平均值)为0.60], 说明模型相关动力学方程与参数能够反映渤海海域的地域化特征。应用模型模拟结果计算了营养状态质量指数(NQI), 与复合富营养化指数(CEI)具有较好的一致性(=0.83), 说明该模型能够反映渤海富营养化状况。

本文构建的三维水动力/生物地球化学耦合模型,对富营养化状态的模拟计算结果能够较真实的反映渤海海域富营养化状态的时空变化情况, 建立的富营养化评估方法具有较高可靠性, 可以用于渤海富营养化评估, 并支撑渤海富营养化综合治理。