潘家口水库热分层期溶解氧垂向分布及影响因素
2022-07-01黄岁樑孙逸璇
高 锐，黄岁樑，孙逸璇
(1.南开大学 环境污染过程与基准教育部重点实验室,天津 300071；2.南开大学 天津市城市生态环境修复与污染防治重点实验室，天津 300071；3.南开大学 环境科学与工程学院 水环境数值模拟研究室，天津 300071)
1 研究背景
受气温变化影响，水深超过7 m的水库、湖泊等天然水体由于垂向上受热不均会产生热分层[1]，该现象在夏季尤为常见[2]。热分层期水体从上至下分为混合层(变温层)、温跃层(斜温层)和滞温层(均温层)[3-4]，其稳定性受水深、水体流动性以及气象条件等因素的协同影响[5-6]。
溶解氧(DO)在水生态系统结构与功能演化中起着重要作用。DO浓度会影响水体氮、磷的主要存在形态及其迁移转化过程，浓度低于6.5 mg/L时会影响鱼类的生长、繁殖等生理过程和空间分布[7-9]。已有研究表明，热分层期DO垂向分布与诸多因素有关，水温是主要因素之一，水体热分层会引发DO垂向产生混合层、氧跃层和滞氧层[5]。其中，氧跃层DO浓度变化剧烈，当热分层对DO影响显著时，层内会出现溶解氧极小值(metalimnetic oxygen minimum,MOM)[10]，MOM区域的缺氧环境(DO浓度<2 mg/L)会导致鱼类及浮游动物等窒息死亡，从而危害水生态环境[11]。此外，浮游生物、光照强度、有机物和pH等均会影响天然水体中DO的垂向分布[11-12]。
目前，国内外对湖、库热分层及其水质的垂向变化关注较多，如Beutel[13]研究了湖泊热分层期滞温层缺氧条件下沉积物中氨的释放过程。Edwards等[14]研究了Erie湖泊初级生产力与水体缺氧现象之间的关系。Noori等[15]研究了Karkheh水库热分层对溶解氧和营养盐的影响。田盼等[16]研究了神农溪不同季节DO与叶绿素a的垂向分布特征。然而，以往研究多集中于热分层对溶解氧、营养盐等水质指标的影响，对导致氧跃层MOM区域形成的原因研究较少。因此需加强对热分层期DO垂向分布、溶解氧极小值等缺氧区域演化过程等方面的研究，这对水资源保护与治理意义重大。
潘家口水库属于大型人工深水湖泊，位于河北省迁西县滦河干流上，是天津市和唐山市的重要水源[17]。为促进经济发展，水库初期大力发展水产养殖产业。随着网箱养鱼规模不断增加，水体中饵料及鱼类排泄物不断积累，致使库区水体富营养化，水质常年为V类或劣V类[18]。近年来虽已全面清理网箱养鱼，但水质仍处于中营养水平，氮、磷等营养盐浓度高，时常有藻类水华爆发[19]，水环境状况及演变趋势受到外界高度关注。本研究针对水库存在的缺氧等水质问题，对坝前水温、DO、浮游植物等指标进行监测，分析库区DO垂向变化特征，并探讨导致氧跃层MOM形成的原因，以期为未来水库生态环境保护与水资源利用提供一定的科学依据。
2 材料与方法
2.1 采样点设置
潘家口水库是引滦入津的主体工程，位于河北省滦河干流，总库容为29.3×108m3。其控制流域面积占滦河全流域的75%，为33 700 km2，多年平均径流量占比53%，为24.5×108m3，正常蓄水位为222.00 m，历年最大水深约为80 m，水库所在区域属温带大陆性季风气候，年平均气温为15 ℃[20-21]。潘家口水库水深库大，流速很小，夏季热分层期各物理量空间分布趋于稳定，为减小水面风力、船舶航行等因素对采样的影响和实现固定位置采样，本研究在距潘家口水库大坝前约1 km处(40°24′00.8″N，118°17′10.7″E)设置固定平台进行现场监测和采样。该点在采样期水深大约为60 m。采样时间为2019年8月，共采样4次，采样点位置如图1所示。
图1 潘家口水库监测采样点位置
2.2 水温、水深、pH、电导率、总溶解性固体(TDS)、叶绿素a(chl-a)、DO监测
水温、水深、pH、电导率、TDS、chl-a、DO利用多参数水质监测仪(YSI-EXO2，A Xylem brand1725 Brannum Lane Yellow Springs,USA)进行现场监测。YSI-EXO2被固定在采样平台上，通过控制转轴将其从表面缓慢投放至库底，再从库底缓慢拉回到水面，采样频率为1次/s。
2.3 总氮(TN)、总溶解性氮(TDN)、总磷(TP)、总溶解性磷(TDP)、高锰酸盐指数测定
利用采样平台固定5 L卡盖式采水器(NISKIN，General Oceanics，USA)采集水样。采样深度从水面起分别为1、3、6、9、12、15、20、25、30、40、50、57 m，共采集12个水样。水样加硫酸冷藏于保温箱中带回实验室分析，TN、TDN、TP、TDP浓度和高锰酸盐指数采用《水和废水监测分析方法》(第四版)[22]测定，藻细胞密度采用藻类实验室分析仪(BBE ALA-03-10)测定。
2.4 水温和溶解氧结构划分
温跃层和氧跃层垂向梯度划分方法在热分层水体中被广泛使用，其划分梯度值依地区而定，范围从0.05 ～ 1.00不等[23-26]。本文分别采用0.5 ℃/m和0.5 mg/(L·m)作为水温和溶解氧各层的划分界限，并采用各层上边界对应水深为各层水深。
2.5 浮力频率
浮力频率常用于表征流体垂直位移的稳定性。当流体质点在稳定环境中受浮力和重力作用做垂向运动时，产生的振荡频率被称为浮力频率[27]，公式如下：
(1)
式中：g为重力加速度，取9.8 m/s2；z为水深，m；ρ为水体密度，kg/m3。当N2<5×10-5s-2时，水体处于混合状态；当 5×10-5s-2
忽略水体泥沙和盐度的影响，水温对应的水体密度由公式(2)计算[28-29]：
ρ=1×10-11T6+5×10-9T5-1×10-6T4+1×10-4T3-
9.1×10-3T2+6.79×10-2T+999.84
(2)
式中：T为水温，℃。
2.6 数据分析
采用ArcGIS 10.7绘制研究区域监测位置。使用SPSS19.0对各项指标进行Pearson相关性及回归分析。利用Origin 2018和Excel对数据进行统计分析。
3 结果与分析
3.1 水体热分层稳定性
潘家口水库采样点最大水深为58 m。2019年8月潘家口水库热分层结构及稳定性特征见表1，水库水温及浮力频率(N2)的垂向分布曲线见图2。
表1 2019年8月潘家口水库热分层结构及稳定性特征
由表1和图2可知，水温呈现混合层、温跃层和滞温层结构。混合层水温垂向分布均匀，平均水温均大于24 ℃，层内N2范围为1.1×10-4～3.5×10-3s-2，分层较为稳定。温跃层水温沿水深方向显著降低，温度梯度最大可达到6.2 ℃/m(图2(c))，层内N2急剧增大，平均值大于5×10-4s-2，分层稳定。滞温层水温较低，N2平均值小于5×10-4s-2，水体处于弱分层状态。总体来看，从月初到月底，混合层平均水温由26.7 ℃降至24.5 ℃，厚度逐渐增大。温跃层水深由4.9 m增大至8.4 m，最大温度梯度由6.2 ℃/m降至4.4 ℃/m(图2)，N2呈波动减小趋势，稳定性逐渐减弱。滞温层水深由13.0 m增大至15.4 m，平均水温由7.7 ℃增大至8.4 ℃。
图2 2019年8月潘家口水库水温及浮力频率垂向分布曲线
3.2 溶解氧垂向分布特征
水体中DO的垂向分布特征如表2和图3所示。由表2和图3可见，监测期DO垂向分布呈“倒V”状，DO浓度混合层>滞氧层>氧跃层。混合层DO平均浓度均大于8 mg/L，浓度较高且垂向差异较小，厚度较小。氧跃层DO梯度显著增大，DO浓度先迅速减小后快速增大，DO平均梯度在0.89 ～ 1.72 mg/(L·m)之间(表2)，最大梯度可达到3.95 mg/(L·m)(图3(c))，水体局部区域处于缺氧状态(小于 2 mg/L)，并出现MOM。MOM在1.58 ～ 3.32 mg/L之间，最小MOM出现在8月10日水深5.1 m处(图3(a))。总体来看，从8月初到月底，氧跃层由1.6 m逐渐增大至6.8 m，其厚度由9.7 m逐渐减小至6.6 m，MOM从1.58 mg/L逐渐增大至2.40 mg/L，氧跃层逐渐下移，厚度呈现波动减小趋势且层内MOM逐渐增大(表2)。滞氧层DO逐渐恢复，厚度较大，DO平均浓度均大于5 mg/L，且沿水深方向缓慢减小。
图3 2019年8月潘家口水库溶解氧DO垂向分布曲线
表2 2019年8月潘家口水库溶解氧垂向分布特征
3.3 pH、电导率和总溶解性固体垂向分布特征
水体pH、电导率和总溶解性固体(TDS)垂向分布特征如图4所示。由图4可以看出，水体pH平均值为8.01，总体呈碱性。在垂向分布上，混合层和温跃层pH变化明显，其沿水深方向先减小后增大，且在温跃层内出现极小值，最小可达到7.63(图4(a))。滞温层pH分布逐渐均匀，且层内随水深增大pH略微减小。在时间分布上，从8月初至月底，pH平均值从7.85增大至8.28，极小值从7.63增大至7.92，极小值对应的水深也由5.1 m逐渐增大至10.3 m。
图4 2019年8月潘家口水库pH、电导率和总溶解性固体垂向分布曲线
电导率和TDS在垂向上存在明显分层现象，其变化趋势与pH相反。混合层电导率和TDS变化相对较小，层内电导率和TDS平均值分别为477.0 μS/cm和309.7 mg/L。电导率和TDS在温跃层中显著增大，电导率平均值由482.5 μS/cm增大至570.0 μS/cm，TDS平均值由312.7 mg/L增大至370.2 mg/L。滞温层电导率和TDS分布逐渐均匀，且沿水深方向略微增加。
3.4 氮、磷营养盐垂向分布特征
水库氮、磷营养盐垂向分布特征如图5所示。由图5可见，氮、磷营养盐沿水深方向整体呈增加趋势，TN和TDN浓度范围分别为4.3～6.3和4.1～6.2 mg/L，TP和TDP浓度范围分别为0.02～0.13和0.01～0.11 mg/L。混合层TN和TDN浓度分别在4.3～5.0和4.1～4.5 mg/L之间，氮浓度及其变化相对较小；TP和TDP平均浓度分别为0.025和0.01 mg/L，浓度较小且垂向分布相对均匀。温跃层氮、磷浓度显著增大，TN和TDN平均浓度梯度分别为0.056和0.034 mg/(L·m)；TP和TDP平均浓度梯度分别为0.003 2和0.002 9 mg/(L·m)，最大浓度梯度分别可达到0.009和0.011 mg/(L·m)(图5(d))。滞温层TN和TDN先波动增加，50 m之后逐渐减小，浓度分别在4.5～6.3 mg/L和4.5～6.2 mg/L之间，浓度较高且变化相对较小；TDP和TP在15～50 m变化幅度较温跃层低，50 m之后磷浓度明显增大，TP和TDP最大浓度分别为0.13和0.11 mg/L。
图5 2019年8月潘家口水库氮、磷营养盐垂向分布曲线
3.5 浮游植物垂向分布特征
浮游植物生物量垂向分布特征如图6所示。分析图6可知，混合层、温跃层和滞温层浮游植物平均生物量分别为2.30、0.89和0.61 mg/L。垂向分布上，绿藻(Chlorophyta)主要分布在混合层，生物量较大；硅藻(Bacillariophyta)主要分布在温跃层和滞温层，分布范围较广；隐藻(Cryptophyta)前期主要分布在混合层，后期在温跃层中生物量增大；蓝藻(Cyanobacteria)数量较少，在混合层有少量分布。时间分布上，8月10日浮游植物生物量较小，优势藻种为绿藻和隐藻，分别占生物总量的25%和52%。8月10日之后浮游植物繁殖速度增大，水深15 m以下硅藻大量繁殖，优势藻种转变为硅藻和绿藻，分别占生物总量的56%和37%。
图6 2019年8月潘家口水库浮游植物生物量垂向分布特征
4 讨 论
4.1 潘家口水库热分层结构特征分析
通过判定水体类型来判断水库分层结构在湖泊学中被广泛应用[30]。α判定值(年入库总水量/总库容)和滞留时间(库容/入库流量)是常见的水体类型判定指标，两者将水体分为河流型、过渡型和湖泊型[31]。通过计算可知，潘家口水库α判定值为0.45(α<10)，滞留时间为2.20年[32](大于 300 d)，符合湖泊型水库的特征。
潘家口水库混合层水体受太阳辐射和风力扰动等外界条件的影响，水温较高且垂向混合均匀。其下方水体受太阳辐射影响逐渐减小，水温沿水深方向温差逐渐增大，热分层逐渐形成。热分层阻碍了能量的垂向传递，导致水温梯度急剧增大而形成温跃层。温跃层较大的温度梯度进一步抑制了能量的垂向传递，致使温跃层下方水体由于长期能量摄入较少而处于低温状态(小于 10 ℃)，且水体由于所处位置较深受外界扰动较小而形成滞温层。即潘家口水库热分层呈混合层-温跃层-滞温层结构，这与Effler等[10]和Noori等[15]对水库热分层研究得出的结论一致。
图7(a)为水温和气温随时间的变化趋势，图7(a)显示，表层水温与日平均气温变化趋势一致，混合层与温跃层分界处水深较浅，受气温影响较大，水温较高，同时水体受外界扰动影响较大，该处水温波动幅度较大。温跃层由于存在较大温度梯度，能量的垂向传递受到限制，致使温跃层与滞温层分界处水温受气温影响较小，水温较低，同时分界处水体受外界扰动影响较小，该处水温波动幅度较小。图7(b)为表底温差(ΔT)与浮力频率(N2)的相关关系。由图7(b)可知，ΔT与N2存在显著线性关系，R2为0.95，且N2>5×10-4s-2，说明能量输入差异(ΔT)是导致潘家口水库热分层形成的主要因素，同时该差异也决定着水体热分层的稳定性[11,33]。余晓等[32]对深水水库的研究表明，表底温差不仅限制了水体的垂向掺混，致使热分层现象产生，而且还影响着热分层的季节性交替。
图7 2019年8月潘家口水库热分层期水温与气温变化趋势及表底温差与浮力频率的相关关系
4.2 潘家口水库溶解氧垂向分布特征分析
潘家口水库热分层形成后，DO垂向逐渐产生混合层-氧跃层-滞氧层结构，其垂向分层结构与热分层相似。图8为水库水温、浮力频率(N2)、溶解氧(DO)和浮游植物生物量的垂向分布曲线。由图8可以看出，各层水温与溶解氧的变化趋势基本同步。水温和DO混合层变化趋势一致，深度随时间逐渐下移，厚度分别增加5.2和3.5 m。氧跃层水深较温跃层浅，氧跃层厚度由9.7 m减小至6.8 m，温跃层厚度由8.1 m减小至7.0 m，两者随时间的变化趋势一致。滞氧层在8月初水深为11.3 m，后期随着热分层减弱，其水深逐渐下移至13.4 m。
图8 2019年8月潘家口水库水温、浮力频率、溶解氧和浮游植物垂向分布曲线
表3为潘家口水库监测点水体各项理化指标相关关系。表3中的相关性分析结果表明，DO与水温存在显著正相关关系，相关系数为0.310(P<0.01)，与N2呈显著负相关关系，相关系数为-0.338(P<0.01)，说明DO垂向分布与热分层密切相关。此外，当水库光补偿点水深较浅时，光补偿点下方水体由于耗氧未能得到及时补充，致使热分层期氧跃层深度较温跃层浅。邱晓鹏等[34]对枣庄周村水库的研究表明，温跃层和氧跃层位置受气温和光补偿点位置等条件的影响，氧跃层深度较温跃层浅，其厚度较温跃层大。综上可知，热分层是导致DO分层的主要原因，同时影响着DO层化结构的时空变化特征，这与Zhang等[25]的研究结果一致。
表3 2019年8月潘家口水库水体各项理化指标相关关系
水体各层DO浓度由耗氧和复氧过程共同决定。表层水体复氧主要来源于大气和浮游植物光合作用的补给，深水区复氧主要通过水体对流扩散补充。而耗氧主要来源于水体中有机物质氧化分解及浮游生物的呼吸作用等生物化学循环过程[35]。淡水中不同温度DO饱和浓度表达式[36]为：
(3)
式中：CDO为溶解氧浓度，mg/L；T为水温，℃。由理论模型可知，当水温沿水深方向减小时，DO浓度垂向变化随水深增加而逐渐增大。然而对于潘家口水库等深水型水库，溶解氧受热分层和不同生化过程协同作用，各层呈现出不同的演变特征。
混合层DO受大气复氧和浮游植物的光合作用产氧影响，且水库氮、磷营养盐充足，混合层水温适宜，浮游植物生物量高达2.30 mg/L，DO的饱和度为112%，处于过饱和状态，同时受风力扰动等外界条件影响较大，DO垂向分布较为均匀。
温跃层内上侧DO浓度显著降低，对应位置形成了6.6～9.7 m的氧跃层，并在该层出现了MOM，其主要原因有：
(1)温度梯度限制混合层DO对氧跃层的补给[37]。研究表明，水体温度梯度会显著影响DO垂向扩散通量，梯度越大，扩散越困难[6]。潘家口水库热分层形成后，温跃层水温显著降低，平均温度梯度高达2.26 ℃/m，浮力频率均大于5×10-4s-2，水体处于稳定分层状态，稳定的热分层阻碍了上层水体对该层DO的补给，导致温跃层DO显著降低。将溶解氧与水温、浮力频率垂向分布曲线(图8)进行对比也可以看出，MOM均出现在水温显著降低和浮力频率急剧增大的区域，这与官卓宇[4]对澜沧江中游水库溶解氧的研究结果一致。此外，表4给出了氧跃层水深(OD)和MOM水深(MD)与温跃层水深(TD)和温跃层厚度(TT)的线性拟合关系，由表4可知，OD和MD分别与TD和TT存在正相关和负相关关系，说明温跃层显著影响氧跃层的垂向分布和MOM的形成[34]。
表4 OD和MD与TD和TT的线性拟合关系
(2)温跃层有机物质氧化分解耗氧[38]。监测期间潘家口水库混合层、温跃层和滞温层高锰酸盐指数分别为2.94、2.80和2.71 mg/L，说明混合层和温跃层有机物浓度较高，水体化学需氧量大。水温是影响藻类生长的主要因素之一[39]，潘家口水库水温与浮游植物生物量呈显著正相关关系，相关系数为0.624(P<0.01)。水体混合层水温高，适宜藻类生长繁殖，藻死亡后形成的有机颗粒物较小，沉降速度也较小，而温跃层厚度较大(大于 7 m)，有机颗粒物有足够的时间在温跃层被微生物氧化分解，这一过程消耗氧气，并释放出CO2和有机酸，同时温跃层较大的温度梯度阻碍了上、下层溶解氧对该层的补给，致使温跃层DO浓度持续降低，pH显著降低，水体逐渐发展为缺氧状态，并出现MOM现象。相关性分析结果(表3)同样表明，DO与pH呈正相关关系，相关系数为0.741(P<0.01)，表明两者变化趋势一致。卡尔夫[40]研究也表明，进入温跃层的大部分有机颗粒体积较小且沉降速度很小，10 μm颗粒沉降10 m需要40 d，证实了上层有机颗粒物下沉经过温跃层时需要较长时间。Kreling等[35]对德国Arendsee湖MOM的形成研究发现，上层颗粒有机碳的下沉促进了MOM的发展。官卓宇等[41]的监测结果也表明，MOM的形成与温跃层中溶解氧的大量消耗和难以得到补给有关。此外，热分层期温跃层下方水体DO浓度较低，增加了低氧区域及沉积物中还原性气体的释放量，并影响了还原性气体的垂向扩散，导致温跃层耗氧进一步增大。Mcclure等[42]对Falling Creek水库研究发现，水体热分层改变了水体中甲烷和二氧化碳的垂向分布，导致甲烷在温跃层中积累矿化耗氧，并释放出二氧化碳，进一步促进了MOM的发展。
(3)浮游植物生命活动耗氧[35]。杨正健[30]在三峡库区的研究中发现，真光层内浮游植物光合作用大于呼吸作用，真光层以下的水体，浮游植物呼吸作用大于光合作用，浮游植物生长受限制，导致DO浓度降低。相关性分析结果(表3)表明，浮游植物生物量与TN和TP均呈显著负相关关系，相关系数分别为-0.624和-0.533(P<0.01)，与DO呈显著正相关关系，相关系数为0.537(P<0.01)。其原因是浮游植物生命活动消耗氮、磷营养盐并释放出氧气，水体溶解氧浓度增加。此外，图9给出了温跃层DO与叶绿素a、浮力频率的线性拟合关系，由图9(a)可知，温跃层DO浓度与叶绿素a浓度呈显著正相关(R2=0.44，P<0.05)，这进一步说明温跃层浮游植物光合作用大于呼吸作用，浮游植物呼吸作用耗氧不是MOM产生的主要原因。
图9 温跃层DO浓度与叶绿素a浓度及浮力频率的线性拟合关系
滞氧层与滞温层水深基本一致，DO沿水深方向缓慢降低。Müller等[43]研究发现，滞温层DO消耗存在两种途径：一是沉积物向底层水体中扩散的还原物质消耗氧气。二是沉积物表面的有机物消耗氧气。由图8可知，滞氧层上侧浮游植物生物量较大，以适宜在12 ～ 25 ℃生长的硅藻为主，DO浓度较高，靠近底泥附近DO浓度显著降低。根据潘家口水库监测数据表明，水库氮、磷营养盐浓度较高，这与长期以来水库的网箱养殖活动有关，网箱养殖产生的鱼食残渣沉降到沉积物表面，导致底泥有机物质污染严重，沉积物中有机物质的矿化作用消耗氧气并向水体释放氮、磷营养盐，层内DO浓度持续降低[32,44]。
热分层改变了水体垂向上的热量传递，致使DO垂向出现分层。由图9(b)可知，温跃层N2与DO呈显著负相关(R2=0.51，P<0.05)，其原因是水体热分层形成时，温跃层显著增大的N2限制了溶解氧的垂向补给，延长了耗氧生物和有机物质等在温跃层滞留的时间，使得温跃层氧气消耗加剧，这一系列过程为氧跃层和MOM的形成创造了条件，官卓宇[4]和俞焰等[36]的研究也证实了这一观点。综上所述，潘家口水库水体热分层是导致MOM形成的主要物理因素，有机物质的氧化分解及浮游植物的呼吸作用耗氧对温跃层DO的共同作用是引起MOM形成的生物、化学因素。
5 结 论
(1)潘家口水库8月份存在显著热分层现象，水温呈混合层-温跃层-滞温层结构。表、底层能量输入差异是促进潘家口水库水体热分层形成的主要因素。
(2)热分层期DO垂向呈混合层-氧跃层-滞氧层结构，分层结构与热分层相似。相关性分析表明，DO与水温呈显著正相关，与表征热分层稳定性的N2呈显著负相关，说明热分层不仅是DO垂向分层的主要原因，而且还影响着DO层化结构的时空变化特征。
(3)溶解氧各层呈现不同的变化趋势。混合层溶解氧受大气复氧和浮游植物光合作用产氧影响，溶解氧处于过饱和状态，垂向分布相对均匀。滞氧层受沉积物中有机物等还原物质耗氧影响，DO沿水深方向逐渐降低。受温度梯度影响，混合层和滞氧层对氧跃层溶解氧补给减小，同时氧跃层有机物质氧化分解耗氧增加，致使氧跃层DO显著降低，并出现MOM；MOM水深与热分层参数(温跃层水深和厚度)存在相关关系，DO与pH、浮游植物生物量、氮、磷营养盐呈显著相关关系，说明热分层是MOM形成的主要物理因素，水体有机物质的氧化分解及浮游植物呼吸作用耗氧是导致MOM形成的生物、化学因素。
