石佛寺水库二维水动力及水质数值模拟与分析

胡 艳 海，周 林 飞

(沈阳农业大学 水利学院，辽宁 沈阳 110161)

长期以来，水库水质状况和库区水生态环境受到广大学者的关注，其中以水库水动力条件为基础进行水库污染物质运移消散过程模拟是主要的研究方式。国外关于水质模型的研究起步于20世纪60～80年代，计算机技术的迅速发展促进了双向耦合模型对水质模型变化所产生的影响，与此同时水质模型也伴随着水动力系数、空间变量、主要的外界影响因素(例如温度等)的引进而变得更加完善[1]。20世纪80年代后，多维数值模型求解伴随二维水质模型的成熟得以实现。国内关于水质模型的研究比较晚，20世纪80年代，由于官厅水库遭受严重的重金属污染事件，由此引起了多个领域的研究学者对水库水质问题的关注，并意识到利用水质数值模型解决水环境相关问题的重要性[2]。因此，众多学者针对不同类型的湖泊、湖泊中各种水质问题、富营养化问题以及不同的水质因素构建水质模型，最终实现水质模型可行性验证，同时为进一步实现模型的简化、促进数值模型的应用进展创造条件[3-5]。龚春生等[6]通过模拟玄武湖混合流的水质动态变化过程为底泥污染和浅水湖泊混合流水质变化模拟提供新途径。20世纪90年代，水质模型逐渐被广泛应用到水质优化方案设计方面[7]。21世纪后，水质模型逐渐向多维、多层次方向发展[8-11]。例如，杨具瑞[12]构建的垂直分层的二维水质模型，计勇[8]构建的水动力-水质-底泥模型，廖临毓[13]构建MIKE21水动力-水质耦合模型对COD、TN、TP三种污染物在3种调水方案下的流场和浓度场时空变化情况进行模拟，马宁[14]构建水动力-水质耦合模型模拟水库水动力和水质状况并获得水质改善方式。

石佛寺水库位于辽河干流，承担防洪、供水、生态等重任。因此，掌握库区水动力条件以及水质状况显得尤为重要。本次研究以石佛寺水库作为研究对象，通过总结前人经验构建水库的二维水动力-水质耦合模型，将石佛寺水库来水量作为水库水动力模型的主要影响因素，以水质模型浓度场变化和水动力模型流场变化作为运算基础，模拟运算了库区汛期和非汛期两种工况下溶解氧(DO)、氨氮(NH3-N)、总氮(TN)、五日生化需氧量(BOD5)、化学需氧量(CODCr)、总磷(TN)6种水质指标的运移消散状况。最后，根据分析模拟结果提出水环境综合治理的具体措施，用以指导水库运行，使辽河水流经水库后能够达到地表水Ⅲ类水质标准。

1 研究区域概况及水质监测

1.1 研究区域概况

石佛寺水库位于辽宁省沈阳市沈北新区黄家乡和法库县依牛堡乡，地理坐标为东经123.427°～123.520°，北纬42.144°～42.188°，距沈阳市区约47 km，如图1所示。石佛寺水库具有防洪、供水、生态三大功能，其中最主要的是防洪，为典型的河道型平原水库，是辽河干流上仅有的控制性水利工程。为改善区域生态环境，净化辽河上游来水，实现其生态功能，2009年起在石佛寺水库实施了生态建设工程，并在库区内辽河左岸种植水生植物芦苇、蒲草和荷花，同时进行生态蓄水，蓄水面积16.13 km2，蓄水位46.2 m。

1.2 水质监测

为了了解辽河水流经石佛寺水库后水质净化状况，在水库的入口、中间和出口布设水质监测点位，分别标记为采样点1、采样点2、采样点3，如图1所示。各监测点监测的水质指标为溶解氧(DO)、氨氮(NH3-N)、总氮(TN)、五日生化需氧量(BOD5)、化学需氧量(CODCr)、总磷(TP)等11个指标。DO利用多参数光谱水质分析仪测定，NH3-N采用水杨酸钠法测定，TN采用全自动间断化学分析仪测定，BOD5采用稀释与接种法测定，CODCr采用重铬酸盐法测定，TP采用硫酸钾—钼锑抗法测定。监测频率为每月2次，监测标准以GB 3838-2002《地表水环境质量标准》为依据。

监测时间为2009～2017年。石佛寺水库作为水源地，水标准执行地表水Ⅲ类。

2 基本原理

2.1 水动力数值模拟基本原理

MIKE21水动力模块是在二维数值求解浅水方程的基础上建立的[15-18]。二维水动力数值模拟分析基于以下4项基本假定展开计算：不可压缩性假定、静水压力假定、Boussinesq假定、Reynolds值均布假定。水动力模块控制方程应用沿水深积分的不可压缩的Navier-Stokes方程，其主要包括二维平面流的连续性方程、动量方程[19]，如式(1)～(3) 所示。本文选用单元中心的有限体积法作为数值解法，并采用MIKE21FM的非结构化三角形网格进行模型求解。

(1)

(2)

(3)

图1 石佛寺水库位置及水质监测点位置Fig.1 Location and water quality monitoring points of Shifosi Reservoir

(4)

(5)

(6)

2.2 水质模型基本原理

库湖水质控制方程是全面考虑污染物质水动力学、物理运输和对流扩散等过程并将质量守恒原理作为前提建立起来的。该控制方程的数学表达式如式(7)～(8) 所示[20]。

(7)

S=Qs(Cs-C)

(8)

二维水动力模型作为水质模型的基础，只有建立水动力模型并满足模型参数率定和模型验证，方能加载水质模型，其主要原理在于水动力模型能够提供水质模型水质浓度输出的流速等水动力条件。

3 石佛寺水库水动力及水质数值模拟

3.1 水动力数值模拟

水库水动力影响因素包括：进出水条件、风场、地形、库内水工构筑物以及水库的岸线形态等。其中，水库的地形和岸线形态作为相对固定的因素，来水量为主要因素。

3.1.1模拟范围

本文选用1∶5 000石佛寺水库电子版地形图作为地形资料，利用Auto CAD和南方CASS提取高程点并生成.xyz格式陆地land和水深water数据文件用于网格生成。采用Beijing_1954_GK_Zone_21N投影坐标，控制最大三角形面积20 000 m2、最小三角形角度30°以及最大节点个数150 000个，利用MIKE ZERO网格生成器生成网格，计算网格如图2所示。完成水深water文件插值并利用网格分析器进行网格分析，获得合理准确网格文件并由MIKE Animator生成模拟区域的三维地形，如图3所示。

图2 石佛寺水库计算网格Fig.2 Computational grid of Shifosi Reservoir

图3 石佛寺水库三维地形Fig.3 Three-dimensional topography of Shifosi Reservoir

3.1.2定解条件

(1) 初始条件。模拟时长为61 d，即2017年6月1日08：00：00至2017年7月31日08：00：00，时间步长为3 600 s。由于模型在短时间内计算趋于稳定，故在初始时刻t=0时，库区初始流速u和v设置为0，初始水位为δ=46.200 m。其他数据为实测值。

(2) 边界条件。根据MIKE软件前处理方法，上游边界为石佛寺水库入口处实测流量数据，下游边界为石佛寺水库出口处实测水位数据；水质模型利用岸壁法确定闭边界，即沿法线方向上的浓度值和流速值均为零；主要根据底床摩擦应力确定底床边界条件；采用可避免模型发散的“干湿水深判别法”确定移动边界，即水深大于湿水深0.1 m时按水域处理，水深小于干水深0.005 m时按陆域处理。

3.1.3模型的率定与验证

(1) 模型率定。选择2017年6，7月石佛寺水库入口处实测水位率定模型参数，率定结果为：曼宁系数45 m1/3/s(MIKE软件的曼宁系数与糙率值互为倒数)、涡黏系数0.28。

图4 出口合成流速模拟值与实测值对比Fig.4 Comparison of the simulated value and measured value of the exit flow velocity

图5 入口水位模拟值与实测值对比Fig.5 Comparison of the simulated value and measured value of the entrance water level

3.2 水质数值模拟

3.2.1水质模板选择与水质指标选取

本文选用改进的Eco Lab内置水质模板WQ with nutrients对石佛寺水库水质进行模拟。该水质模板包含状态变量7个，常量39个，作用力6个，辅助变量9个，过程21个。

综合考虑实际水质监测数据的超标情况和生态数值实验室(ECO Lab)内置模板中状态变量的可模拟情况确定水质指标为：DO，NH3-N，BOD5，TN，CODCr，TP，其中TN和TP作为衍生结果输出。

3.2.2定解条件与模型参数

(1) 初始条件。模拟范围、模拟网格系统及模拟时间同水动力模型。考虑模型整体一致，故在初始时刻t=0时，水质指标的初始浓度值设置为0。

(2) 边界条件。闭边界法线方向上设定流速和浓度值均为零；开边界，将水库进出口处水质指标实测浓度生成的时间序列文件作为上下游边界条件。

(3) 模型参数。建立生态模型时污染物的运移消散过程必然会发生，因此水质的扩散系数以及衰减系数的设置显得至关重要[21]。其中，扩散系数在各个方向上的上下限通常根据其与各网格速度分量之间的线性比例关系而确定；降解系数则通常根据原型观测法[22]获得，然后利用环境部门检测结果和相关研究成果进行修正和再次修正，最后确定各水质指标的降解系数如下：DO为0.075 14 mg/d；NH3-N为0.113 1 mg/d；BOD5为0.140 53 mg/d；TN为0.014 95 mg/d；CODCr为0.102 05 mg/d；TP为0.042 12 mg/d。

3.2.3水质模型率定与验证

(1) 模型率定。选择2017年6～7月石佛寺水库主槽区实测五日生化需氧量、化学需氧量、氨氮、总氮、总磷率定模型参数，率定结果为：横向扩散系数5 m2/s、纵向扩散系数5 m2/s，化学需氧量衰减系数0.078 5 mg/d，五日生化需氧量衰减系数0.108 1 mg/d，氨氮衰减系数0.087 mg/d，溶解氧衰减系数0.057 8 mg/d，总氮衰减系数0.011 5 mg/d，总磷衰减系数0.032 4 mg/d。

(2) 模型验证。基于模型参数的成功率定，将对模型进行验证，模型验证时选用2017年5，7月出口处实测生化需氧量、化学需氧量、氨氮、总氮、总磷数据验证模型，验证结果如图6所示。采用统计学中的偏差统计法，得到五日生化需氧量、化学需氧量、氨氮、溶解氧、总磷、总氮实测数据与计算模拟数据之间平均误差分别为4.04%，7.39%，5.10%，3.80%，3.05%，1.40%，均在可接受范围内，表明实测值与模拟计算值拟合程度较好，完全能够满足精度要求，所以本文所构建的水质模型是有效可行的。

4 模拟结果与分析

4.1 水动力模拟结果与分析

利用验证合理的模型进行模拟，将1号、2号和3号监测点作为特征点，用于反映石佛寺水库流速和流场特点，监测点布置如图1所示。

4.1.1模拟方案

本次水动力数值模拟共分为两种工况：即汛期和非汛期。根据石佛寺水库多年运行经验，6～9月为汛期，非汛期不考虑冰冻期，故非汛期仅模拟为10～11月和4～5月。汛期和非汛期初始流速均设置为0，初始水位均设置为46.200 m。汛期和非汛期上游边界为石佛寺水库入口处实测流量数据，下游边界为石佛寺水库出口处实测水位数据。模拟时长非汛期为2016年10月1日～11月30日、2017年4月1日～5月31日，共122 d；汛期为2017年6月1日～9月30日，共122 d。现选取汛期和非汛期水动力模型运行末期时典型流场分布情况进行分析。

4.1.2模拟结果分析

(1) 根据图7可知，汛期水位大部分区域介于46.280～46.256 m之间，非汛期水位大部分区域介于46.190～46.185 m之间，即汛期库区水位大部分仅略高于非汛期。原因是虽然汛期水量较大，但会通过下游闸门控制水库水位，使其保持在46.200 m左右，保障生态需水。

图6 各水质指标模拟值与实测值对比Fig.6 Comparison of the simulated value and measured value of water quality indicators

(2) 图8为汛期和非汛期不同水量条件下库区典型流场流速分布情况，结果表明两种工况下水库具有较好的流动状态且未出现死水区(流速≤0.3 cm/s时为死水区)和回水区。对比可以看出汛期的流速整体大于非汛期，说明不同的来水量水库内不同区域内的流速变化明显。两种工况下，流速较大区域均出现在上游入口区域和下游出口区域，这是因为过水断面面积小从而加速了水体运动；过水断面面积较大的区域，流速比较缓慢，占库区面积也大；流速较小区域出现在左右岸岸边，与岸边地形相关。总体来看，库区内水体流动性良好、水流不急，有利于水体的交换以及污染物的去除。

图7 水位分布(单位：m)Fig.7 Water level distribution

图8 流场流速分布(单位：m/s)Fig.8 Flow field velocity distribution

4.2 水质模拟结果与分析

来水量是石佛寺水库水动力的主要影响因素，因此选择非汛期和汛期两种不同情景对其进行水质模拟分析。具体模拟时间段为：非汛期为2016年10月1日至11月30日，2017年4月1日至5月31日，共122 d；汛期2017年6月1日至9月30日，共122 d。初始时刻水质指标浓度设置为0，采用水质指标浓度值生成的时间序列文件作为上下边界条件。

4.2.1模拟结果

本次数值模拟获得两种情景下水库DO，NH3-N，BOD5，TN，CODCr和TP在不同时间段的浓度分布情况。汛期情景下各污染物在模拟末期时分布图如图9所示，非汛期情景下各污染物在模拟末期时分布如图10所示。

图9 汛期各水质指标浓度分布(单位： mg/L)Fig.9 Distribution of water quality indicators in flood season

图10 非汛期各水质指标浓度分布(单位： mg/L)Fig.10 Distribution of water quality indicators in non-flood season

4.2.2分 析

汛期和非汛期除DO以外各污染物浓度整体分布情况为：入口>中间>出口，原因是石佛寺水库具有净化污染物质的作用，其中在非汛期尽管水流流速减缓，但是水体中的某些污染物质会伴随泥沙发生沉降，污染物质被存储在底泥中，进而实现净化水质的作用。DO为出口>中间>入口，流经水库复氧能力增强，水质得到改善。

在汛期NH3-N，TN，CODCr和BOD5浓度小于非汛期的浓度，原因一是汛期库区内水生植物长势茂盛，微生物具有很强的活性，有很好的水质净化作用；原因二是汛期有较大的水量，流速也较大，自然增强了水体的自净能力。在汛期DO的浓度大于非汛期的浓度，原因是水体复氧途径由原来的的大气复氧，增加了水生植物向水体中输送氧气和藻类等沉水植物的光合作用产生氧气，而且汛期流速快增加了水体的复氧能力。TP的浓度汛期和非汛期变化不明显。

水质从入口到出口逐渐得到净化，说明辽河水流经水库后水质得到了改善。但分析出口水质数据，仍有指标在一定时间段内达不到地表水Ⅲ类水质要求，因此，仍需加强库区内的污染综合控制，可以从以下3个方面进行考虑：

(1) 控制外源污染。进水污染物浓度为导致水库污染的主要因素，因此可通过控制面源污染和点源污染，尽量减少经降雨、下渗、径流过程形成的污染源汇入水库。在外源污染超过水库的净化能力时，可根据石佛寺水库水动力水质模拟结果进行调水，进一步改善水质。

(2) 减少内源污染。水库在非汛期水流流速会减缓，泥沙会携带一些污染物质沉降而储存在水库底泥中，尤其在水流相对缓慢的非汛期，为避免再次污染可以通过对库底底泥疏浚来控制内源的汇入。

(3) 种植水生植物。汛期的流场和水质指标浓度分布图叠加后得出，生长水生植物的区域水流受阻流速变缓，水生植物的净化作用使污染物浓度降低；在无水生植物的区域，水流流速较大的区域水质指标的浓度值相对较小些。因此，应该在水流流速较小的区域种植净化能力较强且阻水较弱的水生植物，在高流速区域种植阻水能力较强的水生植物。

5 结 论

(1) 石佛寺水库水动力数值模拟结果表明：来水量是水库水流状态的主要影响因素。汛期和非汛期两种工况下，汛期的流速整体大于非汛期，说明不同的来水量水库内不同区域内的流速变化明显；两种工况下水库均具有较好的流动状态，未出现死水区和回水区，且水流不急，有利于水体的交换以及污染物的去除。

(2) 石佛寺水库水质数值模拟结果表明：辽河水流经水库后逐级得到了净化，而且净化效果汛期优于非汛期。但在一定时间段内仍有指标达不到地表水Ⅲ类水质要求，需加强库区内的污染综合控制，具体措施如下：控制外源污染和调水净化；通过底泥疏浚减少内源污染；种植具有净化功能的水生植物，在水流流速较小的区域种植阻水较弱的水生植物，在高流速区域种植阻水能力较强的水生植物。