洪水过程下双峰寺水库水动力特征及污染物迁移规律研究

孙博伦，刘 彬*，张一帆，徐志恒，韩雨航，常一帆

(1.河北工程大学水利水电学院，河北 邯郸 056107；2.河北工程大学土木工程学院，河北 邯郸 056107)

2020年夏季，全国多地遭遇因强降雨引发的洪涝灾害，三峡大坝和葛洲坝等水利工程在防洪及水量调度方面起着至关重要的作用[1]。研究洪涝对水库的影响对于洪涝灾害防治、水利工程建设和水生态治理具有重要意义[2]。

一百多年前纳维、斯托克斯等通过研究得出N-S方程，为水动力模型的开始奠定了基础，随后在19世纪中后期Saint-Venant[3]提出的圣维南方程标志着水动力模拟的开始。发展至今，水动力模拟已被熟练应用于河流、湖泊和水库等相关课题中。伴随水动力模型的发展，水质模型自1925年Streeter和Phelps联合建立的S-P模型至今，已有百年历史[4]。通过阅读文献，在洪水与数值模型结合方面，专家、学者的研究多是关于洪水对库区水动力的影响，例如Stoke[5]通过将Saint-Venant方程组用于河道洪水中，得出求解完整Saint-venant方程组的数值模型，并根据其离散程度的不同，分成显式和隐式；Kamphuis[6]和Cunge等[7]选用的显式方法，对河道、水库洪水进行模拟，给出相关的表达式以及分析的结果；Mateo[8]使用水动力模型模拟水库作业对泰国湄南河流域洪水的影响，用来改善水库管理以减少洪水淹没；刘杰[9]通过数值模拟研究溃坝洪水对并联双桥墩冲击过程；姜治兵等[10]将坝体溃决过程与溃坝洪水演进进行耦合数值模拟等，而对库区的各污染物浓度随洪水过程的迁移规律还鲜有研究。本文在此基础上，通过整理洪水过程流量数据，并将水动力-水质模型耦合对双峰寺水库进行模拟，总结洪水过程下库区各污染物浓度的变化规律，为今后应对洪涝等自然灾害提供一定的理论基础。

双峰寺水库在承德市居民生活供水和发电、改善生态环境，防洪等方面承担着不可或缺的角色[11]，其在给生活带来了巨大便捷的同时，也会对生态系统造成一定的影响，使得生态系统平衡性遭到破坏[12]。研究洪水过程对库区的影响显得至关重要。本文针对双峰寺水库在典型洪水年上游来水水质的不同，河道径流变化较大，水质水量难以保证等问题，利用MIKE21构建水动力-水质耦合模型，模拟分析水库在典型洪水过程下水动力变化及污染物迁移规律的特征，为充分保障双峰寺水库发挥城市防洪、供水、生态改善及发电等功能提供理论依据，为水库水量-水质联合模拟及调控提供技术支撑，为有效控制水体污染提供基础保障。

1 研究区概况

滦河流域一级支流——武烈河流域，地处于华北平原，位置在东经117°42′～118°26′，北纬40°53′～41°42′，源于围场县道至沟，干流全长110 km，流域面积2 580 km2。流域上游有石洞子川、鹦鹉川、茅沟川、头沟川(兴隆河、鹦鹉河、茅沟河和玉带河)等4条支流，呈扇形分布，总集水面积占83%；下游区间6条旱河汇入干流，集水面积占17%。气候为温带大陆性季风气候，降水量年内分配极不均匀，多年平均年降水量为537.2 mm，其中近70%集中在汛期，多年平均水面蒸发量约为1 000 mm，年平均气温8.9℃[13]。

双峰寺水库位于4条支流汇合后的武烈河干流，控制流域面积2 303 km2，是防洪、供水、发电并兼顾改善生态环境等综合利用的大(2)型水利水电枢纽工程。水库防洪标准按百年一遇设计，总库容1.373亿m3，校核洪水位为395.11 m，设计洪水位为392.50 m，汛限水位387.0 m，水库正常蓄水位为389.0 m，死水位382 m；调洪库容0.738亿m3，兴利调节库容为0.451亿m3，死库容 0.335亿m3[13-14]。武烈河流域及水库具体位置见图1。

2 MIKE21模型简介

2.1 MIKE21 FM水动力模型

水动力模块作为MIKE21的基础及核心模块，同时也是运行水质、泥沙、粒子追踪等模块的基础，其主要模拟在不同外力作用下的水动力特性变化。控制方程遵循纳维-斯托克斯方程，并服从Boussinesq假设。

图1 武烈河流域及水库地形

2.2 ECO Lab水质模型

ECO Lab作为模拟水体水质、水体富营养化、重金属等运移变化规律的水生态工具，依据传统水质模块的基础发展起来，主要分为水质模块、富营养化模块和重金属模块。本文通过水质模块运算的模拟方程为积分法的耦合常微分方程，一般常用的有欧拉方程、龙格库塔四阶、龙格库塔五阶方法[15]。

欧拉方程：yn+1=yn+hf(xn,yn)

(1)

式中f(xn,yn)——动点(Xn，Xn+1)的平均速度。

龙格库塔法：

(2)

则四阶龙格库塔方程为：

(3)

3 应用研究

3.1 网格划分

通过双峰寺水库的地形，面积等库区特征来构建二维模型，模拟范围为2 303 km2。根据实测库底高程数据，建立区域模型网格，其中网格尺寸35～125 m，网格的平均面积为5 000 m2，共计单元个数4 905个，见图2a。在此基础上，将水深数据插入其中，生成双峰寺水库的水下模拟地形，见图2b。

a)计算网格

b)库底地形模拟

3.2 水动力模型参数设置

MIKE21水动力模型首先设置模拟开始时间、总时间步数以及主时间步长，其中主时间步长通常以秒为单位；其次对地形、克朗值、干湿边界、涡粘系数、降雨和蒸发数据、源和汇、边界条件以及初始条件进行参数设置[16]。水动力参数具体数值见表1，其中未在表中表明的其他参数，一般取模型默认值。

a)来水。根据本文设计内容，观测典型洪水年对水库水质的影响，故选取承德水文站实测1962年典型洪水过程作为工况来水情境，其实测洪峰流量为2 580 m3/s，见图3a。

b)出水。双峰寺水库工程主要由拦河坝和电站等组成，其中拦河坝由4部分：非溢流坝段、溢流坝段、底孔坝段和电站坝段组成。坝前水位为387 m汛限水位以下时出流为1.85 m3/s，以保证下游生态需水，超过汛限水位则按最大出流能力出流，见图3b。

表1 模型主要参数

a)入库洪水过程线

3.3 模型率定与验证

首先采取1994年实测数据进行参数率定，本文中，对模拟结果影响较大的为干湿边界和底床摩擦力，经调试，流量的率定结果相对误差见图4；再者对率定结果进行验证，采用1962年典型洪水过程，得出相对误差结果见图5。从整体上看，模拟的效果较好，相对误差在规定范围内(0～10%)，其中率定期的相对误差小于4.6%，验证期的相对误差小于2.2%，为水质模拟提供可靠的前提条件。

图4 率定的相对误差

图5 验证的相对误差

4 水质模型的构建

随着水资源日益短缺，对湖泊、河流及水库的要求已经不仅仅局限于水量，同时水质是否达标也愈来愈重要。在现代水文学研究中，对河流、湖泊、水库等水体进行水质模拟也是非常重要的一个分支。由于人们的生活、生产等将污水排放其中，在各种条件的作用下迁移转化，造成不可逆转的后果。通过MIKE21模型进行耦合模拟，从而得出水质状况的变化规律，及时对水质进行预测与分析。

4.1 水质模型参数确定

在水动力模型模拟的基础上，与其耦合ECO Lab水质模块对双峰寺水库进行模拟。

在确定状态变量所涉及的过程后，需要对各过程的参数进行分析率定，在参数的经验取值范围基础上，通过对双峰寺水库数据资料的详细分析，进行多次模拟率定，确定参数取值见表2。未在表中表明的其他参数，一般取模型默认值。

表2 状态变量过程参数

续表2 状态变量过程参数

水质模型的初始条件以2020年7月25日经过人工湿地净化后的入库初始浓度，库区水库以上窝铺采样点水质为依据，具体数值见表3。温度设为10°，盐度设为0，以及风速设置为常数(2 m/s)。

表3 洪水过程的入库和库区水质数据 单位：mg/L

4.2 水质模型模拟结果与分析

采用以上数据及参数对双峰寺水库进行水动力-水质耦合模拟，结果以点源的形式输出，选取T1(581 184.4，4 552 872.1),T2(582 874.3,4 551 137.9)，T3(582 255.8，4 548 442.8)3个坐标来对模拟结果进行分析，具体位置见图6。通过模型模拟T1、T2、T3三点在典型洪水过程中BOD、DO、NH4、NO3、PO4等5种变量随时间的变化，得出各物质浓度在洪水过程下的变化曲线，见图7—11。

图6 双峰寺水库T1、T2、T3三点具体位置(m)

图7 洪水过程中BOD浓度变化曲线

图8 洪水过程中DO浓度变化曲线

图9 洪水过程中NH4浓度变化曲线

图10 洪水过程中NO3浓度变化曲线

图11 洪水过程中PO4浓度变化曲线

由图7得出BOD浓度曲线初期迅速下降，浓度降低到1.0 mg/L，下降坡度变缓直至过程结束；其中在T1点下降至1.0 mg/L时，出现反弹现象，原因为洪水涌入过快导致水流紊乱形成漩涡。

由图8得出DO浓度曲线规律与BOD较为相似，但在T1、T2点上，浓度降低至7.8 mg/L时，受水流紊乱产生漩涡的影响较严重，致使浓度上下浮动。尤其在T1点的位置，洪水过程结束后出现浓度反弹现象。

由图9得出在洪水开始时， NH4浓度呈上升趋势，并且远离入库位置的坐标，曲线上升坡度越陡且时间越长，其中在T2点NH4浓度经过上下波动，重新上升到0.13 mg/L。

由图10得出在洪水初期，NO3浓度曲线短时间内未发生变化，而后浓度曲线急速上升至2.5 mg/L上下，并在浓度最高值附近上下波动且呈下降趋势。

由图11得出PO4浓度曲线规律与NH4的较为相似，其中T1点浓度变化甚微，T2点浓度波动较大，在洪水过程结束时PO4浓度达至0.036 mg/L，处于最高位。

5 结论

基于MIKE21模型建立了二维水动力数值模型，对双峰寺水库在洪水过程下的水位变化过程进行了率定与验证。模拟结果表明其误差在计算要求的范围内。在水动力模型的基础上构建水质模型，并通过对点源输出的结果进行对比分析。

a)在典型洪水过程的工况下，库区内流场变化导致污染物浓度波动剧烈，随洪水过程的发展，其中BOD和DO浓度变化曲线相似，洪水初期曲线迅速下降，降低到某一数值时，下降速度突然变缓直至过程结束。

b)NH4和PO4浓度变化曲线相似，在洪水初期三点浓度呈上升趋势，越远离入库位置的坐标，曲线上升坡度越陡且时间越长，浓度上升至最高点后迅速下降。

c)NO3浓度变化曲线较为特殊，洪水初期呈现短暂的平稳趋势，随后急速增加到最大值2.5 mg/L，并在其附近上下波动。

通过此次模拟可预测典型洪水情况下双峰寺水库的水质变化状况，为制定合理的双峰寺水库治理措施，起到关键性的作用。同时对双峰寺水库各污染物浓度进行分析，采取合理的治理对策对库区水质情况进行改善，具有一定的可行性和可操作性，为今后库区水质或流域改善治理起到一定的技术支撑作用。