闽江河口富营养化数值模拟

2022-10-08肖桂荣林金城

海南大学学报（自然科学版） 2022年3期

高伟，肖桂荣，林金城

（福州大学数字中国研究院〈福建〉，福建福州 350108）

富营养化是指由于人类活动，过多的氮、磷等营养物质进入水体，在适当的环境条件下，导致浮游植物大量繁殖，使原有生态结构改变、功能退化，从而引起水质污染现象［1］. 河口地区是海洋与河流之间的连接枢纽，受经济快速发展和人口数量不断增长的影响，河口地区往往负荷着流域沿岸工农业、城镇生活等产生的大量废水，最终流入海洋，导致河口及周边海洋生态系统遭到不同程度的损害，水体富营养化态势日趋严重，对地区的可持续发展造成了不可忽视的影响. 因此，越来越多的学者对河口水环境污染开展了相关方面的研究.

在有关河口水环境问题的研究中，水动力因素起着至关重要的作用. 水动力是物质输运的基础，对污染物的输移扩散起着决定性的作用. 因此，水动力数学模型的建立是研究水体中物质输运的基础. 水质模型则是描述水体中各类非保守物质，如污染物、营养盐等物质的输移、变化规律的数学描述. 数值模拟是研究水体水动力和水质的一种主要方法和有力工具，采用数学模型来模拟特定的物理现象，通过高效的计算机利用数值计算法进行近似求解，以此来模拟复杂的物理现象、复演自然演变的物理过程、预报未来演变的效果. 在有关闽江河口水质方面，已有许多学者进行大量研究. 侯昱廷等［2］对闽江河口营养盐的季节性变化进行了分析；何露露等［3］对闽江河口潮汐沼泽湿地地区的营养盐分布特征进行讨论；Fan 等［4］研究了闽江河口溶解无机氮的来源与排放；林秀珠等［5］对闽江口及其近岸海域水质进行分析评价. 以上研究成果是基于实测数据进行空间分析与富营养化评价，计算方便，数据量要求不高，但欠缺对河流径流、潮汐、地形、气候、生物等多种因素导致的营养盐限制的考虑［6］，且在连续时间尺度上的研究局限性较大. 采用数值模拟的方法不存在缩尺效应，效益高、易于控制与优化等特点，能够考虑多种因素共同作用，可以更综合、更全面、更精细化地掌握目标海域的环境因素和富营养化程度，为赤潮等生态灾害的防御提供有效的服务管理决策.

目前水动力数值模拟计算模式的发展也较为成熟，常用的有美国麻州大学的FVCOM 模式、普林斯顿大学的POM 模式和ECOMSED 模式等. 其中ECOMSED（Estuarine，Coastal and Ocean Model System with Sediments）模式是当前较为广泛使用的海洋模式，以POM 模式为基础发展起来并继承了POM 模式的一些优点，在近岸海域水动力数值模拟方面具有明显的优势［7］.RCA（Row and Column of Aesop）是Hydroqual公司基于WASP 和AESOP 模型开发出的三维水质模型. 该模型拥有完整的富营养化模块，可应用于河口和沿海地区，在国内外应用较为广泛. 如刘广州等［8］将RCA应用于我国珠江口的水质模拟，以及在美国切萨皮克湾［9］、Mamala 湾［10］和缅甸Mass 海湾［11］等地区均有应用. 此外，RCA 与ECOMSED 模型可以相兼容，共用相同的计算网格，并且RCA 还可以直接利用水动力模型输出的流场信息. 张明耦合ECOMSED 与RCA构建大伙房水库三维模型，并取得较好的实验效果［12］.

综上所述，笔者以闽江河口为研究对象，耦合ECOMSED 与RCA 模型，通过对模型进行率定与验证，构建三维水动力水质耦合模型，模拟闽江河口水质及富营养化状况，研究闽江河口的溶解无机氮（DIN）、溶解无机磷（DIP）、化学耗氧量（COD）、溶解氧（DO）的表层浓度平面分布，并对海水富营养化状况进行分析评价.

1 研究区概况

闽江河口位于福建省东部，地理位置在东经119°20 ′～119°55 ′和北纬25°50 ′～26°18 ′之间，闽江是福建省最大的入海河流，其干流全长559 km，流经福建省北半部36 个县、市，流域面积达624 km2，多年平均入海径流量66 960 m3. 闽江河口区域不仅是东南沿海一带居民生产生活重要的饮用水水源地，同时也是福建省以及周边沿海区域重要的水上交通要道之一. 近年来由于闽江流域经济社会的快速发展以及人类活动的不断影响，大量的污染物沿着闽江流入台湾海峡，给闽江河口及周边近岸海域的生态环境带来了极大的压力.

2 研究方法

2.1 ECOMSED模型选取ECOMSED 模式［13］，通过建立三维斜压水动力模型，对闽江河口的水动力特性进行数值模拟，为水质模拟提供基础. 考虑到闽江河口存在大面积的淤泥质潮滩，模型计算区域因潮间带上的水陆边界随着涨落潮的移动而发生变化，但ECOMSED 无法完全适用此类计算，故对模型进行改进，加入动边界处理技术. 通过动边界处理技术的引用，使得模型的每次计算都相应判断一次水深，若发现网格点水深小于临界水深，则判定为干网格，反之判定为湿网格. 在模型运行过程中，只有湿网格参与计算，干网格不参与计算，并且其流速为零. 临界水深Hdry大小可根据潮间带坡度和时间步长来确定，取0.4 m［14］.

2.2 RCA模型RCA是一个较为先进、成熟的三维水质模型［15］. 如图1所示，该模型包括浮游植物、氮循环、磷循环、硅循环、碳循环和溶解氧等多个相互作用的系统，以及26 个状态变量［16］，并针对水体中26 个参量（如表1 所示）建立平衡方程. 该模型能够模拟参量在水体内的对流和扩散过程，以及物理、化学、生物作用和点源、非点源、不连续线源、大气沉降的输入过程，其质量守恒方程可以描述为

图1 RCA水质模型概念图表1 RCA水质模型状态变量

其中，c为水质参量的浓度，单位：mg·L-1；t为时间，单位：s；E为扩散系数，单位：m2·s-1；U为对流速度，单位：m·s-1；S为化学、生物以及物理过程产生的源或汇，表示水质参量的动力学相互作用过程，单位：mg·L-1·s-1；x，y，z分别表示横向、纵向和垂向坐标.

2.3 富营养化状况评价水体富营养化状况评价采用《近海海域环境监测规范》［17］提出的富营养化指数法，其计算公式

其中，DIP为溶解无机磷，DIN为溶解无机氮，COD为化学需氧量，单位：mg·L-1；E为富营养化指数，E的大小代表水质富营养化程度的高低.1≤E＜2 时，水质处于轻度富营养化状态；2 ≤E＜5 时，水质处于中度富营养化状态；5 ≤E＜15时，水质处于重度富营养化状态；E≥15时，表示水质处于严重富营养化状态.

状态变量Salt PHYT1 PHYT2 PHYT3 RPOP LPOP RDOP LDOP PO4T RPON LPON RDON LDON参数说明盐度硅藻绿藻蓝绿藻难降解颗粒有机磷易降解颗粒有机磷难降解溶解有机磷易降解溶解有机磷藻类磷+溶解无机磷难降解颗粒有机氮易降解颗粒有机氮难降解溶解有机氮易降解溶解有机氮单位psu mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1状态变量NH4T NO2+NO3 SIU SIT RPOC LPOC RDOC LDOC ExDOC RePOC ReDOC O2 DO参数说明藻类氮和氨氮亚硝酸盐氮和硝酸盐氮生物硅藻类硅和可吸收硅难降解颗粒有机碳易降解颗粒有机碳难降解溶解有机碳易降解溶解有机碳藻类渗出的溶解有机碳活性颗粒有机碳活性易溶解有机碳溶解氧当量溶解氧单位mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1 mg·L-1

2.4 模型构建与验证

2.4.1 模型构建考虑到闽江河口河网复杂，岸线曲折，滩涂面积大等特点，模型计算范围东西向约为119°23 ′～119°50 ′E，南北向约为25°55 ′～26°18 ′N.ECOMSED 模型采用正交矩形网格，网格分辨率为200 m×200 m，网格数共计222×208个，垂向采用sigma坐标，等间距分为11层，能够较好地反映了区域岸线及地形变化. 模型中水深数据从海军航保部的闽江口海图中提取，水动力模型内模时间步长设置为5 s，内外模分裂比为5.

模型在外海边界处采用分潮调和常数，其源于美国俄勒冈州立大学建立的区域潮汐潮流模型，取S2，M2，N2，K1，P1，O16 个分潮；上游边界条件设置来源于闽江竹岐水文站的实测月平均径流量，模型中的初始场以及开边界条件采用HYCOM 每隔24 h、精度0.08°的全球再分析温盐数据；模型中的气象驱动条件设置采用NCEP 间隔6 h一次的风速、风向以及海表热通量等气象数据. 水动力模型模拟时段为2014年6月11日～2014年12月31日，其中2014年10月～2014年12月的流场数据用于驱动RCA水质模型.

RCA 水质模型时间计算步长为5 min. 河流营养盐负荷主要参考福建省生态环境厅闽安站水质监测数据及前人的相关调查研究成果［2，18-19］，各水质因子（包括溶解氧、氨氮、硝酸盐、亚硝酸盐、磷酸盐、有机磷、有机氮及有机磷）在外海边界的设置主要参考于寒［20］在闽江口的相关研究成果，沉积物通量模型中的颗粒有机碳、颗粒有机氮、颗粒有机磷以及各营养盐含量的初始设置参照前人对闽江口沉积物中有机物含量的研究进行取值［19，21-22］. 水质模型各主要反应参数取值参考文献［20］进行取值.

2.4.2 模型验证用于水动力模型及水质模型验证的实测数据包括：1）河海大学于2014年6月14日～15日在闽江河口进行的水文观测任务，包括2 个潮位观测点、4 个潮流观测点，监测频率为30 min；2）福建省海洋研究所于2014 年11 月13 日进行的水质监测调查，水质监测数据包括DO，COD，DIN 和DIP等. 闽江河口主要潮位站、潮流站及水质测站分布如图2所示，以上实测数据基本满足模型验证的需求.

图2 闽江河口主要潮位站、潮流站及水质测站分布

水动力数值模拟的验证结果见表2～4. 由验证结果可见，潮流观测点的涨落潮平均流速的偏差值较小，相对误差均小于10%，其对应的流向误差小于10°，潮位站平均高、低潮位的偏差值在±10 cm 以内，均符合技术规范［23］对潮流场的验证要求. 结果表明选取的参数合理，构建的模型能够反映出闽江河口的水动力特征.

表2 P1～P4潮流站涨潮时流速、流向验证情况表

在模型构建之初，需要在初始场及开边界处进行温度条件设置，通过对温度进行模拟验证，进一步提高所构建模型的准确性. 闽江河口海域水温验证采用固定浮标黄岐8 201 站点（119.55.31 °E，26.23.15°N）的数据，观测时间为2014年6月～2014年12月，数据由文献［20］提供. 如图3所示，模拟精度良好，表明模型能较好地模拟研究区域水温变化情况.

图3 监测站位温度验证图

水质测站的验证结果如表5 和图4 所示. 从总体上看，各水质测站的模拟值与实测值在大小、空间趋势变化基本上是一致的. 各水质因子的相对平均误差，除DIP与COD误差稍大，DO与DIN模拟效果较好.造成上述现象的主要原因：1）由于河口中水质过程较为复杂，水质模型对其进行了一定的简化，从而导致部分计算值与实测值在一些地区存有较大偏差；2）因为DIP与COD在一些测站的数值验证上就已存在较大差异，这必然导致水质因子总体平均误差偏大. 因此，误差较大的原因主要是缺少完整准确的有机物数据，而且水体生化过程较为复杂，未能完全反映实际生化过程. 同时闽江径流营养盐输入采用给定月均值，边界数据时间分辨率较粗，在一定程度上，致使营养盐的计算精度受到影响. 另外，在源强分析时未考虑陆源面源输入也是造成偏差的重要原因.

图4 DIN，DIP，COD和DO在水质测站表层验证结果

表3 P1～P4潮流站落潮时流速、流向验证情况表

表4 T1和琯头站潮位验证情况表

表5 水质因子在模拟值与实测值之间的相对平均误差

目前，在RCA 水质模型模拟精度的评价标准方面，已有研究表明，当水质因子的相对平均误差小于35%时，模拟结果是可以接受的［16］. 从各水质因子的总体验证结果来看，相对误差均低于25%，表明了闽江口水动力水质耦合模型的模拟值与实测值拟合良好，基本能够反映闽江口水质变化过程中的一些主要规律.

3 结果与分析

河口海域因海陆交汇的原因，导致不同区域可能呈现不同的富营养化特征. 因此，根据空间区域的差异性将研究区域进行分区评价. 参考ASSETS 方法，将研究海域按照盐度差异划分为三类评价单元区域：感潮淡水区（SALT＜0.5），混合区（0.5≤SALT＜25），海水区（SALT≥25），如图5 所示. 基于ECOMSED 模型模拟研究区域的盐度变化，且由于河口区盐度随涨落潮时间变化显著，因此取一日盐度平均值作为划分单元标准.

图5 闽江河口盐度分布及评价子区域划分

3.1 水质因子模拟结果分析基于上述已验证的RCA水质模型的模拟结果，分析闽江河口的DIN，DIP，COD和DO的表层浓度平面分布（如图6所示），各分区的水质因子浓度平均数值见表6.

图6 闽江河口DIN，DIP，COD，DO和E的浓度场分布

水中的无机氮DIN主要包含硝酸盐氮、亚硝酸盐氮和氨氮3种营养盐，浮游植物可以吸收氨氮和硝酸盐氮并用于自身的细胞生长. 因此，当营养盐含量过高时，水体会呈现富营养化，出现藻华、赤潮等现象.由图6a 可见，在闽江口海域，DIN浓度场从感潮淡水区的1.26 mg·L-1向海水区递减至0.8 mg·L-1左右，平均值为0.87 mg·L-1，越靠近外海其浓度越低. 闽江河口的无机氮浓度超标严重，基本处于劣四类标准海水水质. 其中，感潮淡水区与混合区的DIN表现出较高的浓度值，研究证实近年来福州社会经济的快速发展，闽江口近岸区域大量NH+4输入闽江发生硝化反应导致DIN含量升高［24］. 在空间分布上，可见浓度场等值线偏向东南方向递减，且闽江北支水道的DIN含量低于南支水道，主要是由于北支水道的流速较快，水体交换能力较强，对其产生了积极的扩散输运作用.

磷是水体中藻类生长所需的重要元素，过量排放磷酸盐同样会导致浮游植物大量生长. 如图6b 所示，闽江河口磷酸盐浓度在0.010～0.02 mg·L-1之间，平均浓度在0.011 mg·L-1，符合海水水质第一类标准的要求. 磷酸盐的浓度场空间分布与无机氮浓度分布近似，其随闽江径流排入闽江口海域，均为由河口内向东南方向的外海逐渐降低. 整体上，磷酸盐含量表现为由河口内向外海逐渐递减的趋势，可见闽江河口中的磷酸盐浓度受海水稀释作用较为显著.

化学耗氧量COD与水体受有机物污染的程度密切相关，COD含量的高低可以反映水体受有机物污染程度的大小. 由图6c 可知，COD浓度变化范围在1.9～3.6 mg·L-1之间，平均浓度为2.64 mg·L-1，基本处于二类海水水质标准范围内. 闽江河口海域COD的浓度分布，总体上呈现出河口沿岸较高而远离河口近岸较低的分布特征. 由于COD主要来源于闽江径流，且在感潮区内水道狭窄，流速较快，水体交换条件好，溶解氧充沛，因此表现出浓度较混合区低的情况. 此外，在闽江河口滩涂区域，COD浓度值偏高，一方面由于闽江河口海域的滩涂区域广阔，水体流速慢，致其水体迁移扩散能力弱，从而导致污染物滞留时间长，易堆积；另一方面，由于模型没有考虑敖江径流营养盐负荷，因此如图7所示，在b 和c 区域易导致有机物浓度较高，COD也就越高.

在衡量区域水体自净能力时，溶解氧DO含量可作为一项重要的参考指标，因为其既是水体中的重要营养元素，同时也参与水中绝大多数的生物化学过程. 闽江河口溶解氧含量的变化范围为5.2～8.1 mg·L-1，平均值为6.62 mg·L-1，符合海水水质第一类标准要求. 如图6d所示，闽江河口DO含量较为充沛，结合图2 中各监测站点位置能够看到各站位的变化幅度较小，整体分布相对均匀，呈现向外海逐渐递减的趋势.

3.2 富营养化评价在前文研究基础上，根据富营养化指数法对闽江河口富营养化状况进行评价. 评价结果如表6所示，富营养化指数E范围在4～11之间，表明闽江河口基本处于重度富营养化状态.

表6 各分区的水质因子浓度平均数值

如图6e 所示，从分区上看，感潮淡水区与混合区属于重度富营养化状态，海水区富营养程度较轻，属于中度富营养化状态，与盐度25 psu等值线分界一致. 结合图7可以看出，富营养化指数较高的区域大多出现在闽江河口的主要工业区与浅滩区域. 尤其在图7c区域，富营养化指数明显高于周围海域.

图7 富营养化指数E较高区域

造成这种状况的原因一方面是该区域为浅滩地区，平均水深仅为1.2 m，水体交换能力差；另一方面与沿岸人类活动影响有关，即生活以及农业污水的排放造成的富营养化指数偏高［24］. 在整体空间分布上，富营养化指数呈现由河口内向河口外海区域逐渐递减的阶梯状分布，且近岸水域指数大于远岸水域指数. 出现上述分布情况，与闽江河口独特的五口入海地形和潮流特征影响密切相关，同时闽江径流污染源的输入也不容忽视. 再结合水动力模拟结果来看，闽江河口潮流基本趋势以径流为主要动力，以落潮余流为主，潮流整体偏向东南处流动. 由于东北部海域在涨落潮时流速相差无几，且流速较慢，因此，受流场的作用，污染物更容易向东南方向海域输运，这也是DIN，DIP，COD以及富营养化指数在东南部海域均表现出较高的数值的原因. 此外，闽江河口拥有较高含量的溶解氧，这在一定程度上可减轻该海域的富营养化程度.

4 讨论

采用数值模拟的方法，主要是基于易控制与优化，能够考虑多种因素共同作用等优点. 虽然通过空间分析等方法也可用于水质研究，但该类方法较高依赖实测数据. 考虑到实测数据的缺乏，通过构建模型以达到更好的适应性. 水动力及水质模拟验证结果均在合理的范围内，表明本文模型是可靠的.

根据无机氮的模拟结果可以看出其整体变化趋势呈现从河口内向外海逐渐递减，与叶翔等对闽江河口无机氮的浓度变化特征研究结果相一致［25］. 此外，DIN含量严重超标，初步分析是受沿岸农业施肥和工业排放等的剧烈影响所导致. 磷酸盐主要来源于沿岸工业与城镇的污水排放，从DIP模拟结果来看，DIP的含量基本维持在一个较低的水平，其表面浓度分布趋势与DIN近似，这与深圳湾［26］、长江口［27］等近岸河口区域的无机磷分布变化特征的研究结果相似.COD平均浓度值为2.64 mg·L-1，但相较于林秀珠对闽江河口COD的研究调查数据偏高，这是由于研究方法的异同以及水动力环境的复杂变化，造成与前人的研究结果有所差异，但在整体上的变化趋势是一致的［5］. 溶解氧模拟结果的平均值为6.62 mg·L-1，与何露露等对闽江河口的实测调查数据［溶解氧平均值为（6.38±0.53）mg·L-1］基本吻合［3］，这也进一步验证了模拟结果的可靠性.

对于闽江河口的富营养化评价结果显示水体富营养化严重，水质较差，实验结果同样证实了郑小宏等对闽江河口的富营养化评价［18］，模拟结果显示水体富营养化指数由河口内向河口外海区逐渐递减，受到海水稀释作用较为显著，林秀珠的研究结果也同样证明该区域水质的整体变化趋势［5］.

构建的水动力水质模型数值模拟结果不仅与实测数据基本符合，而且与前人的研究结论也基本一致，进一步说明了所构建模型的可靠性，同时也验证了通过数值模拟的方法预测水质变化也是一种有效的手段.