千岛湖极端水位变化对温跃层的影响

龚发露，王裕成，兰 佳，罗纯良，罗潋葱，李慧赟，张如枫

(1.云南大学生态与环境学院，云南 昆明 650500； 2.杭州市生态环境局淳安分局，浙江 杭州 311700；3.昆明市滇池高原湖泊研究院，云南 昆明 650500； 4.中国科学院南京地理与湖泊研究所，江苏 南京 210008)

温跃层是深水湖泊中常见的物理现象，其变化规律和分布特征极大地影响着湖泊生态系统的结构和功能。温跃层会阻碍表层与底层水体之间的物质和能量交换，降低底层水体的氧含量，导致底泥污染物释放，致使水体营养盐浓度升高。温跃层的厚度、深度与持续时间受湖泊形态特征、湖体水质和气象因子等因素的影响，如湖泊深度和表面积大小是决定湖泊是否有条件形成温跃层的重要因素[1]。在湖泊热力分层稳定期间，温跃层深度与表层水温存在显著负相关关系，与透明度和水深存在显著正相关关系[2]。气温升高可导致湖泊温度分层时间延长，湖泊稳定性增强。水位作为湖泊最基本的水文要素之一，对温跃层的结构和持续时间具有重要影响。水位下降后，中下层水体可接收更多的太阳辐射而升温，上下水层温差减小，分层强度减弱[3]。因此，如果人为大量排水致使水位下降，会在短期内迅速降低水体稳定度，可能使水体提前进入混合期[4]。强降雨和入湖流量的异常增多，可导致水位迅速上涨并使稳定的水体分层结构遭到破坏[5]。夏季水位的大幅波动会对温跃层结构产生显著影响[6]。在水位升降对温跃层的影响方面已有很多研究成果，但鲜有涉及极端水位(年内最高水位与最低水位之差)变化条件下温跃层特征变化的研究。受气候变化影响，从长时间尺度分析，我国长江中下游地区暴雨发生次数呈一定的增加趋势，相应会造成湖泊水位波动幅度和频率的变化[7]。因此，研究极端水位变化条件下的温跃层结构特征，对水库水环境综合管理、水利工程调度和湖泊生态系统稳定维持具有重要意义。

水位的监测手段主要包括水位原位监测、卫星测高、遥感提取、模型模拟[8-10]等，其中水文水动力模型和数理统计模型等模型模拟方法运用较为普遍。例如：李育等[11]采用瞬态气候演变模型和湖泊水量平衡模型对东亚及中亚地区全新世以来湖泊有效水分变化进行了模拟，分析了不同区域湖泊水位的主要影响因素；李琼等[9]基于水位-流量转换模型实现了宜昌市沮璋河洪峰水位的实时预报；华祖林等[12]从河流形态、水生生物和水质等多维度出发，计算了秦淮河不同河段的生态水位。这些研究采用的方法为湖泊水位计算提供了良好的借鉴。水动力学模型DYRESM(dynamic reservoir simulation model)具有参数较少且率定较易的优点，在湖泊水位研究中得到广泛运用。例如：Fadel等[13]通过DYRESM精确地模拟了Karaoun水库的水位波动情况，模拟值与观测值之间的误差仅为0.4～0.7 m；在缺少完整入流数据时，Weinberger等[14]基于水量平衡法，通过DYRESM计算出了Ammersee湖1975—2007年的逐日入流量；Gu等[15]利用DYRESM模拟了海南省松涛水库2016年3月至2017年3月的水位与水温，发现其能很好地反映水库的水位变化与温度分层情况。此外，DYRESM还被广泛应用于水体温度垂直分层的模拟，如加拿大Clearwater湖[16]、美国Mendota湖[17]和我国太湖[18]等。

目前对千岛湖温度分层及其影响因素的研究较多[19-21]，但主要集中于气象强迫因子(气温与降水)和湖体特征(深度与透明度)对温跃层的影响，很少考虑极端水位变化情况下的温度分层规律。因此，本文以研究千岛湖温跃层对极端水位变化的响应为目的，采用DYRESM模拟不同水位变化情景下的水温分层情况，并通过湖泊分析程序Lake Analyzer(LA)计算并分析温跃层特征量和湖泊稳定性参数在不同水位变幅情况下的变化规律，以期为千岛湖水环境保护提供参考。

1 研究区概况

千岛湖(29°22′N～29°50′N，118°34′E～119°15′E)位于浙江省杭州市淳安县境内，是拦蓄新安江而形成的人工深水湖泊。自1960年完成蓄水至2016年，千岛湖年均水位96.67 m(黄海基面)，年平均水位变幅为9.91 m，其中2007年水位变幅最小(4.29 m)，1973年最大(22.57 m)。2000年以来年平均水位变幅为7.01 m，最大水位变幅出现在2014年，为11.49 m[22]。多年平均出入湖水量分别为91.7亿m3与94.5亿m3。千岛湖入湖河道共有30多条，其中新安江、东源港、武强溪、枫林港和云源港5条为主要入湖河道(图1)，共占千岛湖总入湖流量的83.57%[23]，千岛湖主要出流为下游的千岛湖大坝。流域内降水集中于4—8月，多年平均降水量为1 733 mm，年均气温为17.3℃[24]，年平均日照时数为1 951 h，年蒸发量为1 362 mm[22]。

图1 千岛湖流域及湖体水质和气象监测点位Fig.1 Map of Lake Qiandao catchment and sites of water quality monitoring and weather station

2 研究数据与方法

2.1 数据来源

研究所用数据来源及其属性见表1，主要包括千岛湖出入湖河道流量、水位和气象数据。

表1 研究数据

2.2 研究方法

2.2.1 水动力学模型DYRESM

DYRESM是由西澳大利亚大学研发的一维水动力模型，可模拟湖库中温度、盐度及密度的垂向分布，其输入数据包括湖底高程、出入湖河道水量和气象数据，输出参数为不同深度水层的温度和盐度(适用于海洋及咸水湖)，可与水生态动力学模型(computational aquatic ecosystem dynamics model，CAEDYM)耦合，进行水质和生态系统模拟。

模型含有的文件包括：物理参数(.par)、配置(.cfg)、入流(.inf)、出流(.wdr)、气象(.met)，其中出入流和气象文件中均为逐日数据。物理参数文件中包含的参数见表2。气象文件由平均气温、太阳辐射量、云量、日累计降水量、平均风速和平均水汽压组成。

表2 DYRESM主要物理参数

a.模型入流数据。因缺乏河道流量监测数据，本文采用新安江模型[31]计算。张倚铭等[32]采用渔梁水文站2006年1月1日至2012年12月31日和2013年1月1日至12月31日两个时间段的逐日监测流量数据分别对新安江模型进行了率定与验证，率定期与验证期的纳什效率系数分别为0.98与0.95，水量平衡误差百分比均为1%，说明新安江模型在渔梁站及以上子流域内模拟精度较高。而渔梁站以上子流域的坡度、土地利用情况、气候条件等与千岛湖其他子流域非常相似，因此本文直接使用通过渔梁站监测数据率定后的新安江模型，计算了2006年1月1日至2020年12月31日新安江、东源港、武强溪、枫林港和云源港5条主要入湖河道的逐日入湖流量，入流水温则为实测值。

b.模型验证。为确定DYRESM对水温和水位的模拟精度，采用大坝前国控点2018年1月1日至12月31日的实测水位和水温数据对模型进行验证，对比分析水位和水深0.5 m、10.5 m、20.5 m的水温模拟结果与实测结果(图2)，模拟值和实测值之间的相关系数分别为0.99和0.99、0.99、0.92，均方根误差分别为0.29 mm和1.08℃、1.10℃、1.47℃，说明模型对千岛湖水位和水温的模拟精度较高。

图2 千岛湖大坝前国控点水位和水温模拟值与实测值对比Fig.2 Comparison of simulated and measured daily water level and water temperature at Dabaqian

2.2.2 湖泊分析程度LA

LA是计算湖泊(水库)稳定性参数和热力分层指标的程序。程序输入不同监测频率的水温、平均风速、水深等指标，输出温跃层表层深度(MetaT)、温跃层底层深度(MetaB)、湖泊数(Ln)、韦德伯恩数、施密特稳定度(St)、浮力频率(N2)等参数。本文主要采用Ln与St作为水体稳定度参数，其中St反映水体从现状态转换成完全混合状态所需要的能量，在一定程度上可以度量温跃层的强度，该值越大表明湖泊分层越稳定，越接近0表明湖泊越趋于混合状态；Ln表示外界因素(如风力、出入流和人为活动)对水体分层稳定性的干扰程度，其值越小代表水体越趋于混合。本文所用温跃层厚度(MetaTH)为MetaT与MetaB的差值。

3 模拟结果与分析

3.1 情景设置

降雨和蒸发虽然能在一定程度上影响湖泊水位的变化，但引发千岛湖水位变动的主要原因是入湖河道流量和出库流量的波动，因此在配置模拟情景时，主要根据千岛湖出入湖河道历史流量平均值(代表多年平均情况)与水位波动幅度较大情况下的实际值(代表极端水位变化情况)进行设定，共设置两种情情景：情景1代表千岛湖极端水位波动情况，发生于2013年7月1日至2014年7月31日，该情景中5条主要入湖河道的流量通过新安江模型计算得到；情景2中5条主要入湖河道的流量为2006—2020年的多年平均实测入湖流量。两种情景均采用相同的初始水位(104.40 m)，模拟时段均为2013年7月1日至2014年7月31日，气象数据均为该时段内的实测数据，时间步长设为1 h。

情景1水位最大变幅达到11.49 m，2013年7月1日水位为104.40 m，经历231 d水位下降至95.01 m，其中2013年7、8月水位下降速度较快，下降速度分别达到0.09 m/d与0.07 m/d，再历经152 d于2014年7月18日达到最高水位106.50 m；情景2最高水位与最低水位相差5.48 m，从2013年7月1日的104.40 m下降到2014年2月16日的98.92 m，再历经153 d上涨3.58 m达到102.50 m。

3.2 模拟结果

图3、图4与图5分别为两种情景水温随水位的逐日变化、水体稳定度参数和温跃层范围的模拟结果。由图3(a)、图4(a)和图5(a)可知，情景1温跃层在7—9月相对稳定，10—11月逐渐消失，12月水体完全混合。2013年秋冬季，MetaT与MetaB逐渐增大，平均值分别为21.59 m与25.18 m；而MetaTH整体较小，平均厚度仅为2.53 m，低于同时期情景2平均厚度(4.75 m)，最大厚度出现于2013年8月(5.02 m)；St与Ln总体呈下降趋势，平均值分别为5 252.00 J/m2与87.00。2014年1—6月，情景1中水位逐渐上升，而MetaT基本不变，维持在水深6 m左右，MetaB不断增大，5月MetaTH成为该时段的最大值，为8.45 m；7月温跃层位置明显加深，位于水深19.8～24.5 m；2014年下半年St和Ln平均值分别为2 821.00 J/m2与25.90。

图3 千岛湖大坝前国控点水温随水位的变化Fig.3 Chang of water temperature with depth at Dabaqian

由图3(b)、图4(b)与图5(b)可知， 2013年7—8月情景2温跃层较浅，仅存在于水深12～20 m之间，9—12月MetaT与MetaB的平均值分别为25.90 m和29.19 m，温跃层位于水深19.0～35.0 m之间；2013年7—12月，MetaTH逐渐变小，7月平均厚度最大(7.19 m)，12月最小(0.2 m)。2014年1—7月，MetaT从3.85 m变为19.50 m，MetaB从6.10 m增至23.73 m，其中，6—7月二者变化最明显，分别增加9.95 m和7.63 m。

图4 千岛湖大坝前国控点水体稳定度参数Fig.4 Water stability parameters at Dabaqian

图5 千岛湖大坝前国控点温跃层范围Fig.5 Thermocline ranges at Dabaqian

情景1温跃层平均水位98.86 m，St与Ln的平均值分别为4 070.26 J/m2和55.93；情景2温跃层平均水位为101.10 m，St与Ln的平均值分别为4 189.45 J/m2和51.04。情景1与情景2相比，平均水位低2.24 m，St均值低119.19 J/m2，可见极端水位变化情况下，水体稳定度也相应降低。

3.3 典型时段温度垂直分层差异

3.3.1 温跃层持续时间差异

分析两种情景温跃层的形成与消失时间可知，情景1温跃层从模拟开始时间(2013年7月1日)持续到2013年11月20日，在2014年3月23日重新出现，水体混合期122 d；情景2温跃层自模拟开始时间持续到2013年12月5日，于2014年3月16日重新出现，水体混合期102 d。由St与Ln计算结果可知，2013年11月至2014年3月(水体混合期)情景1平均水位比情景2低4.09 m，St平均值比情景2小71.70 J/m2，说明冬季水位降低时湖泊混合相对较强，稳定度较差，而处于高水位且水位波动幅度较小时稳定度较高，温跃层相对稳定。

3.3.2 温跃层结构差异

模拟时段内，2013年7月20—30日，情景1与情景2平均水位分别为99.76 m与103.41 m；2014年7月20—30日，情景1与情景2平均水位分别为106.35 m与101.98 m，水位差异达到最大。选择两情景在上述时段内的水体分层情况进行对比，分析水位差异最大时段内的温跃层结构特征，结果如图6所示。2013年7月20—30日(图6(a))，情景1温跃层出现在水深10～35 m，MetaT与MetaB处水温相差18.07℃；情景2 MetaTH较小，水层的温度变化梯度较大，从水深10～27 m水温下降15.6℃。2014年7月20—30日(图6(b)),情景2从表层到水深15 m的水层温度梯度变化较小，水深15～35 m出现温跃层，水体温度下降15.2℃，而水深35 m以下水体无明显温度变化，此时水体为单一、稳定的温跃层结构。情景1中由于大气降水和河道径流的汇入，水位上涨较快，使得表层水体发生了一定程度的扰动，水深15 m以上的水层水温明显升高，且在水深30～40 m出现15℃左右的“均温层”，水深40～50 m水体温度梯度为0.29 ℃/m，温跃层再次出现，即水位上涨之后MetaB相对情景2有所增大，且温跃层被水深30～40 m处的均温层切断，水体呈从表层到深层依次为“混合层、温跃层、均温层(30～40 m水深)、温跃层、均温层”的双温跃层结构。

图6 模拟水温随深度的变化Fig.6 Changes of simulated water temperature with depth

4 讨 论

千岛湖是一个典型的亚热带单混合水体，正常水位情况下，水体从3月下旬出现分层，至12月初分层消失，温跃层持续时间较长。本文所设置的两种情景中，水位与St、Ln均呈显著正相关关系，情景1相关系数分别为0.81和0.41(样本量为396)，情景2相关系数分别为0.79和0.53。因此，水位越高，MetaTH越大，湖泊水体越稳定，分层越明显。极端水位变化情况下，对温跃层持续时间的影响主要体现在水体稳定度降低和混合期提前，缩短了温跃层持续时间。情景1秋季温跃层比情景2提前15 d结束，而春季温跃层开始日期仅相差5 d，说明水位对水体分层的影响在秋季更明显，这与Wang等[4]对广东省流溪河水库的研究结果一致，进一步说明了本文所采用的模型预测的可靠性。

在国内外具有相同气候背景的其他湖泊中，温跃层持续时间具有与千岛湖类似的特征。例如：韩国Yongdam水库(最大水深70 m)，其温度分层时间为4—11月[33]；云贵高原乌江流域的百花湖水库(最大水深45 m)分层时间为4—10月[34]；流溪河水库(最大水深73 m)温跃层出现在4—12月[4]。但是，包括千岛湖在内的这些湖泊，受气候变化的影响，随着暴雨频率和强度的增强，极端降水事件出现的概率会增大，意味着湖泊径流量和水位波动会更加激烈和复杂化[7,35]，对温跃层的冲击将不可忽视。而在这些极端事件的发生过程中，很难依靠人工监测获取水位和水温垂直分布的信息，因此依靠浮标探头进行自动监测和模型模拟将是获取这些信息的有效手段。相比于探头自动监测技术，数值模型具有经济和高效的特点，且所获数据具有时空尺度上的连续性。因此，在湖泊物理特性的研究方面(例如水温和水量)，模型模拟技术是对人工监测和自动监测的有效补充。

模型模拟结果表明，千岛湖水位较高时水体出现双温跃层结构。通过对影响温跃层的各种可能因子的主成分分析可知，影响MetaTH的第一主成分为入湖河道水温和气温，贡献率达44%，两者与MetaTH相关系数分别为0.78和0.71(样本量均为396)，而入湖流量、降水量和云量为第二主成分(贡献率18%)。可见千岛湖双温跃层结构的形成不仅与气候条件有关，也与水文条件有关。这种双温跃层结构其实普遍存在于深水湖库中，如海南松涛水库[36]和二滩水电站水库[37]均有发现；Liu等[21]研究表明，强降雨条件导致短期水量大量涌入，改变了表层水温，导致水体出现暂时的双温跃层结构。双温跃层结构形成期间，河道径流汇入与其密度相近的水层并发生混合和热量交换，这在一定程度上使得深层水体能得到表层富氧水体的补充，为浮游动植物的生长生活提供溶解氧。同时分层强度的削弱也有利于水体中营养盐的扩散，降低了局部水域营养盐的浓度。在千岛湖，除暴雨过程外，其他气候条件基本稳定，虽然气候变化可能导致气温缓慢上升，但这种效果在短期内几乎可以忽略。因此由暴雨事件、千岛湖发电和水利调度造成的水文条件的急骤改变，进而造成的对千岛湖水体垂直温度结构的破坏以及对水质的影响，值得更进一步深入研究。

已有研究表明，重要的水动力过程(如极端水位变化)会引起上覆水和沉积物的理化性质的变化，进而影响甚至改变底泥污染物的释放规律[38]。当温跃层存在时水体垂直交换受阻，DO的水平分层结构稳定，底部水体氮磷营养盐、金属离子及有机物等污染物在厌氧环境中释放[39]，并可能在沉积物-水界面呈现一定程度的营养盐累积效应。极端水位变化情况下，冬季温跃层提前消失，不稳定的水体可能导致底层污染物向上层水体扩散，上层水体中营养盐富集可能加快浮游植物的生长，并增加水华暴发的可能性[40-41]，给饮用水源安全带来影响[42]。因此，对类似于千岛湖的大型深水湖库，水文条件剧烈波动所造成的水环境的改变和对生态系统的冲击不能低估。随着自动监测技术的愈发先进和模拟技术的愈加成熟，这种影响和冲击的定量评估会变得越来越容易，能给相关管理部门提供制定有关极端气候事件应对策略的科学依据。

5 结 论

a.水位剧烈波动条件下，温跃层持续时间较正常水位波动情况下明显缩短。情景1水体混合期为2013年11月20日至2014年3月23日，情景2则为2013年12月5日至2014年3月16日，情景1相对于情景2温跃层持续时间缩短20 d。

b.高水位时，情景1湖体呈“混合层、温跃层、均温层(30～40 m水层)、温跃层、均温层”的双温跃层结构，情景2则为“混合层、温跃层、均温层”单温跃层结构。

c.情景1平均水位较情景2低2.24 m，St均值低119.19 J/m2，表明极端水位变化条件下出现的低水位情况能显著降低水体稳定度。