变化条件下官厅水库的水量平衡过程研究

2021-01-22秦政雷坤黄国鲜孙明东代丹程全国

环境工程技术学报 2021年1期

秦政，雷坤，黄国鲜,3，孙明东，代丹，程全国*

1.沈阳大学环境科学与工程学院 2.中国环境科学研究院 3.青海大学三江源生态与高原农牧业国家重点实验室

官厅水库是永定河流域最重要的大型控制性水库，水库上游社会经济的迅速发展使需水量日益增加，导致官厅水库入库水量从20世纪50年代的20.40亿m3a骤降到21世纪初的1.13亿m3a，近年来略有回升，达2.00亿m3a，入库水量锐减导致官厅水库持续多年低水位运行。官厅水库曾是北京市的重要供水水源地之一，由于水环境污染严重等原因，1997年被迫退出北京市生活饮用水体系，随着北京市水资源供需矛盾的日益突出，2007年其被重新启用为北京市的备用饮用水源地，但只保留了对北京市局部地区工业和农业的有限供水[1-2]。为缓解官厅水库及其下游水资源紧缺问题，自2003年起通过永定河上游水库陆续向官厅水库调水，在一定程度上改善了官厅水库蓄水量不足的状况，同时提升了水库水质；但官厅水库的蓄水水位仍较低，水库下泄流量不能满足连续下泄和下游河道最小生态流量的需求，使永定河下游的北京—天津段仍长期处于断流状态[3-4]。自2016年起从黄河万家寨调水，使流域获得一定的附加水资源[5]，永定河上游桑干河源头区—三家店区域都得到一定程度的补水。

为了更好地解决官厅水库及其下游水资源供需矛盾和可持续发展，同时维持官厅水库及其下游河道生态环境健康，迫切需要定量求解官厅水库水资源的收支过程和平衡，开展不同调水方案和生态流量下泄条件下的湖库水量过程和平衡问题的研究[6]。目前在官厅水库水量平衡方面的总量研究比较缺乏，已有研究主要基于统计学方法，围绕实测水量和降水等水文特性，进行长序列和年内变化的分析，或对官厅水库流域的降水、蒸发和径流等变量的变化特点进行研究[7-8]，缺乏利用平衡方程对官厅水库不同水文条件下水量平衡的整体分析。笔者采用2002—2016年的气象、水文和水库特征等数据，基于水库水量平衡方程建立模型，求解官厅水库蓄水量和水位的日均过程，在模型中有效考虑水库出入流、蒸发、降水和下渗等因素的具体过程和影响，分析官厅水库的水量平衡状况；定量分析官厅水库在不同收水量和生态流量下泄条件下水库蓄水和水位变化情况，并定量分析黄河万家寨调水能否满足官厅水库最小生态流量下泄的要求，以期为未来官厅水库水资源的合理调度和改善官厅水库下游河道的生态环境状况提供技术支撑和决策参考。

1 研究区概况与研究方法

1.1 研究区概况

官厅水库位于华北沿海平原向中温带干旱的内蒙古高原过渡的地带、河北省怀来县与北京市延庆区交界处(115°34′02"E～115°49′30"E，40°13′46"N～40°25′42"N)。官厅水库建于1954年，水库总库容为41.6亿m3，控制流域面积为43 402 km2，正常蓄水位为479 m，死水位为471.47 m。流入官厅水库的河流主要有洋河、桑干河和妫水河，其中洋河与桑干河在怀来县夹河村汇流进入永定河，永定河向东8 km注入官厅水库，而妫水河自东向西流经延庆区后注入官厅水库(图1)。永定河下泄水流经官厅水库的拦蓄和调节，沿永定河中下游进入北京和天津段。区域气候特征为夏短冬长，年平均气温为6～8 ℃，多年平均降水量为405 mm，降水主要集中在6—9月，占全年的64%～76%；库区多年平均蒸发量为1 200～1 400 mm，水面蒸发量大于降水量[9]。

图1 研究区及各控制站分布Fig.1 Study area and location of control stations

1.2 数据来源

根据模型计算要求，需要收集库区降水量与蒸发量，进、出流流量，调水量和坝前下渗过程的水文数据，水库水位-面积-库容曲线和水库特征水位等参数，数据来源见表1。

表1 水量平衡计算数据来源

由于官厅水库的底泥淤积较为严重，淤积量达6.51亿m3，而水库缺少近年的水位-面积曲线数据，因此基于1979—2014年的长序列遥感影像水体面积解译分析结果[13-15]，结合当时水库对应的水位，通过分析水库水位与对应面积的拟合曲线，求解得到水库的水位-面积曲线，结果如图2所示。由于所选用的分析资料大部分与本研究计算的时段(2002—2016年)重叠，因此得到的水位-面积曲线具有较好的代表性，可用于模拟时段内水库水量计算。官厅水库水位-库容关系采用1980年水位-库容曲线(图2)。

图2 官厅水库的水位-面积-库容曲线Fig.2 Water level-area-capacity curve of Guanting Reservoir

1.3 研究方法

根据质量守恒原理，利用官厅水库日时间尺度数据，建立水库水量平衡方程。公式如下：

(1)

式中：QIn(i)、QOut(j)分别为第i条入库河流、第j个水库出口的日平均流量，m3s；NI为入库河流条数；NO为水库出口数量；TsecDay为时间换算常量，表示1 d内包含的秒数，取值86 400；P为库区日降水量，mmd；E为库区日水面蒸发量，mmd，由直径20 cm蒸发皿的实测数据(E20)及折算系数(k0，取值0.60)(E=k0×E20)计算得出；A为不同水位下水库的水面面积，m2，可根据水库的水位-面积曲线求解得出；ΔV为水库蓄水量的日变化，ΔV=Vt+1-Vt(t为时间，d)；WMan为人工引水水量日均过程；Sp为库区渗透率，主要包括水库床面渗透率(Sp1)和水库坝体有无截渗墙条件下的水库渗漏损失率(Sp2+Sp3)。Sp的计算公式为：

(2)

式(2)第一项为根据土壤渗流达西公式计算得到的库底渗漏量。其中K为地基的渗透系数，md，由于官厅水库河道及库区底部细颗粒泥沙沉积厚度平均约为20 m，因此K取1×10-6～1×10-7cms[16]；Zs为水库水位，m；Zsa为库区地下水水位，m；ΔL为河床平均高程与地下水位的垂直距离[17]，m。式(2)第二项为铺设防渗墙的土坝的渗漏率。其中B1为防渗墙的宽度，m；K0为黏土地基渗透系数，md；T为透水地基厚度，m；Le为覆盖面材料的等效长度，m；Δh1为覆盖面材料的末端的水头损失，m。式(2)第三项为未铺设防渗墙的土坝的渗漏率。其中B2为未铺设防渗墙的土坝宽度，m；ΔH为上下游水位差，m；M1为坡度系数[18]。

计算模型主要采用Visual Studio和Inter Visual Fortran 2013组合平台下Fortran 95语言编写，计算时间步长为1 d，为了保证计算精度，在同一时间步长内采用有限次(最大200次)迭代算法进行水库水量平衡求解，水库接收的有效降水量、库面蒸发量和下渗量等主要与相邻时间步长内水库水面面积平均值有关，因此需要迭代计算的间接变量包括有效降水量、蒸发量和下渗量，最终停止迭代计算的决策变量是相邻迭代步内计算得出的水位、水面面积和水体体积差足够小。在迭代求解过程中，通过保持各出、入流流量不变条件下，迭代更新求解式(1)中每个时间步长内水库水面区域内的有效降水量、蒸发量和下渗损失量，最大循环200次，在迭代过程中发现绝大多数条件下一般经过3次迭代，相邻时段内计算的水位差绝对值就可以控制在1.0×10-8m以内。具体计算结果表明：迭代后相邻迭代步的水面面积差绝对值小于1.1×10-3m2，水体体积差绝对值在2.0×10-2m3以内，每日计算产生的误差较小。若不进行迭代计算，每步水位计算误差绝对值在9×10-5m以内，面积误差绝对值在1.2×104m2以内，体积误差绝对值在7.7×104m3以内，水位平均误差绝对值为3.4×10-6m，面积误差绝对值为34.7 m2，体积误差绝对值为186.7 m3；而迭代后的水位平均误差绝对值为5.2×10-9m，面积平均误差绝对值为5.4×10-2m2，体积平均误差绝对值为0.31 m3。说明从模型计算角度迭代计算是必要的，通过迭代计算能显著降低相邻时段内水位、面积和体积的计算误差，进而减小库区水面蒸发量和下渗量的计算误差。同时应该指出：库区水面蒸发量和下渗量的计算受到物理过程、计算公式类型和主要参数取值等多个环节的影响[19-20]，迭代求解只是减少计算误差的一个方面。

该模型计算所需输入的数据包括库区气温，水面蒸发量与降水量，水库进、出口各水文站数据，水库水位-面积-库容曲线数据，水库下泄的最小生态流量月过程数据，黄河万家寨调水月流量过程，水库各特征水位，水库历史实测水位和水库渗漏流量等。

1.4 平衡因子分析

1.4.1入库径流量

由于永定河流入官厅水库尾闾段的8号桥站没有长时间观测的水文站，所以选取位于洋河上的响水堡站和位于桑干河上的石匣里站的流量之和作为永定河入官厅水库径流量。洋河、桑干河、妫水河3条主要入汇河流多年入库径流量如图3所示。

图3 官厅水库各水文站2002—2016年入库径流量变化Fig.3 Variation of annual inflow water volume of each hydrological station of Guanting Reservoir in 2002-2016

1.4.2库区降水量及蒸发量

官厅水库库区水面的降水和蒸发对水库水量有一定程度的影响。利用延庆和沙城等气象站的日降水量和日蒸发量数据，采用反距离加权空间插值法求得官厅水库库面中心点的日降水量和日蒸发量。其中上述参考气象站的水面蒸发量采用直径为20 cm的小型蒸发皿实测蒸发数据乘以蒸发折算系数确定[21-23]。通过统计得到2002—2016年官厅水库多年水面降水量平均值为0.319亿m3a，蒸发量平均值为0.535亿m3a，降水量和蒸发量年变化如图4所示。

图4 官厅水库库区2002—2016年降水量与蒸发量变化Fig.4 Varation of annual precipitation and evaporation in Guanting Reservoir area in 2002-2016

1.4.3水库下渗量

官厅水库下渗主要分为水库库底下垫面的渗透、铺设防渗墙的大坝或两岸绕坝坝体的下渗、水库大坝与未铺设防渗墙的两岸绕坝坝体的下渗。经过多年运行，官厅水库库底已形成了一定厚度的淤积泥层，其主要由细沙和中沙组成，且分布较为均匀，渗透系数取1×10-7～1×10-6cms[24]；水库拦河坝为黏土心墙坝且坝前铺设水平的防渗盖，因此坝前的渗漏量较小；左右肩坝坝坡经过回填防渗土并与黏土防渗心墙连接等工程措施形成厚达5 m的防渗体，防渗效果较好；水库左右岸绕坝未铺设防渗体坝体，因此按照未铺设防渗墙条件下的渗漏公式计算下渗量。

2 结果与讨论

2.1 模拟结果分析

2.1.1水库水位模拟验证

在原有来流和调度作用下，模型计算的2002—2016年官厅水库水位变化如图5所示。由图5可知，水位的模拟值与实测值相关系数为0.933 3，二者吻合良好。但在水库低水位运行期，模拟值和实测值出现一定差距，主要是由于水库严重淤积，模型采用的水库水位-面积-库容曲线与实际曲线分布可能存在一定的差异。在水库水位较低条件下(如小于470 m)，官厅水库的永定河干道入汇口附近长期淤积形成的拦门沙将官厅水库分成2个相对隔离的子水库，此时实际参与调节的子水库(主要为永定河入汇子水库)的库容较小，在调度过程中，妫水河入汇的子库容不能自由地参与调度调节，因此整个水库的调节功能呈削弱状态，从而使水位的模拟值与实测值有一定差异。

图5 官厅水库水位模型验证结果(2002—2016年)Fig.5 Water level verification results of Guanting Reservoir in 2002-2016 by model

2.1.2下渗量模拟验证

根据模型计算的水库下渗量模拟值与1998年水库下渗量实测值对比如图6所示。由图6可以看出，模型计算的下渗量与水库水位相关曲线和1998年对应的实测水位与下渗量相关点群基本吻合，因此本研究采用的下渗量算法是合理的。1998年官厅水库下渗量随水位增加而增加，其变化幅度为0.050～0.065 m3s，量级远小于库区的来流量、库区降水量和蒸发量，因此可以判断下渗量对库区的水量平衡影响较小。

图6 官厅水库下渗量验证(2002—2016年)Fig.6 Verification of infiltration discharge of Guanting Reservoir in 2002-2016

2.1.3现状调度条件下的水库水量收支过程

现有调度条件下2002—2016年官厅水库水量平衡如表2所示。由表2可知，2002—2016年官厅水库多年平均入库径流量为1.58亿m3a，其中从桑干河、洋河和妫水河控制站直接入库的多年平均径流量为1.26亿m3a，占总入库量的80%；而库区多年平均降水量为0.32亿m3a，占总入库量的20%。多年平均出库水量为1.45亿m3a，其中多年平均蒸发量为0.54亿m3a，占总出库量的37.2%；水库多年平均下泄流量为0.9亿m3a，占总出库量的62.2%；多年平均下渗量为0.009亿m3a，仅占总出库量的0.6%，可见下渗量对水库整体运行水位影响较小。水库多年平均蓄变量为0.133亿m3a，总体上看，2002—2016年官厅水库属于正收支平衡状态，但2006年之前蓄变量多为负收支，之后多为正收支，水位呈逐渐升高态势。

表2 2002—2016年官厅水库水量平衡分析统计

虽然2002—2007年永定河流域的平均降水量为424.20 mm，属于多年偏枯年份，但由于上游对官厅水库的调水作用，入库径流量仍为水库来水的主要部分。为保证一定的下泄流量，官厅水库年均水位从2002年的473.48 m降至2007年的470.29 m，为多年来的最低水位。2007年后水库开始减小下泄流量，并且上游的调水水量减小直至2010年停止调水，但由于2008—2016年永定河流域多年平均降水量增至535.26 mm，流域降水对径流的影响具有一定的滞后且上游各水库的蓄水延长了滞后的时间，官厅水库入库径流量逐年增加，加之下泄流量的减少，官厅水库水位逐渐抬升。

将官厅水库各出入库项变化作为一个整体平衡过程(来水项为正，耗水项为负，下渗量因占比较小未列入)，得出水库平衡分析及各出入库项在平衡中的占比，结果如图7所示。由图7可知，官厅水库来水量构成比例并不稳定，年际变化较大。以洋河响水堡站为例，该站来水量在平衡中的占比最小，仅为7%，在水量充沛时可达31%。入库径流对水库蓄水量影响较大，各站径流量总和约占水库来水的38%～80%。在水库耗水项中，蒸发量的占比较为稳定，为12%～27%；水库下泄流量年际波动十分明显，为7%～51%，且整体呈递减趋势，下泄流量在水库耗水项中占主导地位。入库水量的不稳定导致水库下泄流量年际波动较大，这给下游水生态环境带来了一定风险，也给水库水资源调度管理带来一定的困难与挑战[25]。

图7 2002—2016年官厅水库水量平衡 Fig.7 Water balances of Guanting Reservoir in 2002-2016

2.2 水库下游生态流量计算

由于官厅水库上游来水量存在明显不足，加上水库水面蒸发量较大，水库下泄流量逐年减少，下泄流量不足将会直接影响下游的生态功能。目前官厅水库下游河段存在严重的水资源不足、河道断流和河道沙化等生态退化等问题，保障官厅水库下游生态流量对于维护永定河生态功能健康有着重要意义。

依据Tennant法计算官厅水库下游永定河河道生态需水量[26]。Tennant法以下泄流量与多年平均流量的之比为基准，将河流流量分为几个等级：最大允许极限、最佳范围、极好状态、良好状态、一般或较差状态、最小状态和极差状态[27-28](表3)。

表3 河道不同生态环境状况与流量对应情况

永定河流域降水年内主要集中于6—9月，占全年降水量的76%，10月—翌年5月降水量占全年的24%。根据1956—2009年系列月天然径流，计算出官厅水库下游(坝下)河道的生态需水量为1.84亿m3a，其中汛期为0.81亿m3a，非汛期为1.03亿m3a。基于该生态需水量数据，结合2002—2016年下游河道实测水量，计算官厅水库汛期与非汛期实际生态需水量亏缺情况，结果如表4所示。由表4可知，官厅水库下游河道生态水量亏缺较为严重，汛期多年平均亏缺量为0.50亿m3a，亏缺率为61.73%；非汛期多年平均亏缺量为0.60亿m3a，亏缺率为58.25%；官厅水库下游河道多年平均生态水量亏缺率为60.43%。可见，造成下游河道生态水量亏缺严重的原因之一是官厅水库下泄流量逐年减小。虽然多年来通过上游水库调水以及减小下泄流量，官厅水库蓄水水位有所升高，但下游河道的生态需水量没有得到相应的补给。因此在水库现有调度条件下难以满足下游河道的生态需水量，需要流域外的水资源调度来补充[29-30]。通过山西省万家寨引黄北线工程，连通桑干河与黄北干线向永定河补水(简称引黄调水)，是解决永定河生态水量亏缺的重要措施。引黄北线线路全长164 km，规划引水量为5.6亿m3a，可利用供水能力为2.64亿m3a[31]。

表4 2002—2016年官厅水库下游河道实测水量和生态需水量亏缺情况

2.3 不同调度方案下模型对水位的预测

利用官厅水量平衡模型，设置不同调水和生态流量下泄控制情景，模拟官厅水库水位和下游生态流量的关系。为了检测调水对官厅水库正常运行水位和水库下游生态需水量的影响，设置2种情景进行讨论：1)情景1，不进行调水，保持原有水库调度条件，只考虑满足下游最小生态流量限制；2)情景2，考虑每年3—5月从万家寨引黄调水1亿m3a，官厅水库收水0.5亿m3a，同时满足下游最小生态流量限制。针对以上2种情景进行模拟，同时与基于现状真实水量(不考虑补水调水，也不考虑下游生态需水量)的模拟结果进行比较，如图8所示。

图8 2种情景下官厅水库水位模拟结果Fig.8 Water level simulation results of two scenarios of Guanting Reservoir

由图8可知，情景1时，官厅水库为了满足下游最小生态流量限制，水位将会持续下降，年均水位从2002年的473.48 m跌至2014年的466.11 m，远小于官厅水库正常运行水位。尽管在降水量较大，且来流充沛的情况下还能保证一定的下泄流量，但在降水量减小条件下，官厅水库并不能满足下游最小生态流量的需求。情景2时，官厅水库水位保持稳定并在降水量较大的年份略有升高，即使在没有其他调水且在降水较少的年份也能保持一定的水位运行，既保证了官厅水库的水位维持，也保证了下游的最小生态流量。