大型分层水库翻库特性及溶解氧响应研究

杜彦良，刘小蔚，刘晓波，刘 畅，王世岩，赵仕霖

(1.中国水利水电科学研究院 水生态环境所，北京 100038；2.上海水环境模拟与水生态修复工程技术研究中心，上海 200233；3.长春工程学院，吉林 长春 130000)

1 研究背景

垂向水温分层导致的密度分层现象普遍发生在各种大型湖库中，基于年内分层和混合，早在1957年就有学者提出双循环型、单循环型、局部循环型及混合型的湖库分类[1]。分层水体温跃层界面抑制上下层水体的混合，产生了各种水质指标的分层，衍生生物地球化学分层[2]和微生物分层，形成复杂的生源物质和微生物群落功能的垂向分异机制，同时涉及到氧化-还原热力学、界面平衡动力学、微生物营养生理与代谢动力学等多方面因素的耦合作用。垂向稳定分层结构的变弱和打破被称为“翻库”。翻库期间一些水体表现为颗粒态浓度增加、透明度下降[3]，有时伴有难闻气味，水体上层水质变差[4-6]，溶解氧降低，甚至造成鱼类与底栖动物数量下降[7-8]。湖库翻转也是表层和底层氮磷营养盐垂向再循环的重要路径，尤其为春季藻类提供了早期生长的磷，当湖泊秋冬季翻转缺失后，湖区表层磷的急剧消耗，使得浮游植物春季水华减弱甚至消逝[6]。翻库是内陆湖库底层水体复氧的主要途径，也是温室气体CO2、CH4、N2O等的重要排放方式之一[9-10]。水体翻转是季节性种群控制的主要机制[11]，随着气候变暖，导致一些深水湖库底层充分复氧频次减小，强度减弱，导致有毒藻类增加，对少泥质湖库，翻库相比富营养水平对深层水体的溶解氧(DO)及水生态影响更大[12]。

关于湖库分层及终结的时间，通过高频水温监测可观测得到，不同文献有不同的水温判定阈值：(1)表层水温(Tsur)与底层水温的温差(Tbot)小于0.4 ℃，即 |Tsur-Tbot| ≤ 0.4 ℃(或垂向水温梯度 ≤ 0.01 ℃/m)[13]；(2)|Tsur-Tbot| ≤ 1.0 ℃[14]。当垂向数据不足时，表层水体外观变化，或表层水质、水动力数据分析也可得到相应时间[15]。Blanton通过对多个水库的数据分析，对于某一水库，平均深度似乎是影响翻转唯一重要的形态参数[16]。Lewis等[17]认为纬度、海拔高程和深度并不是影响湖泊混合的唯一因素。Nuernberg[18]提出夏季底部滞温层的平均水温对秋季翻转日期具有相关性，并建立简单三变量模型公式预测湖库翻转日期，通过引入平均深度和经纬度对模型进行改进。水动力变化是翻库的原始动力，但是系统全面的水动力监测并不常见[19]。通过一维[5，12，19]、二维[20]及三维[19]模型研究翻库动力及环境指标变化、驱动要素及时间节点等，逐渐成为主要方法。采用模型分析影响翻库的主要因子及贡献占比，大致归结为气温持续引起的垂向对流、风驱动、入流扰动等，对微咸湖泊可能蒸发作用占较大比重。中小型或弱分层湖库，上述外界的扰动即使不大也可能造成热稳定性丧失，翻库成因似乎更复杂，成为大量文献研究对象。对于大型深水湖库，有更强的热分层稳定性，在分层将结束的时段内受外界影响将导致提前或滞后，总体上翻库的规律性较强。

本文在北方大型深水双循环水库潘家口水库进行原位监测基础上，采用垂向二维CE-QUAL-W2模型，对水库年内的水动力、水温、水质、藻类、溶氧进行模拟，分析秋冬表层及底层水温及DO翻转在水库中的演进过程，通过不同年型模拟，认识水环境响应特征，分析翻库与水深及稳定性的相关关系，揭示大型深水水库翻库规律。

2 材料与方法

研究区域为潘家口水库，位于河北省北部(39°10′N—42°40′N，115°30′E—118°45′E)，水库总体上为狭长的河道形水库，水库北部有滦河和柳河汇入，中游瀑河支流汇入。2000年以来，受养殖渔业发展水库氮磷浓度不断升高。2016年后水库进行了养殖网箱的全面清理，2017年后水质有逐渐改善的趋势，TP下降显著，基本满足Ⅲ—Ⅳ水标准，TN为Ⅳ—Ⅴ类水，库区为轻度富营养化状态。对水库原位垂向监测和常规水质监测数据的收集分析，研究区域及监测断面位置见图1(a)。

图1 研究区域、断面位置及模型计算网格

采用垂向二维开源代码的CE-QUAL-W2模型，空间离散网格见图1(b)。模型入流分别有北部入流和东部的瀑河入流。南部为水库大坝，由于建造时间较久远，没有分层取水措施，水库的发电取水位置位于大坝底部。

模型基本方程为：

(1)

(2)

(3)

(4)

式(1)为连续方程；式(2)为水平方向动量方程；式(3)为垂向动量方程；式(4)为水面方程。式中：U为水平流速，m/s；W为垂向流速，m/s；g为重力加速度，m/s2；q为单位体积通过的流量，m3/s；B为水面宽度，m；α为底坡与水平夹角，rad；P为压力，Pa；τxx和τxz分别为x和z方向的切应力，N；ρ为密度，kg/m3；η为水面高程，m；h为水深，m。

水温及溶质输运的通用方程为：

(5)

式中：Φ为水温(Tw)或溶质(如DO等断面宽度)平均成分浓度，℃或mg/L；Dx为温度或水质成分纵向扩散系数，m2/s；Dz为温度或水质成分垂直扩散系数，m2/s；qΦ为横向流入或流出的物质或热流量，C/S或g/(m3·s)；SΦ为横向平均源汇项，g/(m3·s)。

模型耦合计算水动力、水温、水质及藻类动力学，其中水质主要计算指标为DO、悬浮固体、磷、氨氮、硝氮、DOM、POM及藻类等，藻类分别模拟了硅藻、蓝藻和蓝绿藻三种藻类。模型输入气象数据采用实测高密度小时气象数据，并对缺失数据进行补插。水动力的边界条件数据为逐日出入库流量，水质输入采用常规逐月监测数据，计算边界进行线性内插加密。

模型率定数据为水库常规监测的逐日水位、水温及逐月常规水质监测数据，计算时长为2018年 1月1日至2020年12月31日，采用2018年及2020年全年观测数据做模型的率定验证。垂向水温率定采用2020年间多次的YSI-EXO原位监测的结果。本文主要展示水温、DO等率定验证部分结果。2018年和2020年水库坝前实测水位见图2。

图2 2018年与2020年水位与逐日实测水位比较

由图1可见计算与实测水位吻合度较高，R2=0.999，不仅再现了水库的进出水量的平衡，确保出入流产生的水库热量收支账。在不同入流及调度过程下，2018年和2020年水库水位有较大差异，年平均水位分别为221.5 m和207.3 m，年内水位变幅分别为8.1 m和12.1 m，2018年有较高的运行水位和较小的水位年内波动。

2018年和2020年模型计算坝前表层水温与逐日实测值对比见图3，模型计算坝前浅表层DO与月频次常规监测数据对比见图4。

图3 2018年与2020年坝前表层模型计算与实测水温比较

图4 2018年与2020年坝前计算与实测DO比较

通过误差分析，模拟表层水温，R2=0.92，数据吻合度好，模型计算表层DO的R2=0.87。由图4可见，潘家口水库逐月常规监测结果显示2020年在11月表层DO出现异常低值现象。据管理部门观测，在某些年份的冬季，水库翻库带来持续若干天的水体浑浊，氮磷浓度升高，水质变差现象。2018年的监测数据该现象却不显著。

项目组在2020年对水库进行了6次的全库区的原位YSI-EXO垂向监测，选择水温初步分层(5月)、夏季及夏末稳定分层(7月和9月)、分层消失(11月)合计4次实测结果用于模型垂向结果验证，点位为潘家口(PJK)断面，垂向指标分别为水温(T)和DO，见图5。

图5 库中潘家口点位垂向水温模拟与实测对比

垂向水温误差分析结果R2=0.96。更多不同断面及不同日期的模型计算结果与监测数据的对比参见文献[21]。2020年燕子峪(YZY)与坝前(DAM)断面的垂向DO变化过程见图6。图中黑线为变温层的下缘位置，即从底层向表层扫描，|ΔT|≥1 ℃/m的位置。

图6 2020年燕子峪和坝前断面垂向DO年内变化

由图5和图6结果可见，水库在5月水温开始分层，但DO并无显著分层，随着分层结构越来越显著，温跃层变厚，底层DO逐渐变小，水库燕子峪水域温跃层下缘的上方DO略低，坝前不明显。至11月7日表层与底层的垂向水温基本一致，水库的密度分层已基本打破，但垂向的DO和底层的DO存在差距。该结果与文献[22]中千岛湖2月垂向水温和DO监测形态相似，即当垂向分层被打破后，表层与底层的DO均匀混合状态滞后于温度或密度混合。

3 分析与讨论

3.1 年内表层与底层水温及DO变化为了解水库表层及底层的水温及DO的年内变化，以及变化的特性，通过模拟分别做出2018年及2020年水库坝前水面及底部的水温与DO的变化过程，见图7，图中“Sur”表示水面，“Bot”表示库底。潘家口水库显示出双循环水库特性，垂向水温在气温最冷的2月及春夏季均出现温度分层。当水温分层，表层与底层的DO也开始发生变化，水库的上层，大气对水体的复氧条件好，DO浓度较高，受气温影响，冬季高于夏季，底层DO主要受水库底泥耗氧影响，水温分层后，DO浓度逐渐下降，当水体垂向均匀混合，底部DO迅速上升，翻库使得水库底部溶氧水平基本恢复。

图7 2018年及2020年坝前水面及底部水温和DO变化过程

坝前表层DO受大气复氧及藻类的光合作用，维持在较高水平，冬季低水温有较高的DO浓度，夏季DO略低。底部水体冬季DO在垂向上均匀混合，当水体开始分层，底层水体与表层的对流扩散减弱，滞温层水体内还原物质的不断耗氧，底层逐渐低于表层，在夏季分层时间长，底层DO逐渐耗尽，2018年及2020年底层DO浓度低于2 mg/L的天数分别达到66 d和40 d。

受水库的运行调度、气象条件的影响，水库的水温和DO过程不同。2018年相比2020年有更高的运行水位，底层水温维持在较低水平，越高的水位导致翻库日期推后，2018年的表层DO的变化小于2020年的变化。2020年监测及模拟结果均显示出表层水体的低DO浓度现象，而2018年却不太明显。

水库翻库时段低水位运行，沉积水平较高的水库可能产生较为剧烈的水质污染现象，带来水生态负面影响。在南方的某水库我们也观测到翻库的日期通常比北方滞后，监测氮磷营养盐不高的情况下，当翻库时，水库运行水位过低，在早春季节引发水华。但对某些大型深水水库，在气候变暖背景下，过高的运行水位也可能导致翻库作用减弱，底部复氧的不充分，同样引发生态环境问题。

3.2 稳定性分析深大水库分层通常持续较长的时间，在稳定期间短暂的强对流天气基本不能打破分层结构，水库变温层深度及厚度在全库区的分层水域基本一致。选择水库中的潘家口断面，对2018年及2020年水库年内的垂向水温过程进行描述，见图8。

图8 2018年与2020年潘家口断面的垂向水温年内变化

水库垂向水温的变化(见图7)表明，2018年相比2020年，7月下旬上层温水层厚度明显变厚，变温层(垂向水温梯度≥1 ℃/m)深度大，2018年变温层最大水深与2020年比下移了近10 m。原因是综合的，但水位的变动是主因。2018年7月中旬水库的水位开始抬升，而2020年水位在8月中旬开始抬升，虽然2018年相比2020年有较大的年内水位变幅，但影响夏季分层结构的是夏季混合循环期间的运行水位及动力学变化。

为定量描述分层水体的稳定性，1955年Lorenz[23]最初推导了流体的总势能(重力和内部势能)，很多学者通过可用势能指数(Available Potential Energy，APE)研究分层流体偏离参照密度分层结构位移偏差下的势能[24]，APE计算方程见式(6)，其中密度计算为式(7)。

(6)

(7)

式中：APE为可用势能指数，J/m；z为所在高程，m；ρ对应水深处的水体密度，kg/m；ρr为参考位置的水体密度，kg/m；zx为参考总高柱，m；ρT为水的密度随温度的函数，kg/m3。

水库垂向均匀混合状态为参照状态，表层水体与底层水体密度基本相同，APE指数表征了水体完全混合所需要的势能，越高则代表水体密度分层结构越稳定，反之说明稳定性较弱，与水深、温跃层厚度、温度梯度有密切关系。选取参考位置为底层水体，通过模型计算结果，进一步计算不同断面的APE年内变化过程。2018年及2020年的不同点位APE变化见图9。

图9 2018年和2020年不同断面APE年内变化过程

因潘家口水库的纬度较高，在1—2月垂向形成表层水温低底部水温略高的结果，但APE不大，在3—4月垂向混合后，开始夏季的分层演进，APE升高，在APE增大的过程中，全库区APE变幅不大，在7—9月先后达到最大值。水库不同断面的最大稳定性不同，坝前及深水水域的稳定性高，库尾稳定性较低，越深的水域，APE最大值出现得较晚。不同来流及运行调度水位下，APE存在差异，2018年运行水位高，年内最大APE大，2020年水位低，对应较低的APE值。2018年靠近库尾的清河口断面的年内最大APE值大于2020年坝前值。2018年表温层厚度大，产生了较大的APE指数。总体上APE与运行水位的相关性程度更强。

3.3 翻库日期与稳定性关系翻库的日期的预测，对水库管理及科学研究均有重要意义。本文采用|Tsur-Tbot| ≤ 1.0℃[21]定义，并定义DO完全翻转时|DOsur-DObot|≤ 1.0 mg/L。通过模型反演水库水温、水质及DO的变化，对于狭长型的河道水库的翻库，并不是在全库同时发生的，而是由库尾浅水段逐渐行进至坝前，在坝前水深相差不大的水域，基本同时发生。在表1中，分别列出水库断面位置及最深处高程，在2018年和2020年水温及DO的翻转日期，该日期分别用月/日和该年的计算日数(Julian Day Number，JDN)表示，儒略日数是一种不用年月只从起点累计日数的方法，每年的起点天数始于1月1日。各断面位置见图1。

表1 模型计算水库不同断面翻库日期及DO混合日期

2018年从库尾清河口到坝前，翻转过程为57 d，2020年该过程为98 d，2018年比2020年翻库经历的时间短了41 d。2020年库尾由于水深较浅，在气温较高时期就发生混合，到坝前时间反而更长些。在2年中，燕子峪、潘家口及坝前水域几乎是同时翻库，即水表及库底水温基本一致。

2018年水库水位较高，DO的翻转与水温相比推迟了10 d，2020年水位较低，DO翻转只比水温推迟了7 d。高水位运行，表层和底层DO交换垂向空间大，因此时间也更长，表层DO的变化不似2018年剧烈。潘家口水库纬度较高，底部富含高有机物质，库区水流缓慢，冬季底层水体的复氧水平高。

通过模型可得到比常规监测更为精细的数据，进一步研究APE与翻转时间的相关性。因水库年内发生2次翻转，对于秋季的翻转的有效时间开始于春季均匀混合后。计算2018年和2020年夏季分层后的各断面平均的水深及APE值，与表1中的翻库日期做相关性分析，纵坐标分别为平均水深和APE值，横坐标为年的计算日数，见图10。

图10 翻库日期与APE关系

参考文献中对翻库的预测成果，Blanton[15]结论认为翻库日期与水深相关，分层期间的平均水深与翻库时间关系见图11(a)。Nuernberg[17]推荐计算翻库的公式为：

图11 翻库日期与各断面平均水深和Nuernberg公式计算时间关系

(8)

latad=latx+ax/100×alt

(9)

水库分层后，水库储存越多的热能，翻库的日期越滞后。断面在稳定性较差的分层末期，可能一场风或降温或入流脉冲等导致水体的翻转，具体的日期或提前或滞后。分层期间内的平均APE的大小与水库的翻库日期呈正相关。同年中，库尾断面平均APE较低，翻库日期早，坝前水体翻库迟。虽然2018年与2020年运行水位及水温结构变化过程不相同，但可将不同年份不同断面的翻库日期统一在以APE为指标相关关系中。APE是主导水库翻库日期的主要因子。

4 结论

湖库的翻库具有生态环境学意义，通过水库的原位监测及模型手段，对大型深水水库进行研究，过低和过高的水库运行水位均会通过翻库对生态产生影响。

大型深水水库的翻库过程始于库尾或浅水水域向深水水域演进，狭长形潘家口水库的翻库过程可持续2或3个月。翻库时水库的水温先在垂向上达到基本均匀，DO垂向上的均匀混合通常会滞后于水温。强烈翻库使得水库底层缺氧区复氧。当水位较低时，表层水体受翻库影响更大，表现为较短时间内DO迅速下降等。

关于翻库日期及分层稳定性，水库的运行水位是关键，夏季混合循环中的水位变动，影响温跃层的深度及上层温水层的厚度，也决定了水体稳定性指数的大小，可通过可用势能指数(APE)量化水库分层稳定性。水库不同位置断面的翻库日期与夏季循环期的平均APE指数大小有较好的相关关系。

研究成果将支撑水库优化运行调度、水源地生态环境保护和水环境科学管理等的实践。