辽东湾浮游植物分布的生态系统动力学模型

王昆，李楠，宋广军，关春江，胡超魁，杜静，邵泽伟，吴金浩*

(1.辽宁省海洋水产科学研究院，辽宁 大连 116023； 2.国家海洋环境监测中心，辽宁 大连 116023)

辽东湾位于渤海的最北部，总水域面积约为10 000 km2，水深较浅，地形多变，潮汐类型复杂。多年来，人们针对辽东湾海域进行了大量的水文动力、泥沙输运和水环境保护等方面的数值模拟研究[1-4]，特别是在2014年以来国家开展的渤海综合治理过程中，相关部门在辽东湾水环境保护研究方面取得了较多的成果。

水体中的浮游植物含量是海洋初级生产力的主要贡献指标，是水体中贝类、牡蛎等生物饵料(单细胞藻类)的主要贡献体，在海洋生态系统和环境保护中起着至关重要的作用。而叶绿素a的浓度往往用来定量表征海水中浮游植物的丰度，在海洋初级生产力的计算与评价中，叶绿素a的浓度是主要的决定性因素，故很多研究中常将叶绿素a浓度作为主要水体研究指标，并探讨其与水动力的响应关系[5-7]。相关研究中，徐文龙等[8]基于光学观测数据，采用高精度反演和数理统计手段对南海东北部海区夏季叶绿素a浓度的垂向变化特征及其对水动力过程的响应进行了研究；He等[9]基于南海布放生物地球化学浮标数据，指出南海中央海盆存在一定深度的叶绿素极大值分布区。在辽东湾海区的相关研究中，Pei等[10]于2013年夏季，在辽东湾进行了一项实地调查研究，以确定该海湾浮游植物的动态及其初级生产力受环境制约的程度，结果显示，辽东湾的区域生产受到温度和光照限制的综合控制；2016—2017年，田思瑶等[11]在辽东湾中部近岸海域，采用单因子标准指数法和综合指数评价法对其COD、石油类、叶绿素分布特征及富营养化状态进行了评价，结果显示，调查海域内COD和石油类指标分布呈河口及沿岸浓度较高的特点，存在明显的季节变化特征，在监测海域中、北部以轻度、中度富营养化状态为主，富营养化程度较高的站位为河流入海区部分站位；2018年，裴绍峰等[12]采用观测与模拟分析相结合的手段，对辽东湾海区的叶绿素a和初级生产力进行研究，结果显示，叶绿素a低值区出现在双台子河的河口，其原因是河水中过量泥沙悬浮物降低了该海区水体的透明度，从而导致浮游植物生长受到光的限制。综上可以看出，前期大部分研究都是基于监测数据，采用数理统计方法来评价水质指标在水体中的时空分布特征，即使考虑多指标的水质模型，也假定为各指标间不发生任何生化反应的保守型问题，但实际水体存在温度、盐度和光照等因素的影响，各个指标间会发生复杂的生物化学反应过程，若按单一保守型的变化响应来理想化描述，其结果的呈现往往偏离实际物理规律的本质。为了解决这一问题，本研究中在辽东湾水文动力数值模型构建的基础上，引入水体中叶绿素a、溶解氧、氨氮、亚硝酸氮、硝酸氮和活性磷酸盐各指标间的生化反应过程，建立辽东湾海域以叶绿素a浓度为主要未知量的生态系统动力学模型。该模型结合辽东湾海区的水动力条件，基于实际观测信息，综合考虑可能影响到整个生态系统的物理、化学和生物变化过程，模拟预测多工况下叶绿素a浓度的分布结果，并用2020年辽东湾的观测数据对该模型进行验证，以期为辽东湾海区的生态系统动力学过程分析提供合理方案，为水体富营养化和赤潮预测预报研究提供更为科学有效的技术手段。

1 水动力控制方程

利用二维水动力模型对辽东湾海域的水动力场进行模拟，以期数值再现和预测该海域水体运动过程，为浮游植物在水体中的分布模拟提供背景场[13]。基于不可压缩流体Reynolds平均Navier-Stokes方程，将其水平动量方程和连续性方程在总水深[0,h]范围内进行积分后，可得如下的二维深度平均浅水方程组。

连续性方程：

(1)

动量方程：

(2)

(3)

2 边界条件

给定两种边界条件，即闭边界和开边界条件。

1)闭边界条件。闭边界即水陆交界边界，一般由海岸线及海岛确定，在闭边界处法向流速为0。

2)开边界条件。本模型中开边界处采用如下的分潮调和形式计算：

(4)

其中：m为考虑的分潮个数；ωi为第i个分潮的角速度；fi、ui分别为第i个分潮的交点因子和迟角修正；Hi和gi为调和常数，即分别为各分潮的振幅和迟角；V0i为第i分潮的天文初相位。

3)干湿动边界。对于潮滩，水陆交界的位置随着潮位涨落而变化，本模型中考虑了动边界内网格节点的干湿变化[14]。

4)初始条件。计算开始时采用“冷态”启动，即各物理量初值均为0。

5)离散方法。模型对计算区域的空间离散采用有限体积法，对不同的计算区域采用非结构三角形网格剖分形式，大大增强了系统对岸线和地形变化的适应性，提高了计算精度。

3 生态系统动力学模型

1)描述叶绿素a浓度的生态系统动力学模型。控制方程为

(5)

其中：u、v为水动力方程中的垂向平均流速；Ex、Ey分别为x、y轴两方向的紊动扩散系数；S为源汇项；F为生化反应项。

(6)

其中： Chlacs为浮游植物生长过程对叶绿素a的产生过程项；Chlahx为呼吸过程对叶绿素a的改变过程项；Chlasw为死亡过程对叶绿素a的改变过程项；Chlacj为沉降过程对叶绿素a的改变过程项。

2)影响叶绿素a分布的4个主要生态过程。

浮游植物生长过程引起的叶绿素a浓度变化的表达式：

Chlacs=Gf×a1×a2×GN。

(7)

其中：Gf为光合作用过程函数；a1、a2为浮游植物生长速率[1](表1)；GN为营养限制因子[15]，且

GN=Amax×e(-1.6dz/θgs)×Φd/dz。

(8)

式中：Amax代表海区的正午最大产氧量[16]；dz为水层厚度；θgs为光衰减系数[17]；Φd为太阳辐射强度的日变化函数[18]。

呼吸过程引起的叶绿素a浓度变化的表达式：

Chlahx=Of×a1×a2×GN。

(9)

其中，Of为浮游植物呼吸过程的需氧量函数，且

(10)

式中：CDO为浮游植物呼吸过程中所需的溶解氧半饱和浓度；sfz为浮游植物呼吸速率[19]；a5为温度依赖系数；T为水温(℃)。

死亡过程引起的叶绿素a浓度变化的表达式：

在相同作业温度下，更高的催化剂活性能够提升苯乙烯单体的产量。高选择性则能生成更高比例的苯乙烯单体，降低副产物的生成，如苯和甲苯等。当前可用的苯乙烯单体催化剂只能在高活性和高选择性之间二选一。这是因为在超低SHR作业条件下，催化剂活性的上升通常会导致选择性下降。

Chlasw=a3×Chla。

(11)

其中，a3为浮游植物的死亡率[1]。

沉降过程引起的叶绿素a浓度变化的表达式：

Chlacj=a4/dz×Chla。

(12)

其中，a4为沉降速率[20]。

3)主要生态过程求解需要的6个子系统[21]。

生化需氧量子系统1：

(13)

其中：JBOD为BOD的降解函数，与BOD的降解速率、温度依赖系数和溶解氧半饱和浓度有关。

氨氮子系统2：

(14)

其中：SN为氨氮释放BOD的衰减函数，与BOD衰减释放的氨氮比例与降解有关；FfN、FxN分别为浮游植物和细菌对氨氮的吸收函数，与水深、吸收的氨氮量及溶解氧半饱和浓度有关；Sx为硝化过程中氨氮转化为亚硝酸盐的硝化速率函数，同样与相应过程的衰减速率、温度依赖系数和溶解氧半饱和浓度有关。

溶解氧子系统3：

(15)

其中：Fat为大气复氧函数，与大气复氧速率和溶解氧浓度有关；Fx为需氧量的硝化速率函数，与单位质量的氨氮到亚硝酸盐、单位质量的亚硝酸盐到硝酸盐硝化作用的需氧量，硝化作用中20 ℃时的一级衰减速率，以及亚硝酸盐到硝酸盐衰减过程的温度依赖系数等有关；Fsed为沉积物的需氧量函数，由底层碎屑的半饱和浓度、单位面积沉积物的需氧量及沉积物衰变的温度依赖系数控制。

亚硝酸盐子系统4：

(16)

硝酸盐子系统5：

(17)

其中：K2、K3分别为温度20 ℃时水体的一级硝化和反硝化速率；Cxo为硝化作用中所需的溶解氧半饱和浓度；θx、θfx分别为硝化和反硝化的温度依赖系数。

活性磷酸盐子系统6：

(18)

其中：PrB为每克溶解态BOD中的磷含量；Pfz为磷系统中每克浮游植物磷的吸收量。

将以上微分方程及相关数学表达式(6)～(18)联立起来组成一个七元一阶的微分方程组，利用Euler法或四阶Runge-Kutta法[22-25]求解。然后，将生化反应项代入方程(5)，以方程(1)～(3)求得的水动力场为背景场，便可得到叶绿素a浓度的时空分布。

4 辽东湾海域浮游植物分布特征模拟

浮游植物的丰度分布用叶绿素a的浓度来表征，采用上述生态系统动力学模型对辽东湾海域的水质和叶绿素a浓度的分布状况进行模拟分析。

4.1 计算域及参数设置

计算域范围为整个渤海范围，整个模拟区域内由25 888个节点和47 197个三角单元组成，最小空间步长约为80 m，时间步长为30 min。地形取自自然资源部网站，插值得到研究海域的水深，岸界取自海域水深地形的平面图。计算中所需常量的取值见表1。

表1 计算中所需的常量取值Tab.1 Constant value required in calculation

4.2 模型验证与模拟结果

为了获取准确可靠的数值模型，需要利用观测值对模型进行验证。分别于2020年5月和8月，在辽东湾海域布设4个水文和14个生态监测站位，对辽东湾海域的水动力和生态状况进行了监测分析。水文观测站位置坐标见表2，站位位置见图1。

表2 潮位和海流观测站位位置Tab.2 Location of tide level and current observation stations

图2给出了监测站位H1、H2实测潮位值和计算值的比较，可以看出，监测站位处潮汐属于不规则半日潮，即一天24 h内出现两次涨潮和两次落潮过程，涨落潮历时大致相同，计算和实测潮位过程的高、低潮位及过程线均吻合良好。说明本研究中的水动力模型模拟的海域潮波运动与天然潮波运动基本相似，采用的边界条件合理，能够真实地重演和反映海域内潮波的传递和变形。监测站位v1、v2的流速和流向模拟值与实测值对比见图3、图4，两个监测点的流速和流向模拟过程线与实测值吻合良好，呈现出海区典型的NE-SW往复流分布特征，v1站位略显旋转流特征。综上可见，本研究中水动力模型的模拟结果具有较高的精确性，能够满足后期生态系统指标浓度值模拟计算的需要。

本图基于自然资源部标准地图服务网站GS(2019)3266号标准地图为底图制作，底图边界无修改。The figure is based on the standard map GS(2019)3266 in the standard map service website of Ministry of Natural Resources of the People’s Republic of China, with no modifications of the boundaries in the standard map.图1 监测站位位置图Fig.1 Location map of surveyed station

图2 H1和H2站位的潮位验证Fig.2 Tidal level verification of H1 and H2 stations

为进一步增强模型的可靠性，本研究中对生态系统动力学模型的生态指标浓度模拟结果进行了验证，并基于2020年5月和8月的生态指标监测值来进行分析。表3给出了5月和8月2个航次14个站位叶绿素a浓度监测值和模拟值的对比，可以看出，模型的模拟结果与实测值吻合良好。同时，为了完整表达本研究中生态系统动力学模型的模拟结果，以及与保守水质模型[25]相比，体现该模型的精确度与优势度，图5、图6分别给出了2020年5月3个生态指标的浓度变化过程及一个潮期内4个典型潮时的叶绿素a浓度场的平面分布状况。

图3 v1站位的流速验证Fig.3 Tidal current verification of v1 station

图4 v2站位的流速验证Fig.4 Tidal current verification of v2 station

表3 叶绿素a浓度实测值与模拟值对比Tab.3 Comparison of chlorophyll a concentration results between measured and simulated values

图5 辽东湾湾顶海区生态指标浓度模拟值随时间的变化Fig.5 Change in simulated values of ecological index concentrations with time at top sea area in Liaodong Bay

图6 辽东湾4个典型潮时的叶绿素a浓度模拟场分布Fig.6 Distribution of simulated field of chlorophyll a concentration at four typical tides in Liaodong Bay

模拟结果表明：叶绿素a的分布具有明显的时空分布特征，辽东湾靠近辽河入海口、葫芦岛和长兴岛附近海区的叶绿素a浓度明显高于其他海区；8月的平均浓度明显高于5月，二者的比值在1.6左右(表3)。叶绿素a浓度的模拟值在5月中旬达峰值，为13.4 μg/L(图5)，之后由于藻类呼吸和沉降过程，又引起浓度的显著下降；而在8月，随着气温和溶解氧的升高，藻类的量会快速累加，叶绿素a的浓度则达到一年内的另一个峰值，结果与历年观测的浮游植物丰度呈现双峰结构的分布趋势一致，与实际情况基本吻合。

5 讨论

5.1 生态指标浓度的时间分布特征

近年来，以叶绿素a为主要表征指标的水体中浮游植物分布密度的研究广泛开展。本研究在借鉴前人基于监测数据数理统计分析的基础上，综合考虑实际水体的流动性特征，真实展现可能影响到整个生态系统的物理、化学和生物变化过程，在水动力数值模型基础上引入各生态指标之间的耦合作用，建立了辽东湾海域的生态系统动力学模型。本模型的数值模拟结果表明，其指标浓度分布具有明显的时间分布特征。提取辽东湾湾顶海区叶绿素a浓度的时间序列曲线图发现，在5月中旬浓度达到峰值13.4 μg/L，到5月15日基本是平稳期，而后由于藻类呼吸和沉降过程，5月16日之后浓度出现了显著的下降趋势，5月26日到6月初浓度趋于稳定。8月叶绿素a又出现另外一个峰值，平均浓度明显高于5月，二者的比值在1.6左右，表明叶绿素a浓度即浮游植物分布存在明显的季节性分布特征。而无机氮浓度呈现出5月初到6月初持续递减的趋势；活性磷酸盐浓度在5月6日达到最高值，然后随时间变化至5月18日(12 d左右)时浓度大幅度下降，到5月23日(5 d左右)，活性磷酸盐的浓度开始缓慢降低，最后基本以极低的浓度值稳定在6月初。与2013年建立的水质模型[1，25]模拟结果相比较，本模型在整个模拟过程中，生态指标浓度量值变化明显，规律趋势分辨率更高，模型计算结果在详细反映实际的内部生态机理变化特征方面更有优势，准确性和可靠性更高。

5.2 生态指标浓度的空间分布特征

诸多的研究成果中，采用断面观测和大面观测的方式来反映水体中指标浓度的变化状况，其全面性尚待提高。数值模型的典型特征就是在重演历史时期有限数据的基础上，通过模型验证和调参过程，进而预测未来时间段整个空间平面或立体的数据信息，使结果分析更加具有说服力。为了较全面地展现生态系统动力学模型的研究结果，本研究中，在给出各生态指标时间序列曲线结果的同时，也给出各生态指标模拟结果在辽东湾海域整个空间的分布状况。叶绿素a浓度的模拟平面图表明，辽东湾湾顶海区的叶绿素a浓度明显偏高，而且经统计分析发现，高、低潮时的浓度高值区范围均大于涨急、落急时刻的浓度高值区，辽东湾辽河口入海区、葫芦岛和长兴岛附近海区的叶绿素a浓度明显高于其他海区，这也与历史上辽东湾夏季赤潮易发区分布位置相一致。分析原因，可能是由于夏季辽河口处于丰水期，排入河口入海区的营养盐比较丰富，使湾顶海区的藻类大量繁殖，而且叶绿素a浓度高值区的位置分布也与辽东湾湾顶海区及葫芦岛一带海区为产卵场与索饵场的渔场特征相吻合[26-27]，体现出辽东湾叶绿素a浓度明显的空间分布特征。

6 结论

1)辽东湾海区叶绿素a的浓度分布具有明显的时空分布与季节性分布特征，本研究中建立的浮游植物分布生态系统动力学模型，对生态指标浓度的变化规律性展现出较以往保守型模型具有更高的分辨率，模拟结果与历年观测的浮游植物丰度分布趋势的双峰结构特征相吻合，说明模型具有较高的精确性。

2)本模型模拟的生态指标浓度分布规律与辽东湾历年的赤潮易发区和渔场分布特征一致，说明模型中考虑的整个生态系统的物理、化学和生物变化过程，能够反映出模拟时间段海水中的客观自然变化规律，证明文中建立的模型在反映内部生态指标间相互影响机理方面表达更准确，具有灵敏性，说明模型具有较高的可靠性。

3)对辽东湾海域浮游植物分布模型的模拟结果，可为根据叶绿素a高值区的时空分布位置科学推测和跟踪渔业资源密集区，以及提高捕捞效率或选择增殖放流海区提供重要参考；叶绿素a模拟结果对与浮游植物密度相关性极强的浮游动物，如毛虾集群效应等研究，也具有预测和预警作用，进而也可为规避电厂取水口堵塞风险、保障产业的安全稳定运行提供技术支撑，说明该模型具有较好的实用性。