密云水库垂向水温模型研究陈永灿 张宝旭 李玉梁(清华大学水利水电系)摘 要 本文根据密云水库自身的特点，采用垂向一维水温模型对密云水库水温进行预测.模型中充分考虑水面热交换、入流、出流、热扩散、热对流等影响因素，利用 1991 年实测资料对模型参数 Dz、η 进行校定，根据校定的水质模型对 1992 年的水温进行预测，得出全年水温随水深分层变化并得到相应实测分层资料的良好验证.此模型可作为库区其它水质参数分析预测的基础.关键词 水库，水温模型，水质参数密云水库位于北京的东北方向，是一个山谷型的半封闭型的水库，全库最大蓄水量 43.75 亿 m3，相应水面面积 188km2，最大水深 43.5m，有明显的热分层现象.水库热分层现象是水库水质模型的重要特征之一.大多数水库都有热分层现象，只是有强有弱.根据水库水温分层情况的强弱，可分三种类型，即：1、混合型，2、稳定分层型，3、介于两者之间的过渡型.稳定分层型分层情况最强，而混合型分层情况最差.判断水库水温分层类型一般采用 α 指标法［1］：α=入库总库容/总库容 当 α＜10 时，为稳定分层型；当 α＞20 时，为完全混合型.此处可用 α 指标法判断密云水库的分层状况.已知密云水库 1991 年入库总流量为 8.3499 亿 m3，总库容为 25.0 亿 m3，计算得：α=8.3499/25.0=0.334<10由此可知，密云水库属于稳定的分层型水库.除用 α 指标法外，还可利用实测资料对密云水库的分层情况进行判断.密云水库管理处 1991 年、1992 年对密云水库部分月份进行了温度分层的观测，这些资料基本上反映了水库水温的分布特性.就整个水库而言，水温的变化主要反映在垂向，尤其在夏季可达 10℃以上，而水平向水温相差不大.可以认为，密云水库库区水温主要随时间和深度变化.因此可选择垂向一维水温模型对水库水温进行预测.利用数学模型来研究水库水温分层是 60 年代末从美国开始的.美国水资源工程公司(WRE,Inc)的 Orlob 和 Selna 及 MIT 的 Huber 和 Harlemen 分别独立地提出了 WRE 模型和 MIT 模型［1，2］，实现了水库的垂向温度分层模拟.这两种模型都得到过实测资料的良好验证，现在应用仍很广泛.在我国，水库水温模型正处于发展阶段.工程中大多采用经验公式法，这种方法虽具有简单实用的优点，但不能反映短时段的变化，并缺乏理论依据.本文将利用垂向一维模型对密云水库水温进行预测.1 密云水库水温数学模型1.1 模型的建立 垂向一维模型的基本思想是把水体划分为如图 1 所示的一系列水平薄层，忽略水平薄层中的温度变化，假设热交换只沿垂向进行，水平面温度均匀分布，可对水平薄层写出其质量和能量守恒方程.1.1.1 质量守恒方程在水体中任取一单元，其质量守恒方程可表示为：=Qv,j-1 - Qv,j + Qi,j-1 - Q0,j + Qa (1)式中 Vj 为第 j 层的体积，除表层外其它各层的 Qa 为表面降雨及蒸发的净值，除表层外其它各层 Qa=0；Qv,j 及 Qv,j-1 为第 j 层及第 j-1 层的垂向流量；Qi,j 及 Q0,j 为第 j 层水平向的进流和出流流量.1.1.2 热平衡方程对水体内各单元，其热平衡必须考虑水平向进、出流的热量，垂向移流的热量和离散的热量，太阳短波辐射热量以及由这些因素引起的单元体内热量的变化.对第 j 单元其平衡方程为： (2)图 1 水库垂向分层示意图式中 Cp 为比热；ρ 为密度；Hsz=Hs*e-ηz 是水深 z 处的太阳短波辐射热量.η为短波在水中的衰减系数，一般为 0.1m-1 到 0.2m-2 之间；A 为垂向的过流面积；Dz 为垂向混合系数.式(2)适用于除表面单元外的其他单元，即 J＜N 的各层.对 J=N 的表面单元要考虑水面交换的问题，式中的太阳短波辐射项需用水面热交换量与单元 N 底部的太阳短波辐射热量之差来代替.即表层应按下式计算：(3)1.2 密云水库水面热交换热通量的计算要得出水库的水温结构，必须先计算水面热交换热通量.一般来说，水面热交换包括辐射、蒸发、传导三个方面，具体地，通过水面而进入水体的热通量φm 为：φn=φsn+φan-φbr-φe-φc (4)式中 φsn——太阳短波辐射与水面对短波辐射的反射；φan——大气长波辐射及水面对长波的反射；φbr——水体的长波返回辐射；φe——净蒸发；φc——热传导.(1)太阳短波辐射与水面对短波辐射的反射 φsnφsn 一般可引用现场或邻近主要气象台站所测得的太阳辐射量值，扣除水面反射部分后求得.φsn=φs*(1-γs)， (5)式中：φs 总辐射量；γs 代表反射率，参考其它水库［6，7］，密云水库取 0.1.日照总辐射 φsn 经过水面反射后，部分进入水库水体，其中一半左右在水面被吸收，剩余部分按指数衰减进入水体深处.计算公式如下：φy=(1-B)*φsn*exp(-η*z)式中：B 为水面吸收率，参照其它水库，密云水库取 B=0.5；η 为衰减系数.(2)大气长波辐射及水面对长波的反射 φanφan 值须根据气温及云量观测间接计算，公式为：φsn=(1-γa)*σ*εa*(273+Ta)4 (w/m2) (6)式中 Ta 是水面以上 2m 处的气温，单位为℃；γa 为长波反射率，取0.03；σ 是 Stefan-Boltzman 常数，为 5.67×10-8W/m2*K4；εa 是大气的发射率，它和温度有密切关系.晴天的大气发射率 εac 可用 Idso 及 Jackson 公式算出：εac=1-0.261*exp(-0.74×10-4Ta2);多云天的大气发射率，可用 Bolz 公式算出：εa=εac*(1+KC2)，式中 C 是云层覆盖比例.K 是云层高度确定的，美国田纳西工程管理局推荐其平均值 0.17.(3)水体的长波返回辐射 φbr水体吸收的大气长波辐射能量会向大气进行返回辐射，是水体热损失的很重要的一部分.当把水体作为绝对黑体看待时，φbr 可由 Stefan-Boltzman 四次方定律来计算：φbr=σ*εw*(273+Ts)4 (7)式中 Ts 为水面温度，单位为℃；εw 为水面的长波发射率，它是一个常数.由于水体并非绝对黑体，εw 略小于 1 为 0.97.(4)水面净蒸发热通量 φe估算蒸发的方法很多，其中大多数是经验性的.蒸发的热转换公式通常为：φe=f(w)(es-ea) (W/m2) (8)式中 es 为相应于水面温度 Ts 紧靠水面的空气的饱和蒸发压力：es=exp［20.85-5278/(Ts+273.3)］ (mmHg)ea 为水面上空气的蒸气压力，单位 mmHg.f(w)为用风速表示的风函数.一般来说风函数包括了自由对流及强迫对流两者对蒸发的影响.可按下式计算风函数：f(w)=9.2+0.46W2z (W/m2*mmHg)式中 Wz 为水面以上 10m 的风速，单位为 m/s.(5)热传导通量 φc当气温不等于水温时，水汽交界面上会通过传导进行热交换，热传导率正比于两种介质之间的温度差.类比于蒸发热损失计算式，有：φc=0.47f(w)*(Ts-Ta) (W/m2) (9)现将密云水库 91 年、92 年各种热交换数值列入表 2 中：表 2 密云水库水面热交换计算(W/m2)1991年 1 月 2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月φsn 95.94126.1168.4199.9 236.6 234.4 197.8184.7172.2135.1 97.9 83.5φan 212.8220.5243.2276.7306.6 325.1 333.8327.0303.4273.8241.8218.9φbr 305.5305.5328.5338.1370.6 426.2 445.5435.5442.5402.6368.0305.5φe 41.3634.9043.9231.8537.65 100.1 77.2 55.83150.1106.289.0233.57φc 28.2913.573.69 -40.5-39.52-2.6864.88 3.32 35.2 35.8650.1215.891992年 1 月 2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月φsn 95.94126.1168.4199.9 236.6 234.4 197.8184.7172.2135.1 97.9 83.5φan 212.8220.5243.2276.7306.6 325.1 333.8327.0303.4273.8241.8218.9φbr 315.5311.8321.0340.0371.6 399.8 439.6439.6439.6399.8357.8333.3φe 51.3040.9434.9834.7039.57 38.33 61.6066.03142.0100.072.9062.47φc 41.5321.83-6.14-37.9-38.4 -27.0 0 6.65 32.6633.1838.9949.94由表 2 可以看出，φan、φbr 量值较大，φsn 次之.所以大气长波辐射及水体返回辐射对水库的水温值有较大的影响，另外由于日照辐射将有一部分按指数衰减规律进入水体深处，所以它对于水库水温结构有很大影响.1.3 水库入流、出流流速场的计算对于入流、出流的流速分布，因缺少实测资料，以往的研究均假定为高斯分布.据介绍，日本在水库中的观测表明，库内的流向与流速分布极不规则，取绝对平均后，近似于均匀分布［4］.当入流水温低于水体表层水温时，入流下沉进入水温相等的层面而形成入流层.如图 2 所示.图 2 水库入流、出流模拟入流层厚度为：(10)式中 q1 为入库单宽流量:q1=Qin/Bi,Qin——入库流量,Bi——入水口宽度,入水口中心线高程处密度梯度 当水库有多条河流入流时，按上述方法逐个计算，然后叠加.由于水温(密度)分层的影响，出库水流只在水库中一定厚度的范围内流动.出库流动层的中心高程为出水口中心高程.对于出水层厚度δ0：(11)式中 q1 为出库单宽流量:q0=Qout/B0,Qout——出水流量,B0——出水口宽度,出水口中心线高程处密度梯度 当水库有多个出口时，按上述方法逐个计算，然后叠加.1.4 垂向混合系数的确定垂向混合系数包括垂向紊动扩散和用垂向一维平均化的方法来描述物质运动引起的离散.垂向混合系数 Dz 是随时间、地点而变的，在进行温度计算前必须确定.据文献［7］，可根据密云水库的风速及水深，采用下式计算 Dz：D=Az×10-4+5×10-4*W*e-0.46y (12)式中：W——水面 10m 以上风速(m/s),y——水面以下深度(m),Az——待定系数.1.5 热对流水体表面热通量为正值时，表层处于升温状况，水温较以下各层为高，密度分层稳定.反之，在降温状况，表层水温则可能低于以下各层，形成不稳定状态，这时上下层发生热对流直到不稳定状态消失为止.图 3 水体热对流模拟水温分布示意图模型在各个时段完在前述计算之后，即应对所得温度分布进行检验.如发现存在不稳定状态，假设即刻发生热对流，将上层冷水与下层热水均匀掺混.计算时，沿深度向下逐层掺匀，直至掺匀后温度与该层原温度相等处为止.(13)热对流模拟完成后，所得水温分布即作为该时段末水温值.以此为根据，再进行下一时段的计算.1.6 边界条件的处理要求解此水温一维模型方程，必须先定出边界条件.水面边界条件已由方程(4)给出.定库底边界条件时，认为库底是绝热的，得：Zb。