磨盘山水库温度场的数值模拟

伍悦滨，徐 莹，田 禹，张海龙

(1.哈尔滨工业大学城市水资源与水环境国家重点实验室，哈尔滨150090，ybwu@hit.edu.cn; 2.哈尔滨工业大学土木工程学院，哈尔滨150090;3.哈尔滨供水工程有限责任公司，哈尔滨150076)

水库蓄水带来了防洪、发电、灌溉，供水等综合效益，同时也引起了河流水文、泥沙、生态等各方面的环境影响，随着社会经济的发展，大量工业废水和城市污水未经处理排入河流、水库，超过了水体的自净能力，河流，水库污染日益严重，影响着人类的生存和发展.水温是水质的一个重要影响因素，因此进行水库水温分析和预测，是水库环境评价、规划、治理和管理的基础.

磨盘山水库位于拉林河干流上游五常市沙河子乡沈家营村的上游1.8 km附近，拉林河是松花江右岸的一级支流，流域总面积19200 km2.除供水外，磨盘山水库兼有防洪、灌溉及其他综合利用效益.水库的建成将改变库区及大坝下游河段的水文和水环境状况，进而对水库下游河段的灌溉和生活用水产生重大影响.

1 数学模型

EFDC(Environmental Fluid Dynamics Code)是由美国国家环保署资助开发，用于模拟河流、湖泊、水库、海湾、湿地和河口等水环境系统中流场、温度场、浓度场及生态过程的三维数值计算模型［1－2］.对水库水温进行模拟的数学模型主要采用Mellor等［3－5］提出的控制方程.

(1)连续方程

(2)曲线正交坐标中x向和y向动量方程

(3)温度

式中:u和v分别是沿着正交曲线水平方向x向和y向的流速分量;mx和my是平面坐标转换因子，m=mxmy.

其中，z为无量纲垂向坐标，w*为物理垂向速度;水深H=h+ζ，h代表实际垂向的底高程;ζ代表水位;动量方程中f表示科里奥利(Coriolis)参数; Av为垂向紊动粘性系数;Qu和Qv表示动量方程中的源汇项;θ为水温;在温度方程中的Qθ项是外源输入输出项.连续方程中的QSS和QSW表示河床与水体之间的泥沙及水的交换通量，其正值表示通量由河床入水体中，反之为负值.三角函数δ(0)表示这类通量进入水体的最底层.

2 数值计算及分析

磨盘山水库为典型的河道型水库，干流回水长约13 km.地形总趋势为东南高、西北低，呈不对称的“U”字型，地面坡降约为3.8%.磨盘山水库坝高45 m，正常蓄水位318 m，相应的水库总库容355.92×106m3，调节库容为292.95×106m3.汛期限制水位为 317.00 m，重复利用库容为27.27×106m3.设计洪水位为318.77 m，相应的设计总库容为375.98×106m3.校核洪水位为320.89 m，相应的校核总库容为443.02×106m3.

2.1 网格划分

库区地形资料是进行水库模拟计算的基础，首先根据水库地形高程图，按照模型要求划分计算网格.计算中时间步长为10 s，运行周期为1 d，即86 400 s，模型运行时间设定为30个周期.控制方程的数值解法为有限体积法和有限差分法的结合，对于模拟对象，沿流动方向采用正交曲线网格离散，而在垂直方向上则采用分层方法进行离散，进而引入了层厚的概念.根据磨盘山水库的地理特点，将其水平方向上划分1 052个矩形网格，每个单元格大小为85 m×65 m.考虑温度分层效应，垂直方向上分为8层，占水深比例分别为12.5%，12.5%，20%，20%，20%，5%，5%，5%.水库概化及计算区域单元格见图1.

图1 水库概化及计算区域单元格

2.2 初始条件

拉林河流域内现有沈家营、五常等13处水文站，根据磨盘山水库工程地点和资料情况，选用的主要分析测站是拉林河的沈家营测站.测站观测项目有降水、蒸发、气温、气压、日照、温度和风速等以及磨盘山库区的日气温变化.初始水位水深根据地形图得到，初始水温根据水库管理区测量资料以及《水利水电工程水文计算规范》(SDJ214－83)(试行)条例中给出.地面坡降为3.8%，初始计算流速设为零.

2.3 边界条件

根据库区管理站测量的2007年入库流量和出库流量作为上、下游边界条件，入库流量过程线见图2，出库流量过程线见图3.

图2 入流流量过程线

图3 出流流量过程线

2.4 模型验证

通过对磨盘山水库枯水期、平水期和丰水期的模拟，得到了水库出流温度变化曲线，并与实际出流温度变化进行了比较，结果如图4～6所示.模拟结果与实测值基本吻合，说明计算模型及方法的可信性.

图4 2007年4月坝前水温实测与模拟值对比

图5 2007年6月坝前水温实测与模拟值对比

图6 2007年9月坝前水温实测与模拟值对比

2.5 模拟结果及分析

图7 坝前断面的温度分层

本文分别模拟了2007年枯水期(四月)、平水期(六月)、丰水期(九月)磨盘山水库的水温变化规律.坝前水温在枯水期、平水期、丰水期的温度场计算结果分别见图7.水温在4月份为升温期，入流水温和表面热通量逐渐增加，上层水体水温增长迅速;在4月末，表层水温已达到10℃，但是靠近底部的水温几乎没有变化，仍保持低于6℃的低温，在水面附近出现了较大的温度梯度，形成了较明显的水温分层.6月份为高温期，入流水温和气温都达到最高.在表面很薄的水体中，由于水气热交换形成了表层温跃层［6－7］，表层水温保持在15℃作用，底部的水温缓慢上升了2～3℃，但仍然保持低于8℃的温度，在整个垂向断面上的温差达到了接近10℃的水平，在这一时期，水体上部分保持稳定的温度分层现象.9月份为降温期，气温和入流水温逐渐下降，水体向大气散失热量，冷水下沉使表层温跃层消失，水库在整体垂直方向形成水温分层，并开始出现层与层之间的垂向掺混，即秋季翻滚现象［8－10］.

水库的这种分层结构对水质具有重要的影响.对于磨盘山水库而言，由于其温度分层作用，不同深度的水质状况将有所不同.在枯水期，水库水温混合基本均匀，导致水中的营养物质沿垂向均匀混合.又加之冬季水中生化反应缓慢，消耗较少的营养物质，至春季时沉积了大量营养物质，为藻类快速生长提供了物质基础.平水期处于夏季，水库上部分水体水温分层较为明显.此时温跃层在垂直方向上具有“屏障”作用［11－13］，水中的营养物质在垂直方向上的交换量被大大的削减.在丰水期，汛期来临，洪水挟带着冲起的河底淤沙，涌进水库，使大量的亚磷酸盐进入水流.水库水温的稳定分层被破坏后，“屏障”消失，水库底层的沉淀物和有害物质会翻滚到表层，对水体造成污染.

图8是枯水期、平水期和丰水期在水库输水洞断面的模拟水温分布剖面图.磨盘山水库的输水洞位于右岸，进口位于右岸坝轴线上游约250 m处.水体自输水洞进入长直管线，向净水厂输送原水.由模拟结果可知，水库在枯水期和平水期沿深度方向上水温分布稳定，水温较低，可以直接选取水面以下10 m的水库水作为饮用水水源;在丰水期，水温结构发生变化.随着气温降低，表层水温也下降，相应的水密度增大，导致表层水开始下沉，底层水上升，在整个深度范围内表层和底层发生了混合，这种秋季翻滚现象可能会引发水质问题［14－15］.

图8 输水洞断面水温分布

3 结论

1)通过实测值与模拟值的对比验证，应用EFDC建立的水库水动力学模型，能够很好的模拟水库中水温的时空分布.因此基于该模型进行水体中流场、温度场、浓度场及生态过程的数值模拟是进行水库环境保护、管理及运行的有力工具.

2)磨盘山水库属于分层型水库，水温出现季节性分层现象，库表分层现象最为明显，至库底分层减弱，库底水温变化平缓.枯水期垂向分层现象不明显，平水期形成表层温跃层，对于垂向传质具有一定程度的屏障作用，丰水期出现秋季翻滚，对水质影响较大.

3)水库取水口深度应该随季节调整.在枯水期和平水期沿深度方向上水温分布稳定，大部分水体水温较低，水深10 m以下均可布置取水口;而在丰水期出现的秋季翻滚现象，形成水体垂向的自循环流动，可能引发水质问题.因此，丰水期至冬季时段，取水口位置应有所调整，并建议净水厂强化水处理工艺.

［1］HAMRICKJ M.A Three-dimensional environmental fluid dynamics computer code:Theoretical and computational aspects［M］.Gloucester，Massachusetts:Virgina Institute of Marine Science，the College of William and Mary，1992:1－63.

［2］李云生，刘伟江，吴悦颖，等.美国水质模型研究进展综述［J］.水利水电技术，2006(2):68－73.

［3］MELLOR G L，BLUMBERG A F.Modeling vertical and horizontal diffusivities with the sigma coordinate system［J］.Mon Wea Rev，1985，20(5):851－875.

［4］MELLOR G L，YAMADA T.Development of turbulence closure model for geophysical fluid problems［J］.Rev Geophysical Space Phys，1982，20(6):851－875.

［5］ARAKAWA A，LAMBV R.Computational design of the basic dynamical processes of the UCLA general circulation model［M］//CHANG J.Methods in Computational Physics.New York:Academic Press，1977:173－265.

［6］梁瑞峰，李嘉，李克峰，等.立面二维水库水温模型的并行研究［J］.四川大学学报，2009，41(1):30－33.

［7］张士杰，彭文启.二维水库水温结构及影响因素研究［J］.水利学报，2009，40(10):1254－1259.

［8］OUGOLNITSKY G A，USOV A B.Water quality control in watercourses［J］.Water Resources，2003，30(2): 226－232.

［9］申满斌，陈永灿，刘昭伟，等.岸边排放污染物浓度场三维浑水水质模型研究［J］.水力发电学报，2005，4(3):93－98.

［10］郭磊，高学平，张晨，等.北大港水库水质模拟及分析［J］.长江流域资源与环境，2007，16(1):11－16.

［11］WOOL T A，DAVIE S R，HUGO N R.Development of three-dimensional hydrodynamic and water quality models to support TMDL decision process for the Neuse River Estuary，North Carolina［J］.Journal of Water Resources Planning and Management，2003，129(4):295－306.

［12］JIN K R，HAMRICK J H，TISDALE T.Application of a three-dimensional hydrodynamic model for lake okeechobee［J］. Journal of Hydranlic Engineering ASCE，2001，126:13－14.

［13］陈景秋，赵万星，季振刚.重庆两江汇流水动力模型［J］.水动力学研究与进展，2005，12:829－835.

［14］JOHNSON B H，KIM K W，HEATH R E，et al.Validation of three-dimensional hydrodynamic model of chesapeake bay［J］.J Hyd Engrg，2003，119:2－20.

［15］李锦绣，杜斌，孙以三.水动力条件对富营养化影响规律探讨［J］.水利水电技术，2005(5):23－28