张 晨,宋迪迪,廉铁辉

(1：天津大学水利工程仿真与安全国家重点实验室,天津 300350)(2：天津市引滦工程于桥水库管理处,天津 301900)

引调水工程是解决水资源时空分布不均、缓解资源性缺水问题、优化水资源配置的重要举措[1]. 在跨流域引水系统中，调蓄湖库一定程度上能够缓解水资源供需矛盾，但其供水水质也会受到引调水、当地地表水等不同来水的影响[2]. 总磷(TP)是调蓄湖库水质优劣的显著性指示指标[3-4]，与湖库水生态之间有着密切联系，弄清哪些因素影响湖库磷滞留量是保障供水水质安全的前提. 前人研究成果分析了外部营养负荷、单位面积磷负荷率、单位面积水力负荷、平均湖深、水力停留时间、冲刷率等一系列因素对磷滞留的影响[5-7]. Vollenweider模型[8-10]是众多学者针对这一问题得到的典型成果，该模型基于湖泊总磷负荷质量平衡原理，提供了一种基于湖库形态测量和水力数据预测湖泊水体磷滞留的工具. 模型指出外部TP负荷、水力停留时间、冲刷率是磷滞留的重要影响因子[3]. 后续研究证实了该观点，外部营养负荷决定了湖泊的总体水质，但沉积物对磷的内部循环起着核心作用[11-12]. 水力停留时间通过流动的快慢影响营养物和有机物的供应和流失，冲刷率为水力停留时间的倒数[13]，它们影响湖库中TP浓度[3,14]. 一般认为，湖泊的水力停留时间越长、冲刷率越小，将有更多的磷滞留，湖水TP浓度越低[15-16]. 反之，在一些调蓄水库中，水力停留时间较短(冲刷率较大)，磷往往未参与内部循环便被带出. 但是，Wu等[17]对太湖、巢湖和鄱阳湖研究发现，预测TP浓度随着水力停留时间的增加(冲刷率的减少)而升高，这是因为沉积物中的磷还原释放到水中导致TP浓度升高.

“引滦入津”工程是我国早期建成的大型跨流域调水工程，解决了天津城市生活及工业用水需求. 南水北调通水前，于桥水库是天津市唯一饮用水水源地，也是“引滦入津”工程的重要调蓄水库，年均受纳滦河水系引水5.4亿m3(2001-2014年). 2014年底，南水北调中线一期工程通水后，天津市的水源结构发生变化，主要饮用水源从滦河水变为长江水(引江水不调入于桥水库). 引水结构的变化导致外调滦河水的需求发生改变，进而影响于桥水库受纳引水量. 对于桥水库而言，引水量的变化影响外部营养负荷的改变，水库运行方式和出入库流量的变化也影响水库TP浓度. 已有研究表明，于桥水库处于中贫营养水平、沉水植物为主的清水状态，但沉水植物物种单一，清水稳态易发生转变(破坏)[18-19]. 同时，根据监测结果，水库TP浓度自2013年之后有升高趋势，若一旦超过清水转换阈值将突变为藻型浊水状态[20]，修复富营养化将十分困难. 最新成果表明，城市的用水需求、引水量、土地利用和气候变化均将深刻影响于桥水库的水质状况[21]. 综上，在关注于桥水库TP浓度时，弄清引水结构变化引起的入库水量、外部TP负荷、水力停留时间改变对磷滞留的影响机制是重要科学问题.

本研究基于于桥水库2001-2018年实测资料，利用Vollenweider模型计算TP浓度动态过程，分析入库水量、外部TP负荷、水力停留时间三个因素对南水北调通水前后水库磷滞留量的影响，并进一步估算双水源新情势下的生态水量，为保持于桥水库中贫营养状态提供科学的运行管理建议.

1 材料与方法

1.1 研究区域概况

于桥水库(40°00′～40°04′N，117°26′～ 117°37′E)位于天津市蓟州区城东4 km处，自1983年“引滦入津”输水工程通水后，是天津市人民生活饮用及工农业用水的主要水源地，也是“引滦入津”输水工程的重要调蓄水库[22]. 控制流域面积2060 km2，占整个州河流域面积的96%，图1所示. 最大回水长东西约30 km，南北宽8 km，最大淹没面积250 km2，正常蓄水位21.16 m时淹没面积86.8 km2[23]，汛限水位19.87 m，死水位15 m. 浅水库，平均水深4.6 m，总库容15.59亿m3. 水库控制流域属于温带大陆性季风性半湿润气候，境内雨量充沛，多年平均降水量为748.5 mm，主要集中在6-9月，水库多年平均径流量5.06亿m3[24-25]. 于桥水库是一座典型的山谷与平原过渡的水库[26]，库区南岸地势陡峭受原河道深泓控制水深大，北岸地势则相对平坦水深浅[27]. 由沙河、淋河、黎河三大支流汇合而成，其中黎河为引滦输水通道，与沙河汇流果河后随果河汇入于桥水库.

图1 于桥水库流域及水文监测站Fig.1 Yuqiao Reservoir basin and monitoring stations

1.2 数据来源与处理

1.2.1 水量 黎河逐日总水量来源于前毛庄监测站，为黎河自产水量与引水量之和，2001-2018年平均径流量为5.65亿m3，引水量为4.79亿m3. 沙河逐日流量来源于水平口监测站，为自产水量，平均径流量为0.95亿m3；淋河逐日水量来源于龙门口(2012-2018年为淋河桥站)监测站，平均径流量0.16亿m3，由于入库水量远小于黎河和沙河，且年内常处于断流状态，故本研究中忽略淋河. 果河水量即入库水量WIN，为黎河与沙河水量之和，多年平均径流量为6.60亿m3.

1.2.2 TP负荷 各监测点TP浓度逐月数据源自水库管理部门. 沙河水质取样点位于沙河桥，2001-2018年TP浓度平均值89.5 μg/L；入库TP浓度TPIN取样点位于沙河、黎河汇流后的果河桥，平均浓度201 μg/L；黎河TP浓度根据磷质量平衡反算得出；库区TP浓度TPL由TPIN计算得到，详见下文公式(3)，年均值36.8 μg/L. 沙河、入库、黎河、出库各点TP负荷由水量与浓度计算得到.

1.2.3 数据处理 沙河、黎河2013-2018非汛期部分月份缺失的水量数据由皮尔逊III型曲线插值确定，再根据水量平衡，结合已有的引水数据对不合理插值结果进行修正，修正后月数据累加后得到的年水量与水资源公报中年水量相对误差为4%.

1.3 Vollenweider模型

Vollenweider模型用于预测已知营养负荷输入和水力停留时间的湖泊营养浓度，已被广泛用于估算湖泊和水库中的TP浓度，数学表达式为[7]：

(1)

式中，TPL代表水库TP浓度(库心测点)，μg/L；L表示总磷的单位面积负荷率，单位mg/(m2·s)，由入库总磷负荷TPLIN除以库区面积A得到，t/a；其中，TPLIN等于入库流量Q与入库总磷浓度TPIN的乘积；H为平均水深，单位m；Sr为冲刷率，a-1，等于入库流量Q比库容V；σ为沉降率，代表湖泊磷流失的一阶速率系数，a-1. 以入库总磷浓度表示，式(1)化为：

TPL=TPIN/(1+στ)

(2)

式中，τ为水力停留时间，等于V/Q，a; 根据Brett和Benjamin[3]研究成果，本文假设σ=aτb，其中a代表湖内磷衰减的无量纲参数，b代表通过水文冲刷输出磷的无量纲参数，两个参数由式(3)和实测数据统计回归求得(见下文).

TPL=TPIN/(1+aτb+1)

(3)

1.4 库区总磷滞留量

营养物质进入水体受生物、物理及化学过程的影响而截留或通过形态转化损失，从而影响物质向下游的输送通量，这些过程广义上被称为“滞留作用”[28-29]. 沉积物中磷的滞留量常由以下质量平衡的差值计算得出：磷滞留量=磷总负荷-磷损失量-磷转换量[30]. 由于磷损失和磷转换均发生在水库内部，大部分成为内部负荷，结合Vollenweider模型，故本文将磷损失量及转换量一同视作滞留量，并用RP表示，定义为：

RP=TPLIN-TPLOUT

(4)

式中，TPLOUT表示出库TP负荷(t/a). 在湖库混合良好，处于稳定状态的假设条件下，出库TP浓度TPOUT等于TPL[10]，TPLOUT计算为出库流量与式(3)算得的库心TP浓度TPL的乘积. 滞留率为RP与TPLIN的比值(%).

2 结果分析

2.1 南水北调通水前后WIN和TP负荷变化规律

2.1.1 入库水量WIN变化 从WIN年际变化来看，由图2(a)可知，2001-2014年，南水北调通水前，WIN年均值6.88亿m3，其中流域自产水量1.51亿m3，引水量5.37亿m3. 2001-2005年WIN相对较少，2006-2013年略有增加，但波动幅度较小，直至2014年引水量大幅提高，WIN达到峰值11.08亿m3. 2014年南水北调通水后，WIN年均值为5.16亿m3，其中自产水量2.39亿m3，引水量2.77亿m3，且2016年WIN为历年最低，仅有3.04亿m3. 对比南水北调通水前后，年均WIN减少了1.72亿m3，其中，引水量减少2.6亿m3是导致WIN大幅减少的主要原因. 由此表明，引水结构变化致使于桥水库调水需求发生改变，入库水量剧减.

图2 2001-2018年于桥水库入库水量Fig.2 Inflow volume flow into Yuqiao Reservoir during 2001-2018

WIN来源比例方面，如图2b所示，2001-2014年，南水北调通水前，沙河、黎河自产水分别占总水量的11%和10%，引水比重较大，约为79%. 南水北调通水后，沙河、黎河自产水分别占总水量的29%和15%，而引水比重显著下降为56%. 南水北调通水前后，沙河、黎河自产水比重分别提高18%和5%，引水比重显著下降23%. 水量来源方面发生变化，沙河自产水对水库来水的贡献增加.

进一步分析引水期和非引水期WIN变化，如图2(c)所示，引水期WIN年均值为5.53亿m3，引水量4.79亿m3，沙河水量0.44亿m3，黎河水量0.30亿m3. 其中，2004年水量略少. 非引水期，WIN年均值仅为0.97亿m3，不足引水期1/5，沙河水量0.48亿m3，黎河水量0.49亿m3，二者持平. 值得注意的是，2012年沙河水量远超黎河，为历年最大，达到1.89亿m3. 对比引水期与非引水期，沙河水量始终在0.45亿m3左右. 非引水期黎河来水平均占比51%，沙河49%；引水期和非引水期黎河自产水和沙河比重基本持平(图2d). 入库自产水量方面，引水期与非引水期沙河自产水量无显著变化，而黎河非引水期自产水量增加了63%.

2.1.2 入库总磷负荷TPLIN变化 由图3(a)可知，南水北调通水前，TPLIN年际变化相对平缓，整体水平较低，年均值为111.42 t/a. 2004年出现第一个小高峰达到154.70 t/a，2008年降到历年最低42.65 t/a后便逐年上升，并于2014年达到最高峰249.61 t/a. 南水北调通水后，波动较为剧烈，年均值上升至163.26 t/a，较通水前增加了51.84 t/a. 其中2016年达到历年最高254.85 t/a，2017年又骤降为55.91 t/a. 引水TP负荷与TPLIN类似，于2008、2014年分别达到南水北调通水前最小值和最大值，2016年达到通水后峰值，2017年又大幅降低. 由此可见，引水TP负荷是影响TPLIN的主要因素. 为探讨引水TP负荷对TPLIN的影响程度，对引水期和非引水期TP负荷进行对比，如图3(b)所示. 引水期年际变化较非引水期剧烈，年均TPLIN为102.45 t/a，占入库TP负荷的83%. 进一步分析自产水TP负荷，引水期自产水TP负荷年均仅为11.54 t/a，占引水期TPLIN的5%，非引水期为20.44 t/a. 从全年角度来看，TPLIN的74%来源于引水，36%来源于自产水. 可见，本流域内的自产水携带的TP负荷相对较少，外流域引水携带的TP负荷是TPLIN的主要来源.

图3 2001-2018年于桥水库入库总磷负荷Fig.3 Total phosphorus loadings flow into Yuqiao Reservoir (TPLIN) during 2001-2018

于桥水库引水期一般分为汛前和汛后两个时段，大部分年份于4-6和9-12月引水，共持续5～6个月. 结合引水情况分析，2004、2016年均为枯水年，但引水期较短(不足4个月)，短时间高水量导致河床底质扰动增强[31]，加剧了枯水年河流中的TP浓度，故导致TP负荷较大. 2017年仅8月份引水，上游潘家口水库、大黑汀水库采取了清除网箱养鱼等措施，TP浓度大幅减少，有效降低了TP负荷. 除上述原因造成通水后年均TP负荷仍然较大外，综合图2(a)和图3(a)，笔者发现入库水量的减少并未使TP负荷减少，入库水量和TP负荷的关系是怎样的？

图4 入库水量与入库总磷负荷关系Fig.4 Relationship between inflow and total phosphorus loadings

为进一步研究入库水量与TPLIN的关系，对二者进行回归分析，如图4所示，发现TPLIN随WIN呈先降低后升高的趋势. 当WIN较小时，如小于6亿m3，入库水量对TP的稀释作用较为显著，随着水量增加入库TP负荷减少，即外流域引水稀释本流域TP负荷的作用较携带TP更为突出；当WIN大于6亿m3，引水携带TP负荷造成两者显著正相关. 同时，图4回答了前述问题，南水北调通水后年均入库水量为5.16亿m3，TP负荷随水量减少而增加. 虽然该种关系受水库水环境管理措施和年尺度统计数据的影响，在个别年份具有一定的不确定性，但整体趋势如上所述.

2.2 水库磷滞留过程的机制分析

2.2.1沉降系数a、b确定 Larsen和Mercier[32]发现将磷沉降模拟为水力停留时间τ的函数σ=aτb时预测TPL与实测TPL拟合度最高. 随后学者们发现总磷损失率与湖泊水力停留时间的倒数呈正比[33-34]，Brett和Benjamin[3]研究得到a、b取值决定于具体的湖泊. 为得到适用于于桥水库的a、b值，图5(a)绘制了先将Brett和Benjamin研究结果a=1.12±0.08年-0.47，b=-0.53±0.03代入公式(3)计算得到的TPL与实测TPL的对比图，可以看出预测TPL与实测TPL呈显著线性相关，然而，78%的预测值偏向高估. 这是因为当通过τ计算营养物浓度数据时，长的水力停留时间(> 0.4年)会使湖泊的预测TP浓度倾向于高于实测值，而当水力停留时间<0.04年时预测值低于实测值[35]. 于桥水库的τ历年均值为2.33年，因而预测值较大. 基于σ=aτb的假设，利用实测值非线性回归得到a、b取值：a=2.203,b=-0.641，如图5(b)所示. 低估的预测值降至44%；R2值为0.73,表明该假设合理；均方根误差(RMSE)从0.06降低至0.03，表明预测值更接近实测值,参数选取合理. 进一步利用该模型计算库区TP负荷滞留量RP,并与实测推算值对比，R2=0.94(图5c)，说明模型精度较好.

图5 预测库区总磷浓度与实测库区总磷浓度、预测库区总磷滞留量与实测库区总磷滞留量对比图Fig.5 Comparison of predicted total phosphorus concentration and observed total phosphorus concentration，and predicted total phosphorus retention and observed total phosphorus retention

图6 2001-2018年于桥水库库区总磷滞留量Fig.6 Total phosphorus retention in Yuqiao Reservoir during 2001-2018

2.2.2 库区总磷滞留量RP变化 根据图6，RP总体波动范围较大. 南水北调通水前，RP波动范围为30.94～153.59 t/a，年均值76.38 t/a，分别于2008、2014年处于最小值和最大值. 南水北调通水后，年际变化剧烈，2016年达到历年最大值215.86 t/a后2017年又骤降为44.66 t/a，年均值上升至120.47 t/a，RP增加了44.09 t/a. 分析引水期和非引水期RP，发现引水期RP随时间变化情况与年均RP保持高度一致，且基本重合. 图7中，引水期TPLIN与RP的线性正相关(R2=0.91)说明引水期高水平TP负荷流入导致更多的磷滞留，即引水期决定了库区全年TP负荷滞留水平. 非引水期RP除2001、2012、2016年外始终维持在较低水平(接近于0)，此3年汛期(7、8月份)入库水量均较多，入库TP负荷较大(图2b)，故RP相对较高. 另外，非引水期RP多年出现负值，这是因为在非引水期TPLIN显著减少的情况下，出库水量受供水需求依然维持在较高水平，当出库负荷略大于入库负荷，便出现RP为负值的情况.

图7 入库总磷负荷和库区总磷滞留量关系Fig.7 Relationship of total phosphorus loadings and total phosphorus retention

图8 水力停留时间与库区总磷滞留量的关系Fig.8 Relationship between hydraulic retention time and total phosphorus retention in the reservoir

2.2.3 水力停留时间τ对磷滞留的影响机制分析τ与RP之间的相关关系如图8所示. 通常，水力停留时间增加有利于水体自净能力[36]. 从图8可以看出，τ< 0.23年时，RP与τ之间存在显著的正相关关系，该结果与Smith和Kalff[37]指出的随水力停留时间增加，水库对磷的滞留能力增强的结论一致. 但是，随着τ增加，正相关性下降，当考虑全部τ时，RP与τ呈负相关. 这是因为入库水量WIN也可能影响磷滞留.τ<0.35年时，WIN较大，此时大量WIN携带丰富的TP进入水库. 当入库TP浓度较高时，水库底泥主要表现为“磷汇”[36]. 因而较小的τ促进沉积物磷的富集，导致磷被大量滞留在库中. 有研究表明[5,38]，具有较小τ的湖泊含有更多容易与矿物结合的颗粒磷，而这种颗粒磷更容易沉淀. 所以，当τ< 0.35年时，RP较大且与τ表现为正相关. 而τ> 0.35年时，WIN较小，入库TP负荷也随之较少，例如于桥水库入库负荷均值由临界τ(0.35年)前23.74 t/月降低为6.62 t/月. 同时，TP浓度自果河桥至坝前沿程递减，并且在水库进口段存在较大的浓度梯度[23]. 进一步分析TPLOUT发现，出库负荷始终处于较低水平，临界τ前、后均值分别为2.71和2.56 t/月，变化不大. 由此可知，TPLIN骤减决定了RP迅速减少.

2.3 生态水量估算与水库管理建议

夏军等[39]提出生态需水是指维系一定环境功能状况或目标(现状、恢复或发展)下客观需求的水资源量. 水利部海河水利委员会编制的《海河流域水资源保护规划》提出生态水量原则上采用Tennant法[40]计算. Tennant法主要基于水文资料和年平均径流量百分数来计算生态水量，根据Tennant法推荐表得知多年平均径流量的60%～100%为最佳生态水量，平均径流量的200%为极限值. 本文将果河径流量视为于桥水库入库水量WIN，2001-2018年WIN均值为6.60亿m3，则于桥水库入库生态水量最佳范围为3.96亿～6.60亿m3，最大不超过13.20亿m3.

由图6、图7可知，RP受到TPLIN和τ的综合影响，且RP在引水期随τ升高有明显增大的趋势，在引水量较少时随TPLIN降低而减小. 由此可见，引水期RP主要影响因素为τ，非引水期或引水量较少时主要受TPLIN的影响. 同时，出库TP负荷TPLOUT主要受出库水量的影响. 因此，针对不同的入库生态水量，可以通过调节出库水量达到控制磷滞留量RP的目的，从而实现库内TP负荷(内源)的管理控制，即在双水源新情势下防治于桥水库富营养化趋势.

综合以上分析，笔者期望在已知生态入库水量情况下，估算出双水源新情势下的生态出库水量，以保持于桥水库处于中贫营养状态. 根据图5可知，南水北调通水前库区TP负荷滞留率(RP/TPLIN)年均值为69%，通水后上升至77%. 库区总磷浓度表明，2001-2014年TPL年均值为41.4 μg/L，处于较低水平，2015-2018年TPL年均值则升高了一倍. 江雪等[26]研究表明，于桥水库上游径流入库口处从1980s起磷的沉积量迅速增加. 2015年库区部分区域有富营养化趋势，沉积物吸附磷增多，近年来库底淤泥对污染物的沉积能力已经趋于饱和[41-42]. 滞留于沉积物中的磷越多，上覆水TP浓度越低，水质越好；但当磷滞留量过多，在厌氧条件下沉积物又会变成“磷源”，还原释放TP到上覆水中. 为保证沉积物仍为“磷汇”，本文假设磷滞留保持在通水前的水平，即取滞留率为通水前历年平均值计算生态水量. 表1给出不同入库水量WIN时，生态水量WOUT的估算值. 例如当WIN= 6亿m3时，在建立的WIN与TPLIN、τ关系式基础上，可得到TPLIN、τ；根据公式(3)由τ得到TPL=54.5 μg/L；取滞留率为69%，则可求出目标RP=60.12 t/a；再根据TPLIN、TPL与公式(4)求得WOUT= 5亿m3，即于桥水库出库生态水量不小于5亿m3. 同时分析表1发现，当WIN< 5亿m3时，TPL较高，不利于于桥水库维持中贫营养状态. 因而，结合前文Tennant法计算得到的生态水量范围，建议生态入库水量大于5亿m3，此时生态出库水量约为入库水量的80%.

此外，2017年于桥水库前置库建成并试运行，对降低入库TP负荷进一步起到了积极作用. 该前置库位于水库上游果河，总占地面积22 km2，工程设施包括坝高4 m橡胶坝1座、渠道、闸涵桥等工程，种植挺水植物(芦苇和香蒲)和沉水植物(菹草、狐尾藻和微齿眼子菜). 通过前置库单元内的生态系统进行水体净化，库水面积18 km2，设计年净化规模5.82亿m3，年削减总磷58 t(约占入库负荷的48%)，从而减少营养负荷进入于桥水库. 依据前置库的设计总磷消减率，当WIN为5亿～9亿m3范围时，TPL最大值为41.8 μg/L(WIN=9亿m3)，前置库有效降低了水库TP浓度. 因此，本文建议的生态水量配合前置库运行将有效保持于桥水库处于Ⅲ类水状态.

表1 年入库水量WIN下生态水量估算*

*TPLIN通过图3的拟合线计算得出；τ的计算公式为：τ=-0.44 lnWIN+1.2481 (R2=0.67)；TPL根据公式(3)计算得出；RP=TPLIN×69%;TPLOUT根据公式(4)计算得出.

3 结论

本文采用Vollenweider模型分析了入库水量WIN、入库总磷负荷TPLIN和磷滞留量RP南水北调通水前后变化规律. 研究表明，引水结构变化后，于桥水库入库水量减少，但并未使TP入库负荷减少，进而磷滞留量增加. 库区TP滞留量增加也将增加水库内源负荷污染风险. 因而，笔者建议通过调节出入库生态水量以维持磷滞留在南水北调通水前水平，防治于桥水库富营养化趋势. 建议入库生态水量大于5亿m3，出库生态水量约为入库水量的80%，配合新建前置库工程运行将有效保障水库水源功能.

此外，本文在揭示入库水量和TP负荷、水力停留时间和TP滞留量机制关系方面做出贡献. 通过WIN与TPLIN相关关系发现，TPLIN随WIN呈先降低后升高的趋势，WIN较小时，对TP的稀释作用导致TPLIN降低；WIN较大时，引水携带TP负荷作用造成两者显著正相关. 探讨水力停留时间τ对RP影响机制发现，引水期τ(τ<0.35年)和RP之间存在显著的正相关关系，随着τ增加，正相关性下降，引水期τ越长越有利于TP滞留；但同时考虑引水期和非引水期τ时，τ与RP呈负相关，非引水期或引水量较少时RP主要受TPLIN的影响，此时TPLIN随着WIN降低而降低，导致RP较低. 研究成果为于桥水库运行管理提供了科学依据.