水温是湖泊水库水质监测中一项非常重要的理化指标，它与水中溶解氧(DO)、化学需氧量(COD)等水质参数密切相关，也是湖泊水库分层研究中的重要参数，水温的变化还直接影响湖泊水库的富营养化和生态变化^[1-4].绝大多数深水型水库一般在夏、冬季节会普遍发生水温分层，夏季典型的正向水温分层结构从上到下依次为变温层、温跃层和等温层.季节性水温分层是分层水库水质内源污染的主要诱因之一^[5-6].温跃层会阻碍上下水层的物质和能量交换，表层水体中较高浓度的溶解氧很难穿过温跃层而传递到底部水体，底部水体和底泥因各种化学、生物等作用而耗氧，二者共同导致底部水体溶解氧浓度逐渐降低，随温跃层厚度增加和持续时间延长，溶解氧浓度的降低程度增大，底泥中氮、磷、有机物等污染物释放量越大，内源污染越严重^[6-8].当溶解氧浓度低于2mg/L时，水库底泥中的无机和有机污染物会大量释放，形成水库内源污染，在冬春交替季节“翻库”而污染整个水库水体^[9].因此，研究温跃层的形成及其影响因素对进一步研究破坏分层水库水质的演变和内源污染控制技术具有重要的理论和现实意义.

尽管国内外学者开展了大量水库水温预测模型和应用研究^[10-12]，但很多模拟软件均不同程度地存在速度场和能量场耦合难的问题，多数只是注重模拟结果与实际数据的对比，缺乏对影响温跃层形成主要因素的研究与分析，尤其是对于富营养化分层水库.目前关于温跃层的研究主要针对海洋和浅水湖泊的水温观测，如陈希等^[13]、张文静等^[14]、赵保仁^[15]现场考察了净辐射通量、海潮、洋流对海洋温跃层形成的影响，认为净辐射通量的季节性变化是造成海洋跃层产生季节性变化的主要原因之一；由洋流引起的海潮将海面处高温水体带入下层，加速温跃层的形成.Brett等^[16]针对水深较浅的Mirror湖和Swan湖，提出了包含太阳辐射、风速、水深、水体密度和比热容等因素的水温分层参数，但未深入研究温跃层的形成机制.Lap等^[17]根据对浅水湖泊的观测结果，认为水温分层的形成与消失主要是由于风的影响和太阳辐射.Tuan等^[18]发现风速、风向对于浅水湖泊中热量的混合作用会产生一定影响.佘丰宁等^[19]对太湖温跃层形成因素的研究认为，气温的变化会使得水温呈线性变化，而大气湿度的增加会使水温呈指数变化；目前，仅在个别深水型湖泊和水库进行了温跃层的现场观测和初步分析^[20-22].Zhang等^[20]根据千岛湖多年水温监测数据，分析了该湖水温及温跃层的月度和季节性变化规律，建立了温跃层深度、厚度和强度与表层水温和水体透明度的经验模型，认为气候变化明显影响该湖水温结构.学者们对抚仙湖水温分层进行了系统的长期观测^[23]，对深水湖泊中温跃层的相关特征作出了定义，并根据实测数据，分析了温跃层的形成与演变过程.然而，关于深水型分层湖泊水库温跃层的形成特性及影响因素的系统研究尚未见报道.

本文应用商业Fluent软件^[24]并编写水体热交换的用户自定义函数(UDF)^[25]对其进行了二次开发，以平均水深80m的西安金盆水库水体为对象，系统研究不同表层短波辐射强度及藻类浓度条件下的水温结构，旨在研究和分析气象及环境因子对温跃层形成的影响，为控制分层水库水质污染提供参考.

1 数值模拟方法 1.1 水动力学控制方程

采用Fluent软件自带的可压缩传热紊流模型，包括流体连续性方程、动量方程、能量方程和紊流方程，其中紊流方程采用标准的k-ε紊流模型.在浮力传热紊流中，密度和粘度均随温度和压力变化，但压力的影响可忽略，根据不同温度下水的密度和粘度数据，建立水的密度和粘度与水温的函数关系，采用C语言编写UDF，编译并加载到Fluent中.

1.2 水体热交换模型

水库水体内的热交换按性质通常分为对流、传导、辐射和蒸发4种，按影响的区域不同又可分为表面传热和穿透传热两种.大气长波辐射、蒸发和传导均属于表面传热，表面传热一般影响深度不超过2m；太阳短波辐射则属于穿透传热，会作用于更深的水体.实际上，表面0~2m处水温受表面传热和表面短波辐射影响.水温分层是表面传热和水下短波辐射传热共同作用的结果，二者影响范围的有限性导致水温在一定水深处出现跃降，自上而下形成水温变幅较小的变温层、温度梯度较大的温跃层和水温相对恒定的等温层.

根据Hodges的水体热交换模型^[26]，水体表面总的热交换量(Q_S)为长波辐射量(Q_R)、水面蒸发量(Q_W)、水面显热对流量(Q_H)之和，即：

$ {Q_{\rm{S}}} = {Q_{\rm{R}}} + {Q_{\rm{W}}} + {Q_{\rm{H}}} $

(1)

长波辐射量(Q_R)的计算公式为：

$ {Q_{\rm{R}}} = - {\varepsilon _{({\rm{water}})}} \cdot \delta \cdot {T^4}_{({\rm{water}})} + {\varepsilon _{({\rm{air}})}} \cdot \delta \cdot {T}_{({\rm{air}})}^4(1 - {R_{{\rm{t}}({\rm{W}})}})(1 + 0.17{C_{({\rm{cloud}})}}^2) $

(2)

$ {\varepsilon _{{\rm{(air)}}}} = 0.920 \times {10^{ - 5}}{T^2}_{({\rm{air}})} $

(3)

式中，T_(water)、ε_(water)、ε_(air)、R_t(w)(1/C)、δ、T_(air)、C_(cloud)分别为水面0m处温度(K)、水的辐射系数(取0.96)、大气的辐射系数、水面长波辐射反射率(取0.3)、波兹尔曼常数(取5.669×10^－8W/(m²C⁴))、水面上2m处气温(K)、云层覆盖率(取0.17).

水面蒸发量(Q_W)的计算公式为：

$ {Q_{\rm{W}}} = L \cdot {C_{\rm{W}}} \cdot {u_{({\rm{wind}})}} \cdot {\rho _{{\rm{(air)}}}}\frac{{0.622 \times 10\left\{ {\frac{{0.7859 + 0.03477T}}{{1 + 0.00412T}} - 3} \right\}}}{P}({R_{\rm{h}}} - 1) $

(4)

式中，L、C_W、ρ_(air)、u_(wind)、P、R_h分别为蒸发潜热(取2.5×10⁶J/kg)、蒸发传热系数(取0.0014kg/m³)、空气密度(取1.225kg/m³)、水面处的风速(m/s)、一个标准大气压(pa)、相对湿度.

水面显热对流量(Q_H)的计算公式为：

$ {Q_{\rm{H}}} = {C_{\rm{H}}} \cdot {c_{{\rm{p(air)}}}} \cdot {\rho _{{\rm{(air)}}}} \cdot {u_{{\rm{(wind)}}}}({{T}_{{\rm{(air\_r)}}}} - {{T}_{{\rm{(water)}}}}) $

(5)

式中，C_H、c_p(air)、T_{(air_r)}分别为显热传导体积系数(取0.0014)、水的定压比热(取1003J/kg C)、水面上空气的干球温度(K)，其余同上.

水面太阳净短波辐射量(Q_sw(S))和水下短波辐射量(Q_sw(Z))的计算公式为：

$ {Q_{{\rm{sw}}}}(\mathit{S}) = {\mathit{Q}_{{\rm{(sw\_surface)}}}}(1-0.65{\mathit{C}_{{\rm{(cloud)}}}}^2)(1-{\mathit{R}_{{\rm{t(sw)}}}}) $

(6)

$ {Q_{{\rm{sw}}}}(\mathit{Z}) = {\mathit{Q}_{{\rm{sw}}}}(\mathit{S}){\rm{exp\{-}}{\eta _{\rm{e}}}{\rm{(}}\mathit{S-z}{\rm{)\} }} $

(7)

式中，R_t(sw)、Q_{(sw_surface)}、Q_sw(S)、η_e分别为无量纲的短波辐射表面反射率、通过水面的总短波辐射量(W/m²)、透过水面的净短波辐射(W/m²)、衰减系数(m^－1).

1.3 金盆水库及计算网格、初始条件与边界条件

研究对象为西安金盆水库，该水库平均水深70~100m，日供水量80×10⁴m³，为西安市主要饮用水水源.近年来，该水库水体富营养化程度不断加剧，其中氮、磷营养盐浓度阶段性超标，有机质含量上升，藻类季节性高发等水质问题尤为突出，处于中富营养化状态，7 8月表层藻类浓度可高达2×10⁷~3×10⁷cells/L^[4].以坝前主库区的断面水体为研究对象，所研究水域半径约为500m.水库地形采用中海达RTK系统测量.为长期监测金盆水库水质，在主库区设置5个监测点，本文所用数据均取自引水塔和主坝之间的S1点^[4]，该点约位于图 1所示断面宽度的1/3处(近右侧).水温采用美国HACH Hydro-Lab DS5型多参数水质分析仪现场监测，垂向测点间距一般为5~10m；藻类现场取样(1L)，采用鲁哥试剂(15ml)固定后，带回实验室进行显微镜计数分析，垂向测点分别位于水面下方0.5、5、30m，频率1~3次，本文所用藻类浓度为水深0.5m处值.采用GAMBIT构建水域断面的实际地形网格，针对水深月际变幅较小的3、5、7、9和12月，水深约为80~88m，图 1为水深80m(7月)坝前断面的二维网格，共有12835个节点(nodes)，12600个单元(cells).

该河道型水库上游河道长近40km，根据水库坝前特征计算出流量和断面面积，进口采用速度入口边界(velocity-inlet)，计算的断面平均流速约为0.0002m/s，由此计算左上侧等效进水口和右下侧出水口断面的平均流速.水库边壁和底部都采用无滑移固体壁面条件(Wall)，设置为绝热墙，其当量粗糙高度取0.0003m.

根据金盆水库底部水温资料(断面中心处)，将水库水体初始温度设为279.15K，开启能量方程(energy equation)，采用k-ε紊流模型进行计算，通过水体的自然热交换实现分层.表面和水下传热过程采用C语言编写的UDF写入，表面传热项由DEFINE_PROFILE(name，thread，i)宏将方程通过Heat Flux项写入，水下传热通过DEFINE_SOURCE(name，cell，thread，dS，eqn)宏写入水体源项(source terms).

1.4 水库温跃层影响因子及模拟参数

在分层水库地形和水文条件相对固定的条件下，水库温跃层的形成是水库内热力和动力过程共同作用的结果^{[23, 27]}，受各种气象和环境等因素的影响，如风力混合、水体清澈度(藻华堆积)、大气温度、云量、太阳辐射强度等^[27]；水库的运行调度(如进出水流)、原位水质控制设施的运行等也会影响温跃层的形成与破坏；自然水文事件(如夏季汛期暴雨径流和/或冬春季融雪径流等)也在一定程度上影响温跃层的变化^[28].本文以中富营养化的金盆水库为案例，因该地区春、秋季一般均为20d左右，故未将气温作为主要因素进行研究，而是重点研究气象因子(如太阳辐射强度、风速等)以及水质因子(藻类浓度)对温跃层形成的影响，并利用非暴雨径流期间实测水温数据进行模型验证.根据西安市气象局和相关文献^[29-30]，2009年西安金盆水库库区相关气象和水文参数如表 1所示.

根据改进的水体热交换模型，将以上气象条件参数通过UDF导入，直接进行基于Fluent的数值计算，时间步长选取10s；由于实际水库中水体形成分层过程较慢，每个条件下运行30d左右(即259200s)，通过温度云图观察温跃层的形成和变化特性.模型垂向深度间隔约为2m，由于表面传热仅仅作用于水面下1~2 m深度处，故忽略可能产生的能量跃变问题.

2 结果与分析 2.1 模型验证及温跃层演变

根据王银珠等的研究结果^[21]，将垂向的温度梯度大于0.2℃/m的水层定义为温跃层.利用表 1所示基础资料并对照在金盆水库库区断面中实测的垂向水温分布数据^[30]，变化不同的短波辐射衰减系数进行水温模拟，对比不同条件下模拟的水温结构和实测的水温结构，对短波辐射衰减系数进行参数率定，发现不同季节水温结构的模拟结果与实测结果吻合良好，对弱分层、较强分层和混合阶段的水温结果进行对比，发现其相关系数均在0.91以上，说明了这种模拟方法的准确性(图 2).

3月水库处于弱分层期，温跃层约位于水深3~8m处，温跃层内温差约为3℃.7月水库处于强分层期，温跃层约位于水深3~19m处，温跃层内温差约为20℃，水深2.5m内为上部变温层，水深19m以下为等温层.9月水库水温分层不断减弱，温跃层约下潜至水深30~70m处，温跃层内温差约为8℃，水深30m内为变温层.12月，冬季气温较低，上下层水体温差微小(约1.9℃)，基本处于等温状态，水温分层结构消失，不存在温跃层.对照2009年实测数据，温跃层的形成过程更为完整，冬、春季节整个库区水温基本相同，不存在温跃层；夏、秋季节水面和水底开始存在一定温差，极易形成水温分层，尤其是7 8月，水面和底部水体温差最大，水温分层最为强烈，温跃层最为稳定、厚度最大.

对比国内关于深水湖泊的开创性研究成果^[21]，金盆水库温跃层的垂向演变特性与平均水深87m的抚仙湖情况大体类似，只是该湖地处中亚热带半湿润季风气候区，水温常年高于12℃，不同于金盆水库.抚仙湖1、2月份处于完全混合期，其余月份水温自上而下均可分为变温层、温跃层和等温层，但温跃层深度、厚度以及层内温度梯度因季节而异，增温期温跃层位置较浅、厚度较大、温度梯度较小，升温期温跃层位置下移、厚度变薄、温度梯度增大.此外，该湖水面宽广，温跃层特征也在水平方向略有差异，与水面较小的金盆水库不同.

2.2 风对温跃层形成的影响

风在一定水深范围内、一定条件下影响着水面的混合强度和水体能量的耗损，一般风驱动的混合作用使得上部水体水温趋于均匀，在中部形成稳定度很大的温跃层，风对浅水湖库和海洋水温分层影响较大^{[20, 27]}.具体而言，风对水温的影响至少体现在3方面：1)影响表面传热；2)影响波浪和湖流，造成上下层水层交换；3)改变水柱中光衰减物质组成，造成短波辐射衰减系数变化，影响短波辐射传热量.

根据表 1所示金盆水库气象及水文资料，应用根据Hodges的水体热交换模型，分别计算3、5、7、9、12月水库表面各项传热量和总传热量，并绘制于图 3.根据传热量计算结果，如果不计表面的短波辐射传热量，长波辐射传热量约占总传热量的95 % ~100 %，说明长波辐射传热量是决定表面热量散失的绝对主导因素，而水面蒸发和热对流引起的传热量对水面处热量散失的影响极其微小.根据段誉等的研究^[31]，当水面总传热量为负值时，风的影响较小，基本可以忽略；另外，金盆水库实际气象条件(表 1)表明，风速较小，水面蒸发和热对流引起的传热量与风速相关，因此风对金盆水库水温结构的影响可忽略不计.如果考虑表面的短波辐射传热量，则在春、夏季节，水面总传热量为正值，水面将处于加热状态，表层水的密度将随水温增加而减少，温跃层将会变厚；在秋、冬季节，水面总传热量为负值，水面处于散热状态，温跃层变稀薄，等温层加深.在本研究条件下，坝前主库区水面面积和风速均较小，可以认为风的作用对金盆水库温跃层形成的影响较小，不同于浅水环境^[32].

2.3 短波衰减系数对温跃层形成的影响

在本模拟中，水库水面净热交换量为负值，说明水库温跃层的形成主要取决于水体受纳的短波辐射传热量，影响短波辐射传热量的主要参数是水面短波辐射强度和水体短波衰减系数.先以表层和底部水温差别最大的7月为例，此时水体分层较强烈，利用表 1所示的实际气象资料，分别取短波辐射衰减系数为0.5、2.0、3.0和4.0m^-1进行模拟，运行一个月后的水温结构如图 4所示.

随着短波辐射衰减系数从0.5增加到4.0，在表层水深5m范围内，虽然不同条件下的水体水温总体基本在26℃附近；但表层水温仍有略降低的趋势，水深5m处水温约从26.5℃降低到25℃；随着水深的增加，水温随短波辐射衰减系数的增加而降低的趋势更为明显，水深20m处水温约从18℃降低到9.5℃.依据公式(7)，短波辐射传热量随水库水深的增加而呈指数衰减，衰减的快慢主要与短波辐射衰减系数有关，短波辐射衰减系数越大，单位水深内藻类及其他悬浮物等吸收的热量越大，水体实际受纳的短波辐射热量越小，表现为水温增加值越小，太阳短波辐射的穿透深度越小.模拟结果与理论分析相一致.

根据图 4所示的水温结构，短波辐射衰减系数对温跃层的形成有直接影响.如以水温7~9℃作为等温层水温，则当短波辐射衰减系数分别为0.5、2.0、3.0和4.0m^-1时，温跃层位置分别约为水下8~43m、5~35m、4~30 m、3~21m，跃层内温度梯度分别为0.50、0.57、0.63、0.92℃/m.利用表 1所示其他4个月份的基础气象和水文资料，变化短波辐射衰减系数进行水温模拟，得到的温跃层特性也较为相似，即随短波辐射衰减系数的降低，温跃层位置下移、温跃层厚度降低，温跃层内温度梯度逐渐增加.

根据不同短波辐射衰减系数条件下模拟的垂向水温分布，对比实测的垂向水温分布数据^[18]，得出针对3、5、7、9和12月的最适宜短波辐射衰减系数分别为0.5、1.5、0.5、2.0和0.5m^-1.国外对湖泊水库水体短波辐射衰减系数的参考范围较大，从清澈的贫营养湖泊的0.2m^－1到浑浊的富营养湖泊的4.0m^－1[27]，本文校验的短波辐射衰减系数值亦与其相近.将不同月份经过校验的短波辐射衰减系数和实测藻类浓度共同绘于图 5，可以看出短波辐射衰减系数大小与藻类浓度高低的对应关系良好，较好地反映出中富营养化水库不同季节藻类生长情况与短波辐射衰减状况的关系.水库水的浊度很低，非汛期一般为10NTU左右，无机颗粒对短波辐射衰减系数的影响可以忽略，藻类成为影响短波辐射衰减系数的主要物质，故藻类浓度的增加会影响短波辐射的水下穿透量，从而造成短波辐射衰减系数降低，影响温跃层位置.水质监测结果显示，藻类浓度随水深增加，表现出先增加至水深5m处的峰值，然后降低，并在15m左右约降为峰值的1 %.尽管藻类浓度在垂向分布不均，但由于短波辐射热量呈指数型衰减，在夏季净短波辐射较强的时节，短波辐射能量的消耗主要集中在表层水体，故垂向采用统一的短波辐射衰减系数对水温预测结果的影响较小.

2.4 短波辐射强度对温跃层形成的影响

根据表 1所示水库7月的气象和水文资料，取实测短波辐射强度197.4W/m²和经过校验的短波辐射衰减系数0.5m^－1，变化水面短波辐射强度100、400、800W/m²进行水温模拟，虽然所取得的短波辐射强度变化范围较大，超出实测范围，但利用极端条件的水温模拟数据，可更全面、直观地评判短波辐射强度对温跃层的影响趋势，为后续短波辐射强度的相关影响研究提供参考.图 6表示水面短波辐射强度与上下层水体温差的关系.当水面短波辐射强度为100、400W/m²时，上下层水体温差随时间的推移而逐渐增加，但增幅减小，20d后趋于稳定，说明在次模拟条件下，水体热交换达到相对平衡.当水面短波辐射强度为800W/m²时，上下层水体温差也随时间的推移而增加，但增幅并未减小，20d后水体温差仍在继续增加.随着短波辐射传热过程的延续，水面及以下水体吸收的热量增加，导致表层水温增加，上下层温差随之增加.根据公式(7)，在短波辐射衰减系数相同的情况下，短波辐射在水下的衰减速率相同，并呈指数型衰减，可以看出在夏季净短波辐射较强的时节，短波辐射能量在穿透水下40m左右时，能量消耗达90 %，水体受纳的短波辐射传热量理论上主要与水面短波辐射强度有关，从而上下层水体温差随短波辐射强度的增加而增加.

图 6表明，水库可能存在一个水面短波辐射强度极值Q_S∈[400, 500].在较低的水面短波辐射强度(≤400W/m²)下，随着时间的推移，水面短波辐射强度越高，水面处吸收的热量越多，水面处水温增加；由于短波辐射强度在水深方向呈指数关系衰减，水体垂向各层吸收的能量向下递减，当短波辐射强度小于水库辐射强度极值时，表层水温因不断吸收短波辐射而增加，但能量来不及作用到水底，从而上下层水体温差逐渐增大，即上下层水体温差与水面短波辐射强度呈正相关.在较高的水面短波辐射强度(＞400W/m²)时，当短波辐射强度高于水库短波辐射强度极值时，由于水体吸收的热量较多，会影响到底部水体温度，使底部水温也增加，在一定的传热时间内，导致上下层水体温差减小；而随着传热过程的继续，受较高的水面短波辐射强度影响，表层水温增幅会高于底部水温增幅，导致后期上下层水温差增大.上下层水体温差随水面短波辐射强度的增加而呈现先增大后减小的趋势，在太湖地区也曾有所观测：在太湖地区，当短波辐射强度增大到700W/m²时，太阳辐射强度会使得垂向水体受到均匀加热，温差降低，产生短波辐射极值^[32].

对比本文初步得出的水库短波辐射强度极值与太湖地区^[31]报道的值，可以发现水库水温分层的形成与水面短波辐射的关系较为密切，上述变化过程应该还会与水库或湖泊的水深及不同纬度带上的气象与地质要素等有关，具体关系仍需要进一步研究与讨论；但是现有结果表明，水面短波辐射强度对水温结构和温跃层的形成有重要影响，当水面短波辐射强度超过极值时，很难形成稳定的温跃层.表 1资料显示，金盆水库地区短波辐射强度不高(约200W/m²)，水深较大(约80m)，故在季节交替时易形成稳定的温跃层，不同于浅水环境的分层情况.

3 结论

1) 通过模型参数率定，模拟研究了水温结构与温跃层随季节的变化，温跃层从春季气温回升逐步形成，在夏季形成稳定分层、温跃层厚度最大，秋季随着气温下降而逐步下潜，至冬季再次消失.模拟结果对水库水质演变和内源污染控制研究具有重要的参考价值和指导意义.

2) 在本模拟条件下，水库水面长波辐射、蒸发和传导传热量总和为负值，长波辐射量约占表面总传热量的95 % ~100 %，风基本不影响金盆水库温跃层的形成过程；短波辐射是影响温跃层形成的主要因素.

3) 经校验的短波辐射衰减系数值与实测的藻类浓度存在良好的正相关性，二者的年际变化相互一致；随短波辐射衰减系数的降低，温跃层位置下移、厚度增加，层内温度梯度减小.

4) 水面短波辐射强度对温跃层形成的影响较为复杂，不同类型水库具有不同的极限短波辐射强度，温跃层内温差随水面短波辐射强度的增加呈现先增加后减小的变化趋势；但水面短波辐射强度过高时，难以达到热交换平衡，对形成稳定温跃层有重要指示作用.