大型跨流域调水工程对水源地水动力的影响研究

周红玉，刘 操，吴晓辉

(1.首都师范大学 资源环境与旅游学院，北京 100048；2.北京市水科学技术研究院,北京 100048)

0 引 言

南水北调工程是中国有史以来最大的跨流域调水工程，其中线工程将汉江上丹江口水库水由南向北调到北京市团城湖，再由团城湖通过京密引水渠输送至密云水库。密云水库为北京市唯一的地表水源地，南水北调来水后密云水库的水利联系由原来的由北向南来水变为南北同时通水，必然会对其水动力、水质及其水生态产生一定的影响，而对水动力的影响是最基本、最重要的部分。为确保密云水库水环境安全，本文就南水北调来水对密云水库水动力的影响进行模拟研究，为研究南水北调来水对密云水库水质及水生态的影响做准备。除了实际现场调研外，水文模型是研究该类问题较成熟的手段，而且方便成本低。水文模型对河流、湖泊、水库等的模拟已发展成为对水环境模拟、研究、预测与治理的较为成熟的研究手段。MIKE21模型采用有限体积计算方法，在国内外水环境研究中得到广泛的使用，其数值模拟的科学性得到世界公认[1-4]。

MIKE21模型利用二阶精度的有限差分方法对连续方程和动量方程求解，可以模拟各种作用力下忽略分层的二维自由表面流的水位和水流变化[5]。本研究搭建2015年MIKE21水动力模型(HD)模拟密云水库水位、面积、流速、流态等水文信息，通过与实际监测数据对比分析，经过不断调试，找到适合密云水库的模拟参数；通过2016年水文数据验证分析检验水动力模型的可靠性；分析模拟了代表性引水方式、出入流量、不同风场等不同条件影响下密云水库水动力的空间变化特征。本文通过构建密云水库水动力模型，研究南水北调来水后密云水库的水文变化情况，为下一步研究南水北调对密云水库水质影响奠定基础，为密云水库的管理和决策提供技术支持。

1 研究区概况

密云水库位于北京市密云县境内，是京津唐地区第一大水库，也是南水北调中线工程的重要调蓄水库，原具有防洪、灌溉等重要功能，现作为北京市地表饮用水源地。密云水库上游承接潮河和白河来水[6]，下游北京市第九水厂从中取水，流域面积为15 788 km2,占潮白河流域面积的88%[7]。密云水库的经纬度范围大致为北纬40°19′～41°31′，东经115°25′，～117°33′，气候类型为温带大陆性季风气候，冬夏季时间长，春秋季节短，降雨季节性比较明显，夏季多雨，其他季节少雨，年降水量变幅较大。密云水库最高蓄水面积188 km2，最大入库容量43.75 亿m3, 水深40～60 m，但由于近年来气候干旱，水库库容锐减，水库库容在10 亿m3左右，水位在135 m左右[7]。北京市缺水状况一直以来是北京市水环境的重要问题，为缓解缺水严重问题，南水北调实施后将丹江口水库水调到北京团城湖，再由京密引水渠输送到密云水库。

2 密云水库模型的搭建

密云水库有两大支流，一条是白河，源于河北省沽源县，流经赤城县、延庆县、怀柔区，流入密云水库；潮河起源于河北省丰宁县，自古北口入密云水库。模拟区域的边界设为进湖两处(白河和潮河)，白河来水以张家坟水文站每日实测流量控制，潮河来水以下会水文站每日实测流量控制。模拟区域的出口有一处，密云水库下游出水口是第九自来水厂(简称水九)，以水九水文站每日实测水位控制。

南水北调中线工程由丹江口水口引水至北京的团城湖，再由京密引水渠输送至密云水库，2014年年底正式通水，计划在2014年后平均每年为北京带来10～14 亿m3的水资源，平均每年有2 亿m3左右的水进入密云水库。在模型验证时将南水北调来水设置为进水口，以实际监测的每日流量控制。

2.1 MIKE21模型简介

MIKE21模型是MIKE模型中的核心基础模块，是丹麦水力学研究所(DHI)开发的系列软件之一。主要应用领域包括湖泊、河流、水库、海洋、波浪及泥沙等的模拟研究，主要应用于与水有关的工程领域及环境问题方面，包括河流水体等的水动力、水质模拟，水环境功能模拟，海湾波浪、潮汐等的模拟，水环境规划研究等。

MIKE21模型的控制方程如下：

连续方程：

(1)

动量方程：

(3)

式中：ε为水位高程；p、q分别是x、y方向上的流量；h为水深；t为时间；g为重力加速度；C为谢才系数；Ω为Coriol系数；f为风摩擦力；V为风速；Vx、Vy分别为风在x、y方向上的分量。

2.2 二维水动力模型计算条件设置

二维水动力模型需要输入的数据包括：密云水库库区网格数据(库区闭合水陆边界)、库区水深高程数据、开边界流量和水位时间序列数据、模型模拟参数、模拟初始水位及流量、模拟时间和时间步长的设定。经过多次调试得出适合密云水库的二维水动力模型。

2.2.1 模型网格设置

在密云水库水力联系中，潮河和白河从上游来水进入到密云水库，密云水库位于两河下游，系拦蓄两河河水而成，模型将两河入库处设置为开边界，第九水厂出水以点源的形式输出，其余为闭边界。模型采用非结构三角网格划分研究区，三角网格能较好的拟合研究区海陆边界，且能随意调整网格密度和网格大小。由密云水库地形图得到研究区海陆边界线数据(Land.xyz)以及水深数据(Water.xyz),然后将海陆边界线数据(Land.xyz)导入网格生成器Mesh Generator生成研究区网格并设置模型开边界，模型网格数为41 267个，三角网格的边长在100 m之内。然后将水深数据(Water.xyz)导入网格中并插值生成研究区地形图如图1。

图1 研究区地形Fig.1 Topographic map of the study area

2.2.2 水动力模型计算条件的设置

MIKE21 水动力模型需要设置的计算条件包括：模拟时间与时间步长、地形、克朗值、干湿边界、涡黏系数、降雨和蒸发数据、源和汇、边界条件以及初始条件等。

初始条件：初始条件设为2015年年初水面高程133.18 m，水平速度和垂直速度均为0 m/s。

时间序列：模型中需要添加的时间序列数据有白河和潮河两个水文站张家坟、下会每日实测流量数据，水九每日水位实测数据以及库区实测降雨和蒸发数据，时间序列数据由密云水库管理处提供的2015-01-01～2015-09-30的实测数据。两个边界的位置以及第九水厂的位置如图1所示。

模拟时间与时间步长(time step)：水动力模型的模拟时间为2015年1月1日8∶00至2015年9月30日8:00。模拟时间步长(time step)为100 s，时间步数为235 008步。

克朗值(CFL number)：网格分辨率、水深和时间步长决定了模型设置中的克朗值。克朗值在小于1的情况下模型才能正常运行，所以为了保证模型运行的稳定性将CFL值设置为0.8。

干湿边界(Flood and dry)：MIKE21 可以根据每个网格的水深情况调整计算条件，灵活的调用公式参与计算[8]。干湿水深是用来判断单元网格是否参与到模型的计算中来。在该模型中干水深为0.001m,淹没水深为0.05 m,湿水深为0.1 m。

涡黏系数(Eddy Viscosity)：选用Smagorinsky公式，取值0.28 m1/3/s。

底床摩擦力(Bed Resistance)：底床摩擦力可以选择谢才数(Chezy Number)、曼宁数(Manning Number)或者不选，这里选用曼宁系数(Manning Number),取值32 m1/3/s。

其他没有设置的参数采用默认值。

2.3 水动力模型计算结果分析与率定

利用MIKE 21FM水动力学模型对密云水库库区进行模拟，将白河库区的水位结果输出，白河库区位置如图1所示。将白河库区的实测水位与模拟水位进行比较率定，结果如图2所示。由图2可以得出白河库区的模拟结果与实际结果相比，相对误差在0.2%以内，绝对误差在-0.015～0.268 m范围内。从总体上看，模拟的效果比较好，模拟误差控制在模型计算要求的范围内。

图2 白河库区模拟水位与实测水位对比图Fig.2Comparison of simulated and measured water level in Baihe reservoir area

2.4 水动力模型的验证

为了验证水动力模型的可靠性，在率定好的水动力模型的基础上，用2016年1月1日至2016年11月8日水文数据对模型进行验证，同时将南水北调来水以点源的形式通入密云水库，采用每日实际监测数据。模型同样将白河库区水位输出，通过与实际监测数据对比，验证模型的可靠性。如图3所示，白河库区实测水位与模拟水位对比，绝对误差98%在-0.206～0.30 m范围内。通过验证，模拟水位与实际水位吻合较好，模型比较可靠。

图3 2016年白河库区水位模拟值与实测值对比Fig.3 Comparison of simulated and measured water level in Baihe reservoir area in 2016

3 南水北调来水对密云水库水文情势的影响分析

3.1 工况设计

通过对实际调水情况的分析，本模型在综合考虑当地气象条件，在不同的风场变化的条件下以及不同的引排水情况分析不同情景对密云水库水动力的影响。

密云水库地区常年主导风向为SW风，多年平均风速为2.24 m/s。综合考虑现有气象水文资料，本研究设计两种风场情况分别为2.24 m/s SW风和无风，来验证密云水库受不同风场影响下水动力的变化情况。

根据密云水库水力联系及调水情况可知，水体流动情况为北部上游白河和潮河来水与西南部南水北调来水汇兑在东南的第九水厂流出。根据实际调水情景可知，现有调水工程设计为调水量在10 m3/s左右,且冬季(12月-次年2月)不调水。为验证加大调水量以及不调水时不同的流量组合对密云水库水动力条件的影响，分别设计调水口进水量取0、10、20 m3/s。考虑到调水量增大后第九水厂的取水量相应的增加，水九的出水流量取1.5、3.0 m3/s。不同工况组合情况见表1。

表1 各计算工况方案Tab.2 Configuration of different water diversion schemes

表1中，工况1、2主要用于比较在相同引排水量及出入流组合方式下，密云水库在受不同风场的影响下流场的变化情况。工况1、3、4主要用于比较在受相同风场影响的情况下，不同的流量组合对密云水库水动力的影响变化情况。

3.2 水位及水面面积的影响

南水北调来水后对密云水库水动力最直接的影响就是水量的增加引起了水位的增高同时水面面积增大。图4为不同工况计算的水位随时间变化情况，可得，工况1与工况2的水位变化线完全重合，说明风对水位变化的影响不大。表2和3分别列出了不同工况下水位、水面面积的模拟情况以及与不调水相比水位、水面面积的变化情况。与不调水(工况3)相比，调水量为10 m3/s(工况1)时水位升高了2.4 m,变幅为1.7%，水面面积增大了7.43 km2，变幅为8.73%；调水量为20 m3/s(工况4)时水位升高了4.24 m,变幅为3.04%，水面面积增大12.24 km2，变幅为14.39%。

图4 水位模拟结果Fig.4Water level simulation result

工况水位/m最小值最大值差值水面面积/km2最小值最大值差值1135.76142.016.2568.7792.5023.732135.76142.026.2668.7792.5023.733135.76139.613.8568.7785.0716.34135.76143.858.0968.7797.3128.54

3.3 流场的影响

工况1是在风以及入流和出流的作用下形成的水体流动，而工况2是在无风的条件下，水库内水体的流动在入流和出流的作用下，由入流和出流的惯性驱动产生[9]。由图5和图6工况1和工况2计算的流场分布图可知，在调水同为10 m3/s的情况下，风场对湖区整体流场形态产生了较大的影响，在SW风的作用下整个湖区流场比无风时明显的多，而且在局部地区形成许多涡旋。工况2的流速大部分分布在0～0.006 m/s的范围内，工况1的流速大部分分布在0.0019～0.040 m/s的范围内，说明风场对流速的增大发挥了明显的作用。

表3 调水对密云水库水位、水面面积的影响Tab.3 Effects of water diversion on water level and surfacearea of Miyun reservoir

工况2、3、4为在不同引水流量的情况下流场的分布情况，工况2、4与工况3相比，除了在水九附近形成涡旋外，在调水口附近也形成一个涡旋，而且工况4形成的涡旋比工况2的大。从局部来看，调水口受调水的影响，流速明显增大，由原来的0.1 cm/s增加到0.46 cm/s左右。

根据上述不同工况的计算结果可得，南水北调来水后带动了水库水体的流动，水体流动性增大，减小了水体富营养化发生的可能性。同时水位升高，水面面积增大，外来物种的进入，可能会对密云水库原有水生植物、动物的生存环境带来一定的影响，接下来将做进一步的研究。

4 结 语

本文基于MIKE21FM模型建立了密云水库的二维水动力数值模型，对密云水库的水位变化过程进行了验证分析。模拟结果表明，水动力模拟结果较好，误差在计算要求的范围内。并对模型进行了检验，验证结果较好。在验证结果较好的基础上,分4种工况进行模拟分析，在不同的风场变化的条件下以及不同的引排水情况分析不同情景对密云水库水动力的影响，得出的结论如下：风场对水位及水面面积的变化无影响，不同的调水量引起水位升高和水面面积增大的幅度不同；风场的变化对流场的分布及流速的影响较大，在SW风的作用下整个湖区流场比无风时明显的多，而且在局部地区形成许多涡旋，流速明显增大；不同的调水量对流场的影响，与不调水相比随着调水量增大，在调水口附近形成一个涡旋，调水量增大涡旋也随着增大，而且调水口附近流速明显增大。

图5 工况1流场分布(2016年3月8日)Fig.5The distribution of a flow field for the No.1 scenario (March 8, 2016)

图6 工况2流场分布(2016年3月8日)Fig.6The distribution of a flow field for the No.2 scenario(March 8, 2016)

□

[1] COX B A.A review of currently available in-stream water-quality models and their applicability for simulating dissolved oxygen in lowland rivers[J].Sci Total Environ,2003,314/315/316(1):335-337.

[2] 谭炳卿，张国平.淮河流域水质管理模型[J].水资源保护，2001,(3)：15-18.

[3] KELLER V.Risk assessment of “down-the-drain”chemicals:search for a suitable model[J]. Sci Total Environ,2006,360(1/2/3):305-318.

[4] 陈沐松，宋兰兰.Mike21在垃圾发电厂温排水数值模拟中的应用[J].安徽农业科学，2011,39(15)：9 221-9 223.

[5] Danish Hydraulic Institute (DIH).MIKE21 Flow Model:Hydrodynamic Module Scientific Documentation[M]. Denmark :Danish Hydraulic Institute, 2007.

[6] 王晓燕.非点源污染过程机理与控制管理——以北京密云水库流域为例[M]. 北京：科学出版社，2011:22.

[7] 吴晓辉，吴 钢，潘轲旻，等.南水北调来水对密云水库水质和水生生物影响的预测分析[J].北京水务,2015,(6)：4-6.

[8] 李 娜，叶 闵.基于MIKE21的三峡库区涪陵段排污口COD扩散特征模拟及对下游水质的影响[J].华北水利水电学院学报,2011，32(1):128-131.

[9] 王 辉.大伙房水库流场影响因素的数值模拟研究[J].中国水运,2015,15(6)：76-77，81.