运用AEM3D模型分析滇池蓝藻对营养盐变化的响应

杜成栋,潘继征,张国正,李 勇,罗纯良,李清濯,华跃洲 (1.苏州科技大学环境科学与工程学院,江苏 苏州 21011;2.江苏省环境科学与工程重点实验室,江苏 苏州 21011;3.中国科学院南京地理与湖泊研究所湖泊与环境国家重点实验室,江苏 南京 210008;4.江苏纬信工程咨询有限公司,江苏 南京 210008;.昆明市滇池高原湖泊研究院,云南 昆明 60000;6.中国科学院大学,北京 100049)

随着人类活动加剧,点源、面源等外源污染物进入湖体[1],加上自身内源污染等一系列因素,使得湖体水质受到不同程度的影响[2],其中一部分水体出现富营养化现象。在适宜的条件下,水体富营养化通常引起某些特征性藻类迅速繁殖,最终导致蓝藻暴发[3]。根据《2020中国生态环境状况公报》,我国的110个重要湖泊(水库)中,中度富营养化状态占61.8%,重富营养化状态占29%,其中滇池为重度富营养化[4]。

关于蓝藻形成机理的研究,营养盐与藻类生长之间的关系最受学者的关注。氮和磷是藻类水华暴发的关键营养元素[5-6]。已有大量研究表明,磷是淡水生态系统中藻类暴发的主要限制性营养物质[7]。如SCHINDLER[8]对加拿大安大略省西北部湖进行氮和磷的添加试验,发现添加磷的试验几周后出现了蓝藻暴发的现象,表明磷有利于蓝藻生长。氮被认为是海洋生态系统中的限制因子[9]。学者们对不同湖泊数据进行收集和整理[10],发现蓝藻的暴发受氮和磷的共同作用,而不是单一因素发挥作用。从时间角度来讲,在春季磷通常是限制因子[11],夏季和秋季的温度和气象条件有利于蓝藻生长,氮起到了关键作用[12]。国内学者通过对太湖梅梁湾的水质监测和原位营养盐添加实验发现,春季和冬季浮游生物的数量及生长率随着磷含量的增加而增加,到了夏季和秋季氮则替代磷成为主要的限制因子[13]。因此蓝藻的暴发与水体中氮磷营养盐的含量密切相关,通过模型模拟营养盐与蓝藻的变化关系可以为防治蓝藻暴发提供依据。常规的一维模型难以模拟浅水湖泊在垂直方向上水体中营养盐的分层现象,如SALMO(Simulation of an Analytial Lake Model)模型等;二维模型则缺乏构建湖泊的水动力特征,如CE-QUAL-W2(Water Quality and Hydrodynamic in 2D)等;三维模型则可以很好地模拟湖泊流体动力学。

先前有关滇池蓝藻与营养盐的关系报道多以整个滇池(含草海和外海)或滇池外海作为一个整体区域进行研究[14-15]。而滇池外海作为滇池的主体水域,受气象和水文条件的影响,其水质和蓝藻发生规律均具有较大的空间差异性[16]。因此,有必要对滇池外海蓝藻进行分区研究。笔者结合最新的滇池蓝藻水华程度级别界定标准,结合外海水质监测点位布设规律,将外海划分成南、中和北部3个大区,每个大区中又分成8个小区,总共24个分区(图3),探讨滇池蓝藻与营养盐的时空分布特征及营养盐调控对蓝藻的影响。同时通过滇池水质与气象数据构建三维水动力水质模型,采用人工和蒙特卡洛自动率定相结合的方法确定模型主要参数,采用2018年的数据对模型进行率定,再应用2019年的实测数据对模型模拟结果进行验证,以检验模型参数的适用性。使用率定及验证校准后的模型,设置不同的外源氮、磷削减情景(仅改变氮、磷或氮磷共同改变),探究外源营养盐削减与滇池湖体营养盐及蓝藻对应变化的响应关系,为外源营养盐调控及减少蓝藻暴发提供科学支撑。

1 研究区概况与水质数据

1.1 研究区域

滇池流域位于云贵高原中部,属长江流域金沙江水系[17-18],是昆明市主要的社会经济发展区域,滇池流域面积为2 920 km2,占昆明市域面积的13.8%[19],南北长42 km,东西宽12.5 km,平均宽8.32 km,湖岸线长163 km,平均水深4.4 m,最深处(海眼)11.2 m,总蓄水为15.6亿m3[20]。自20世纪80年代以来,滇池水体富营养化严重,影响了周围群众的生产活动[21]。

1.2 数据来源

气象、地形、水文和水质数据来源见表1。

表1 数据来源Table 1 Data sources

根据昆明市高原湖泊研究院提供的10个滇池国控监测点1990—2019年逐月水质监测数据(Chl-a浓度数据从1999年开始)将滇池分为4个区域,具体点位分布见图1。除上述数据外,2019年每季度对滇池湖体进行一次采样。主要的分析指标有:溶解氧(DO)浓度、pH值、透明度(SD)、水温(WT)、总氮(TN)浓度、总磷(TP)浓度、Chl-a浓度。其中SD采用塞式透明度盘测量;DO浓度、pH值和WT采用携式多参数水质分析仪(AP 2000,Aquaread 公司,英国)测量;TN浓度采用碱性过硫酸钾消解-紫外分光光度法测定[22];TP浓度采用钼酸铵分光光度法测定[22];Chl-a浓度采用热乙醇提取法测定[23]。数据采用统计学软件SPSS 25.0和MATLAB 2012a分析处理,图表采用Origin 2017绘制。

图1 人工采样点分布Fig.1 Distribution of sampling sites

2 模型与研究方法

2.1 AEM3D模型介绍

AEM3D模型是在ELCOM-CAEDYM模型的基础上改进而来的,其中ELCOM模型是三维水动力模型,CAEDYM模型是水质模型,两者耦合能模拟受自然因素或人为活动影响的水体中的温度、流速等水动力变化,也能模拟物理、化学、生物综合影响的水质变化。

2.1.1水动力模块

ELCOM是由西澳大学水研究中心研发的三维数值模型,运用水动力学和热力学原理,模拟河口、湖泊与沿海水域的不同深度分层水体的临时状态,可用来预测受到的风应力、表面加热与冷却等环境力作用时水体的速度、温度和盐度的分布。ELCOM模型是在三维笛卡尔网格的计算单元上面实现模拟的。

主要控制方程组如下:

动能传输公式为

(1)

式(1)中,Uα、Uβ为湍流动能,m2·s-2;Uj为浮力;ρ为水体密度,kg·m-3;v1、v2、v3分别为x、y、z方向上湍流系数。

传输标量为

(2)

式(2)中,C为常数;Sc为传输标量常数;k1、k2、k3分别为x、y、z方向上相关系数。

2.1.2水质模块

CAEDYM模型是一个基于过程的水质模型,与ELCOM耦合时,在每个ELCOM设定的时间步长之后更新生化过程。CAEDYM模型的生态配置可以选择不同的水质模块,如无机物模块、沉积和悬浮模块、物质循环(碳、氮、磷、硅)模块、浮游植物模块、浮游动物模块、高等生物(鱼类、底栖类和水母类)模块等。CAEDYM 模型比传统的N-P-Z模型更先进,因为它是一个通用的生物地球化学模型,可以模拟物种或群体特有的生态相互作用。它独特的盐度算法能够模拟各种水生环境。

CAEDYM模型营养循环模拟过程中的常用参数如下文所示。

温度:

(3)

式(3)中,T为藻类生长温度,℃;g为藻类生物组数;kg为常数;νg为藻类生长温度系数;ag为实测温度,℃;bg为藻类生长温度补偿系数。

溶解氧:

(4)

式(4)中,ρDO为溶解氧浓度,mg·L-1;KDOg为藻类分解所需溶解氧的半饱和常数。

沉积物通量释放:

(5)

式(5)中,Sg为藻类生物量浓度,L-1;KDOS-g为沉积物需氧量的半饱和常数;KpHS-g为沉积物酸碱度的半饱和常数;CpH为水体的酸碱度;Zbot为底层深度,m。

再悬浮:

(6)

式(6)中,gsed为沉积层藻类生物量;αg为再悬浮速率常数;τ为切应力;τcg为蓝藻再悬浮临界切应力;Ksed-g为藻类再悬浮控制常数。

有机物分解:

(7)

CAEDYM模型中N、P模块的循环结构如图2所示,分别代表了矿化、生物摄入、分解、沉降、再悬浮等过程。

DONL—活性溶解性有机氮;DONR—难降解溶解性有机磷;PONL—活性溶解性有机磷;DONR—难降解溶解性有机氮;PIN—颗粒态无机氮;DON—溶解性有机氮;PON—颗粒有机氮;IN—无机氮;PIP—颗粒态无机磷;FRP—可溶性活性磷;POP—颗粒态有机磷;IP—无机磷。

应用2018—2019年滇池湖体及出入湖河流的水质数据、昆明站监测气象数据等,以AEM3D模型为基础,构建滇池三维水动力水质模型,采用人工和蒙特卡洛方法对模型主要参数进行校正,基于2018年的实测数据与模拟数据进行对比,找到误差相对较小的一组参数,随后应用该组参数进行2019年的水质模拟,应用2019年的实测数据检验模型的模拟结果是否同样适用。

滇池湖体的地形数据来自云南省生态环境厅,对该数据进行提取后,按照深度将湖体划分成6层,表层深度为1 m,最后一层深度为2.5 m。将滇池分成200 m×200 m的网格,为166×103个,其中有效湖体网格为7 306个。对滇池外海水域进行分区,总体分成南(S)、中(M)和北(N)3个大区,每个大区又分成8个小区,总共24个分区,其命名规则、分区信息和出入流河道位置如图3所示。

图3 滇池模拟分区图Fig.3 Simulated zoning map of Dianchi Lake

汇入滇池的主要河流有35条,其中部分为季节性断流,因此此次模型使用其中31条非季节性断流河流。为了配置模型的入流文件,对上述河流的水位、流量和水质情况进行收集,入流河道断面监测点的经纬度见表2。

表2 模型主要入湖河道监测点经纬度Table 2 The longitude and latitude of the main monitoring points of the river entering the lake

气象文件作为AEM3D模型的边界条件,是输入文件最重要的一部分,该文件包含处理后的时间、气压(AP)、气温(AT)、相对湿度(RH)、降雨(PP)、风速(WV)、风向(WD)等参数数据。其中大部分数据由昆明市滇池高原湖泊研究院提供,少数自测,由于悬浮物浓度(SS)目前监测数据较少,模型中暂定为10,盐度默认为0。

2.2 模型率定

将模拟结果与实际观测数据对比是模型必不可少的评估步骤[24]。对于复杂的水质模型,其生物化学参数众多,若采用传统的人工率定方法来确定水质参数,不仅费时费力,还依赖模型使用者的经验和技术,因此采用蒙特卡洛方法与人工调参相结合。蒙特卡洛方法是一种以概率统计理论为主要理论基础的随机抽样计算方法,优点在于避开了描述研究对象的复杂特性及其对整个系统的影响,直接以运用过程模拟替代系统分析[25]。

研究采用蒙特卡洛随机抽样方法对模型参数进行率定,首先根据文献[26]设定参数的取值范围,每次计算时在取值范围内随机生成一组参数值,运算结束后,比较模拟值与实测值之间的误差,多次迭代模拟后,选择模拟值与实测值误差最小的一组参数值作为模型的最优化参数。最终确认的参数见表3。

表3 AEM3D模型参数取值Table 3 Parameters of AEM3D model

采用均方根误差(RMSE,ERMS)、纳什系数(NSE,ENS)和相对误差(RE,ER)作为模型评价指标,表达式为

(8)

(9)

(10)

2.3 外源营养盐削减情景设置

基于AEM3D模型,考虑到外源营养盐与湖体的相应关系,设置不同程度的营养盐削减比例,探究湖体的响应关系。具体设计方案如下:削减外源TN浓度(10%、25%、50%、75%)、削减外源TP浓度(10%、25%、50%、75%)、同时削减外源TN与TP浓度(10%、25%、50%、75%),与2019年实测值进行对比。对削减情景模拟结果开展分析,探究不同的削减策略对滇池营养盐和蓝藻含量与分布的影响。

3 结果与分析

3.1 率定与验证

滇池水体中TP率定和验证的结果如图4所示,实测值和模拟值存在较高的吻合度。率定期外海TP的模拟效果NSE达到0.93,全湖TP的RMSE为0.005,NSE为0.86,RE为5%,验证期全湖TP的RMSE为0.01,NSE为0.67,RE为11%。率定期和验证期TP的RMSE均在0.011以内,NSE均在0.6以上,RE均在12%以内,TP模型模拟值和观测值之间拟合程度良好,能正确反映滇池湖体的TP浓度变化情况。

图4 滇池水体中TP浓度的模拟与验证对比Fig.4 Comparison of TP simulation and verification concentrations in Dianchi Lake

滇池水体中TN率定和验证的结果如图5所示,TN实测值和模拟值存在较高的吻合度。滇池草海和外海在率定期和验证期的TN模拟值与观测值基本一致,湖体TN模拟的RE在7.3%～19.6%之间。率定期全湖TP的RMSE为0.20,NSE为0.84,RE为8.4%,验证期全湖水温RMSE为0.34,NSE为0.68,RE为13.5%。率定期和验证期TP的NSE均在0.63以上,RE均在20%以内,所以建立的模型能正确反映滇池湖体的TN浓度变化情况,TN模拟值和观测值之间拟合程度良好。

图5 滇池水体中TN浓度的模拟与验证对比Fig.5 Comparison of TN simulation and verification concentrations in Dianchi Lake

Chl-a为湖体藻类浓度最直观的反映指标,滇池Chl-a率定和验证的结果如图6所示,Chl-a实测值和模拟值大小基本一致。湖体Chl-a模拟的率定期RE为13.7%,验证期RE为7.2%,RE均在15%以内;率定期RMSE为0.009,验证期RMSE为0.007,RMSE均小于0.01;率定期NSE为0.90,验证期NSE为0.87,NSE均大于0.6。总体来看,模拟结果能够反映滇池Chl-a的整体变化规律。

图6 滇池水体中Chl-a浓度的模拟与验证对比Fig.6 Comparison of Chl-a simulation and verification concentrations in Dianchi Lake

3.2 营养盐与蓝藻分布特征模拟

图7为2019年滇池TN浓度时空分布模拟图,从图中可以看出,草海的TN浓度较高,随时间变化显著(P<0.05),月均最大值为4.09 mg·L-1,出现在9月的断桥点位,月均最小值为0.72 mg·L-1,出现在6月的草海中心点位。外海的TN浓度较低,且随时间变化不明显(P>0.05)。7月全湖的TN平均浓度最高,为1.67 mg·L-1,而相邻的6月和8月出现全湖氮浓度降低现象,到了秋季全湖TN维持在较高的浓度。

图7 2019年滇池TN浓度时空分布模拟Fig.7 Spatio-temporal distribution simulation of TN concentration in Dianchi Lake in 2019

滇池TP浓度的时空分布模拟如图8所示,TP浓度变化范围为0.05～0.09 mg·L-1,平均值为0.07 mg·L-1。外海TP浓度的月变化相对显著(P<0.05)。草海夏、冬两季的TP浓度较高,模拟结果与该年的实际情况一致。外海的TP浓度从2月开始逐渐升高,到5月浓度达到峰值,夏季的TP维持在较高的浓度。

图8 2019年滇池TP浓度时空分布模拟Fig.8 Spatio-temporal distribution simulation of TP concentration in Dianchi Lake in 2019

图9为2019年滇池蓝藻Chl-a时空分布模拟图。从图9可以看出,一年之中滇Chl-a浓度有着明显的逐月变化规律,春秋季节浓度高,冬季浓度低,从冬季到春季Chl-a浓度逐渐上升。模拟结果显示,Chl-a浓度在5月和9月出现了2次峰值。

图9 2019年滇池Chl-a浓度时空分布模拟Fig.9 Spatio-temporal distribution simulation of Chl-a concentration in Dianchi Lake in 2019

3.3 外源营养盐削减情景模拟

3.3.1单因子营养盐削减

如图10所示,当外源TN削减10%时,湖体TN浓度呈现下降的趋势,平均浓度下降4.3%,同时湖体内部TP也受微弱影响,下降1.3%,Chl-a的年均值几乎无变化;外源TN削减25%时,湖体TN的浓度快速下降,年平均浓度下降11.5%,Chl-a年均值下降了3.4%;外源TN削减50%时,内源TP的变化依旧不明显,TN下降的趋势减弱,下降了16.3%;当外源TN下降75%时,湖体TN的降幅达到31.7%,TP的降幅为4.5%,Chl-a的年均降幅为8.5%。

图10 不同外源TN、TP单独削减下湖体TN、TP浓度及表层Chl-a浓度变化Fig.10 Changes in the corresponding concentrations of TN, TP and surface Chl-a concentration in the lake under the reduction of different exogenous TN and TP

如图10所示,当外源TP削减10%时,湖体TP浓度下降6.8%,同时湖体内部TN受其影响,下降1.5%,Chl-a年均值变化不明显;当外源TP分别削减25%、50%、75%时,内源TP分别下降12.7%、24.3%和34.9%,呈现一定的线性关系,内源TN全程无明显变化,Chl-a年均值从TP削减10%时下降0.9%迅速变为在TP削减75%时下降13.1%。

从图11和图12可以看出,在外源氮、磷削减的前期,Chl-a浓度无明显变化;随着时间的推移,Chl-a年内波动变大;在第120～150天时,Chl-a随外源氮、磷的削减量有明显的变化,Chl-a的峰值无明显下降趋势。春季外源TP削减25%与50%时对Chl-a浓度的影响相似,夏季外源TP削减10%的情况下Chl-a浓度反而比无削减处理高。

图11 不同外源TN削减下湖体表层Chl-a浓度随时间的变化Fig.11 Changes of surface Chl-a concentration with time under different exogenous TN reduction

图12 不同外源TP削减下湖体表层Chl-a浓度随时间的变化Fig.12 Changes of surface Chl-a concentration with time under different exogenous TP reduction

3.3.2氮磷营养盐同时削减

图13为不同外源TN、TP同时削减下湖体TN、TP及表层Chl-a变化,图14为不同外源TN、TP同时削减下湖体表层Chl-a浓度随时间的变化。从模拟结果看,当外源TN和TP输入量同时减少10%时,Chl-a下降2.7%,明显高于单一营养盐的削减效果;当两者同时削减25%时,Chl-a的变化和只削减磷处理的变化几乎一致;当两者削减量大于50%时,Chl-a变化高于单一营养盐削减。Chl-a在春季变化较为显著(P<0.05),在夏秋季节变化不显著(P>0.05)。

图13 不同外源TN、TP同时削减下湖体TN、TP及表层Chl-a浓度变化Fig.13 Changes of corresponding concentrations of TN, TP and surface Chl-a in the lake under simultaneous reduction of different exogenous TN and TP

图14 不同外源TN、TP同时削减下湖体表层Chl-a浓度随时间的变化Fig.14 Changes of surface Chl-a concentration with time under simultaneous reduction of different exogenous TN and TP

4 讨论

4.1 营养盐变化对Chl-a浓度的影响

已有研究表明,湖泊水体中的氮、磷营养盐含量是蓝藻生长的重要限制因子[28],这与该研究模型模拟结果一致(图11)。在沉积物中越冬的蓝藻是第2年春季蓝藻暴发的根源,随着温度的上升、光照的加强,冬季埋藏在沉积物中的蓝藻细胞开始进入到水体中,吸收水体中的氮、磷营养盐,并且快速繁殖生长[29],因此5月Chl-a浓度出现了峰值(图11)。

到了夏季光照强度大、温度高,同时大量降水对湖水中营养盐含量有稀释作用,此时藻类的生长受到抑制。国内学者对洱海藻类与环境因子关系的研究也发现,夏季洱海流域藻类生长被抑制,这种现象是降雨产生的地表径流对洱海水体中的营养盐浓度有一定程度的稀释作用导致的[30]。秋季光照开始减弱,且降水减少,更适合藻类的生长繁殖,Chl-a浓度达到最高值(图11)。冬季的水温为一年中最低,浮游藻类的生长受到抑制或衰亡,导致Chl-a浓度降低。研究表明,藻类本身含大量营养盐[31],由于冬季水环境的变化部分蓝藻会衰亡分解,在其分解过程中部分营养盐会释放到水体[32],导致Chl-a浓度低而氮、磷营养盐浓度较高。

蓝藻的暴发往往是生态系统突变的结果,Chl-a浓度的峰值相对均值更能评估蓝藻风险。根据研究结果显示,滇池Chl-a浓度在5月和9月出现了2次峰值,因此可以在该段时间之前采取调水、打捞、外源削减等综合措施减少蓝藻爆发的面积和次数。

4.2 外源营养盐削减对Chl-a浓度的影响

单因子营养盐削减结果显示(图10),当外源TN削减10%时,Chl-a浓度的年均值无明显变化,当外源TN分别削减25%、50%和75%时,Chl-a浓度年均值分别下降了3.35%、6.49%、8.50%;当外源TP削减10%时,Chl-a浓度年均值变化不明显;当外源TP分别削减25%、50%和75%时,Chl-a浓度年均值分别下降了5.53%、8.97%和13.01%,说明Chl-a浓度对外源磷的响应高于外源氮。马健荣等[32]在富营养化湖泊中发现,外源磷的增加通常会导致水体中Chl-a浓度在浮游植物群落中占优势,这与该研究结果其相似。外源氮、磷削减的前期,Chl-a浓度无明显变化,但随着时间的推移,Chl-a浓度波动变大,在第120～150天,Chl-a浓度随外源氮、磷的削减量有明显的变化。MINAUDO等对卢瓦尔河30 a的研究也发现,减少氮、磷营养盐的输入是降低Chl-a的可行方法之一[33]。春季外源TP削减25%与削减50%对Chl-a浓度的影响相似,夏季外源TP削减10%时Chl-a浓度反而比无削减处理高。考虑不同削减策略对藻类季节浓度的影响,春季和冬季磷削减引起的Chl-a浓度变化比氮削减显著。

从同时削减TN和TP的模拟结果看(图13),当外源TN和TP输入量同时减少10%时,Chl-a下降2.7%,明显高于单一营养盐的削减效果(图15);当两者同时削减25%时,Chl-a的变化与只减少磷的变化几乎一致,表明此时磷对Chl-a起到限制作用;当两者削减量大于50%时,Chl-a的变化高于单一营养盐削减(图15)。

图15 不同削减方式下湖体表层Chl-a浓度变化Fig.15 Changes of surface Chl-a concentration under different reduction methods

Chl-a春季变化较大,夏、秋季的变化不明显,可能是未考虑夏季调水对其的影响导致的。外源营养盐控制是治理湖泊蓝藻暴发的根本,否则任何湖内治理都很难获得长期效果。已有研究表明,同时控制氮、磷营养盐的输入来抑制蓝藻暴发具有很好的效果[34]。对德国康斯坦斯湖的研究发现,尽管温度适宜藻类的生长,但降低外源营养物质输入,蓝藻生物量仍得到大幅降低[35]。针对荷兰Veluwe湖[36]和德国穆格尔湖[37]的研究同样表明,降低外源营养盐输入量可以有效削减蓝藻的优势度。CULLEN等[38]评估了削减外源营养盐对43个湖泊的响应结果,发现当氮、磷削减量在50%～75%之间,55%的湖泊富营养化等级有明显的降低,营养盐与Chl-a浓度显著下降,其余湖泊营养盐与Chl-a浓度有小幅度下降,但不显著。UTTORMARK等[39]对现有的13个湖泊进行评估,发现当外源输入的TP与TN得到控制时,湖泊富营养化状态发生了明显的转变,Chl-a浓度大幅度降低,水质有所提升。华盛顿湖在外源营养盐削减88%时,水体富营养化与水质取得了明显的改善,蓝藻暴发现象减少[40]。因此通过调控外源营养盐减少蓝藻暴发是可行的手段。

5 结论

(1)基于滇池气象、水质及Chl-a浓度数据,利用AEM3D模型构建三维水动力水质模型,模拟结果能够较好的反应滇池实际水质情况。

(2)模拟结果显示,春季Chl-a浓度在不断增长,表明磷在藻类生长初期充当重要的营养盐因子,到了秋季,随着温度和气象条件有利于蓝藻的生长时,氮起到了关键作用。

(3)模拟单因子外源营养盐不同比例削减时发现,外源TP对蓝藻的影响大于外源TN。在仅控制TP时,削减25%较为合理;在仅控制TN时,50%的氮削减能有效地减少湖内Chl-a浓度。

(4)与单一营养盐的削减相比,氮、磷营养盐同时削减对Chl-a的影响更为显著,在同时削减10%或50%以上时,削减效果明显高于单一营养盐。但在同时削减25%时,Chl-a浓度的变化与只减少25%磷的变化几乎一致。