一种基于热力耦合模型的结冰水库水温-冰情反演及预测方法与流程

本发明属于水利工程技术领域，涉及结冰水库水温-冰情反演及预测，尤其涉及一种基于热力耦合模型的结冰水库水温-冰情反演及预测方法。

背景技术：

冰是水资源重要的存在形式，普遍存在于寒冷或高纬度区域，在水文循环中担任着举足轻重的角色。在我国30°n以北的地区每年冬季都会有不同程度的冰情现象发生，约占全国国土面积的70％。众所周知，水库的修建使河流的水力学和水文学特征发生显著变化，与此同时，这些变化直接影响到水库及其上下游河道的热力过程。值得注意的是，随着社会的发展，水利水电工程的开发正逐渐向寒冷高海拔、高纬度的河源区推进，例如：雅鲁藏布江中上游、雅砻江上游、澜沧江上游等，在这些地方修建水库必然会面临冰情问题，如库区结冰和解冻延迟、冰下水温增加、水库下游产生不封冻河段等。

水库冰水的相互作用，主要体现在结冰前的水库蓄热量、水内冰与冰盖的热力生消、水内冰与冰盖的转化、冰下水温变化、“气-雪-冰-水”的热量传递过程等。在本领域内，已有的对水库冰水研究主要集中在经验公式和简化的数学模型上，大多只聚焦于冰盖的热力生消过程，如度日法只考虑关注局部地点冰的热力变化，水库垂向一维水温冰情耦合数学模型不能获取冰-水的空间差异，即使垂向二维ce-qual-w2在内的冰盖模型也均没有考虑冰水之间的作用与转化，使得无法对结冰水库冰、水热力全过程进行反演和有效模拟。

因此，开展高寒地区水库冰水热力的研究，实现高寒地区结冰水库的水温-冰情的有效反演和模拟，对于水库设计、运行管理、防凌减害、环境保护、水资源调度、工程安全等方面具有重要的科学意义和工程应用价值。

技术实现要素：

针对目前缺少对结冰水库冰、水热力过程进行有效反演和模拟，进而影响结冰水库水温-冰情进行有效反演的技术现状，本发明目的旨在提供基于热力耦合的结冰水库水温-冰情反演方法，基于所建立考虑水库蓄热影响、水体过冷却、水内冰生成与输移、初始冰盖生成、冰盖热力增长的宽度平均水库立面二维水动力、水温、冰情的耦合模型，反演、反演和评估结冰水库水温、冰情的热力时空变化与水库对水温、冰情的影响程度和范围。

为达到以上目的，本发明采用以下技术方案来实现。

本发明提供的基于热力耦合模型的结冰水库水温-冰情反演及预测方法，所述热力耦合模型包括水动力学模块、水温模块和冰模块；所述水温模块包括常规水温模型和水体过冷却模型；所述冰模块包括水内冰输移模型、静态冰盖生成判断模型和冰盖热力生消模型；对于给定的时间段，所述结冰水库水温-冰情反演及预测方法包括以下步骤：

s1依据水库边界条件以及当前时刻的水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量初始值，利用水动力学模块、常规水温模型和水内冰输移模型获取当前时刻的结冰水库流场参数、水温分布、水内冰浓度；

s2遍历水库表层，识别遍历位置是否已存在冰盖，若已存在冰盖，进入步骤s7；若不存在冰盖，则进入步骤s3；

s3遍历水库表层，判断遍历位置水温是否小于零，若遍历位置水温小于零，进入步骤s4；若遍历位置水温不小于零，进入步骤s6；

s4通过水体过冷却模型获取遍历位置的水内冰生成情况，进而对遍历位置水温进行修正；

s5根据遍历位置水内冰浓度大小，判断是否由水内冰转化为冰盖；若转化为冰盖，进入步骤s7；若没有转化为冰盖，进入步骤s8；

s6通过遍历位置流速及冰盖生成判断模型判断静态冰盖是否生成，若生成静态冰盖进入步骤s7，若没有生成静态冰盖进入步骤s8；

s7通过冰盖热力生消模型获取冰盖厚度，之后进入步骤s8；

s8在当前时刻基础上，增加时间步长；若增加后的时间不大于设定时间，进入步骤s9；若增加后的时间大于设定时间，进入步骤s10；

s9将当前时刻得到的水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量值作为下一时刻的水库流场参数、水温分布、水内冰浓度和冰盖厚度初始值，在当前时刻基础上增加时间步长后的时间作为下一时刻，返回步骤s1；

s10当前程序结束，得到水温分布情况以及冰盖分布情况。

本发明以结冰水库水域为研究对象，将其沿水流方向和垂直方向划分为若干控制单元，以每个控制单元水温、表层冰盖情况表征其所在位置的水温和冰盖情况。所述水温分布情况包括结冰水库的库区水温时空分布、下泄水温随时间变化过程；所述冰盖分布情况包括结冰类型(例如水内冰，静态冰盖)、结冰与解冻的顺序以及水库表层冰厚的时空分布。所述水库边界条件可以通过本领域常规手段来设定。给定时间段内初始时刻的结冰水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量值可以为测定值或者直接给定。

所述水动力学模块用于模拟结冰库区流场，计入了风应力对表层流速的影响。本发明中，所述水动力学模块为结冰水库宽度平均的水库立面二维水动力学模型，是将三维湍浮力流的时均连续方程、动量方程、k-ε方程沿河宽方向积分得到，其包括时均连续模型、动量模型和k-ε紊流模式模型组；

所述时均连续模型为：

所述动量模型为：

所述k-ε紊流模式模型组为：

其中：b为河宽，m；u和w分别为结冰水库水流x方向和结冰水库水深z方向上的速度分量，m/s；νe为分子粘性系数ν和紊动涡粘系数νt之和，m²/s；νt＝cμk²/ε；ρw为水体参考温度下的密度，kg/m³；p为压强，pa；β为水的热膨胀系数，℃^-1；δt为水温初始值与参考温度的差异，℃；g为重力加速度；k为紊动动能；ε紊动动能耗散率；σk和σε分别为紊动动能和耗散率的普朗特数；cμ、cε1和cε2为经验常数；gk为紊动动能生成项，gb为浮力生成项，

所述水温模块考虑了水气、冰水界面热交换以及冰水相互转化对水温和水内冰浓度的影响。其中，常规水温模块用以模拟结冰水库库区的水温场。水体过冷却模型单独运行，用以模拟水内冰的生消和水体的过冷却，其得到的水温变化值用来修正常规水温模型中的水温，得到的水内冰浓度变化值作为下一时刻水内冰输移模型的源项。该水温模块为冰模块提供了水温边界条件。

所述常规水温模型为：

当水体发生过冷(即tw,s＜0)时，所述水体过冷却模型为：

其中，tw为水温，℃；cp为水的比热容，j/kg·℃；为透过水体表层，穿过z平面的太阳短波辐射通量，w/m²；σt为温度的普朗特数；ρice为冰的密度，kg/m³；lice为冰的潜热，j/kg；cv为单位体积内的水内冰浓度；为单位体积内由于热力生消的水内冰浓度。

进一步，将式(5)展开，并省略tw²的项，水体过冷却模型变换为：

所述水内冰输移模型模拟了水内冰的输移和上浮以及与冰盖的转化。所述水内冰输移模型为：

其中，sice为单位体积内由于水内冰热力生消而产生的源项；kw为水的热传导系数，w/(m·℃)；为单个水内冰颗粒与水体的热交换nusselt数；nf＝cv/v0，为单位体积内水内冰的晶体个数；v0为水内冰颗粒的平均体积，m³；αo＝πdf·de，为正对α轴的单个水内冰颗粒表面积，m²，df为水内冰颗粒α轴的长度，m，de为水内冰颗粒的厚度，m；wb为水内冰在z方向的紊动上浮速度，m/s。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s1中，依据当前时刻的水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量初始值，然后按照以下分步骤得到当前时刻的结冰水库流场参数、水温分布、水内冰浓度：

s101通过水动力学模块获取结冰水库的流场参数；

s102通过常规水温模型获取结冰水库的水温分布；

s103通过水内冰输移模型获取结冰水库的水内冰浓度分布；

s104判断流场参数、水温和水内冰浓度是否收敛，若收敛进入步骤s2，否则返回步骤s101进行迭代计算。

本发明中，流场与水温场进行耦合求解，在迭代计算过程中利用上一次迭代常规水温模型获取的温度值作为下一次迭代步骤s101中输入值，从而对动量模型和k-ε紊流模式模型组进行修正。上述水库流场参数、水温和水内冰浓度的收敛条件是指相应变量的误差余量与结冰水库入口处相应变量值的比值小于或等于设定的容许误差，误差余量指的是通过对相应变量方程离散处理后的方程两边差值。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s2中，遍历水库表层每个控制单元体，依据遍历位置冰盖厚度是否大于0，识别遍历位置是否已存在冰盖，若大于0，表明遍历位置存在冰盖，若不大于0，则进入步骤s3。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s3中，遍历水库表层每个控制单元体，识别遍历位置水温是否小于零。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s4中，通过上述公式(5)或(8)给出的水体过冷却模型得到遍历位置的温度变化值，依据温度变化值可以获取结冰水库遍历位置的水内冰生成情况，例如若温度变化值为正值，存在放热现象，表明水内冰生长；若温度变化值为负值，存在吸热现象，表明水内冰融化。因此，将步骤s1得到的遍历位置的水温加上温度变化值即为遍历位置的实际温度，进而实现对表层遍历位置水温的修正。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s5中，根据步骤s1得到的遍历位置的水内冰浓度大小，判断是否由水内冰转化为冰盖。本发明，当遍历位置的水内冰浓度大于0.8时，水内冰转化为冰盖，否则表明水内冰没有转化为冰盖。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s6中，对遍历位置流速初步判断，若遍历位置流速小于0.07m/s，进一步通过冰盖生成判断模型来识别是否生成静态冰盖；若不满足直接进入步骤s8。所述静态冰盖生成判断模型采用matousek经验判别，考虑了初始冰盖生成时的热力学和动力学条件。冰盖热力生消模型以冰盖的一维热传导方程为基础，考虑了水气、冰水界面热交换，并在冰盖上表面还考虑了雪层对冰气热交换的影响。所述静态冰盖生成判断模型为：

当tws>0℃、tw,s≤tcr时，静态表层初始薄冰盖形成；其他情况均不考虑，直接进入步骤s8。

其中：tw,s为结冰水库表层的判断水温，℃；tws为结冰水库表层的平均水温，℃；为水库表层水气热交换通量，w/m²；uws为结冰水库表层的纵向流速，m/s；b为与结冰水库表层水面宽度有关的系数；wwind为水面上空2m处的风速；tcr为设定的静态冰盖生成的结冰水库表层经验判别临界值，℃。

水库冰盖(包括静态冰盖和水内冰转化成的冰盖)形成后，在“气-雪-冰-水”系统下的热量传递和作用下，冰盖会发生厚度变化。

上述基于热力耦合模型的结冰水库水温-冰情反演及预测方法，步骤s7中，所述冰盖热力生消模型的建立基于以下假定：(1)冰盖的融化发生在冰盖的上、下表面，冰盖上表面的融冰水可通过冰盖裂隙或蒸发等作用及时得到排除；(2)冰盖始终浮于水库表面，并随水库水位的涨落而变动，忽略冰盖对其下水流的阻力作用；(3)冰水交界面处的水温始终处于冰点温度，即0℃；(4)冰盖及雪层内部温度呈线性分布。

所述冰盖热力生消模型为：

当气温ta≤0℃时，冰盖厚度增长期，冰盖厚度的变化主要发生在冰盖底部，其冰盖厚度的变化：

当气温ta>0℃时，冰盖开始融化，融化可在冰盖表面和底部同时发生，冰盖厚度变化：

表面：

底部：

其中，hice为冰盖厚度；φr为冰盖吸收的净太阳短波辐射通量；αice、βice为考虑了长波辐射、蒸发、热传导过程的经验系数；ta为空气温度，℃；ts为冰盖表面温度，若有雪层覆盖则为雪层表面温度，℃；ks为雪层的热传导系数，w/(m·℃)；kice为冰盖的热传导系数，w/(m·℃)；hs为雪层厚度，m；tice-f为水体冰点；hw-ice为冰水热交换系数，w/(m²·℃)，hw-ice＝ρwcwchuws，ch为无量纲的冰水热交换系数。

本发明提供的基于热力耦合模型的结冰水库水温-冰情反演及预测方法具有以下有益效果：

(1)本发明通过对由水动力学模块、水温模块和冰模块构建的热力耦合模型进行处理，将水动力-水温-冰耦合作用考虑在内，能够实现高寒不同类型结冰水库的水温、冰情的热力演变过程，从而实现对结冰水库水温及冰情分布情况的反演及预测，进而为水库修建提供有效的数据支持；

(2)本发明从热力学和动力学角度建立水体过冷却、水内冰生成与输移、初始冰盖形成的冰水转化模型，能够较为全面的模拟水库冰的热力学行为过程；

(3)本发明全面系统考虑“气-雪-冰-水”的冰厚热力生消，更能真实反映冰-水的能量传递过程以及冰盖厚度的变化；

(4)本发明考虑了水体过冷却对冰水之间转化的能量交换过程的影响，从而能够准确获取结冰水库水温的演变过程。

附图说明

图1为本发明实施例热力耦合模型组成框图。

图2为本发明实施例所涉及的“气-雪-冰-水”系统热交换发生和发展过程示意图。

图3为本发明实施例基于热力耦合模型的结冰水库水温-冰情反演方法流程示意图。

图4为本发明实施例通过冰模块获取结冰水库的结冰情况分步骤流程示意图。

图5为本发明实施例对东北fm水库离散成的网格示意图，其中横坐标表示以结冰水库入库边界为起点沿水流方向的距离，纵坐标表示高程。

图6为本发明实施例对东北fm水库反演得到的结冰期2010年12月20日的库区水温二维分布图，其中横坐标表示以结冰水库入库边界为起点沿水流方向的距离，纵坐标表示高程。

图7为本发明实施例对东北fm水库反演得到的水库下泄水温与实测值的对比变化曲线图。

图8为本发明实施例对东北fm水库反演得到的坝前水库冰盖厚度随时间的变化曲线及与实测冰厚的对比。

具体实施方式

结合附图对本发明各实施例的技术方案进行清楚、完整的描述，显然，所描述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施例，都属于本发明。

实施例

本实施例运用本发明所提供的基于热力耦合模型的结冰水库水温-冰情反演及预测方法获取东北fm(丰满)水库在2010年9月2日～2011年5月31日的水温分布情况以及冰盖分布情况。

首先采用不等距矩形网格，在水库大断面、断面间距、水平网格间距分布、垂向网格间距分布基础上进行网格划分。矩形网格大小可以根据水库特征进行设定。

本实施例中fm水库被离散成532×52个矩形网格，纵向(沿主流方向)网格尺寸为20～800m，垂向网格尺寸为(沿水深方向)为0.3～3m，划分后网格如图5所示。每一个网格即为一个控制单元，以每个控制单元水温、表层冰盖情况表征其所在位置的水温和冰盖情况。

本实施例东北fm水库所涉及的边界条件和给定时间段(2010年9月2日～2011年5月31日)的初始条件如下：

(一)结冰水库边界条件

所涉及结冰水库边界条件包括水库上游入库边界条件、坝出口断面边界条件、水库表面边界条件和水库库底边界条件。

(1)水库上游入库边界条件为：

水库入库边界(即进口)水温采用库尾的实测水温或给定值；若有多条支流，将库尾下泄水温、库区支流水温按流量加权作为入库水温边界条件。

水库上游入库边界处的流速(u、w)设定为均匀流速，k和ε按照以下公式计算得到；

k＝0.00375u²(13)；

h0为进口断面水深，m。

(2)坝体和坝出口断面边界条件为：

坝体表面采用无滑移边界、且为绝热边界；坝出口断面为充分发展的湍流，有w＝0。

(3)水库表面边界条件为：

(31)水面采用“刚盖假定”，即

(32)水面无冰条件下，面表层计入水气热交换，其计算方式为：

φwa＝φsnw+φan-φbr-φe-φc(15)；

其中，为水库表层水体净吸收的太阳短波辐射，为净吸收的大气长波辐射，为水体长波的返回辐射，为水面蒸发热损失，为热传导通量；

φsnw＝βwφs(1-γw)(16)；

φan＝σεa(273+ta)⁴(17)；

式中：βw为水库表层水体对太阳辐射的吸收系数，取值0.65；γw为水体表面反射率，取值0.1；σ为stefan-boltzman常数，取值为5.67×10^-8w/(m²·k⁴)；εa为大气长波发射率，可由如下公式计算：

式中：ta为空气温度，℃；cr为云层覆盖率；ea为水面上空气蒸发压力，hpa。

φbr＝σεw(273+tws)⁴(20)；

式中：εw为水体长波发射率，取0.965；tws为结冰水库表层平均水温，℃。

φe＝f(w)(es-ea)(21)；

式中：f(w)为风函数，w/(m²·hpa)；es为空气饱和蒸发压力，hpa；ea为水面上空气蒸发压力，hpa；hum为相对湿度。

风函数f(w)反映了自由对流和强迫对流对蒸发的影响，可采用以下计算公式：

式中：w为水面上10m处的风速，m/s；δt′为水库表层温度与气温的差值，℃。

φc＝0.627f(w)(tws-ta)(25)。

(33)太阳短波辐射在雪层和冰盖内部均按指数衰减，因此冰盖吸收的净太阳短波辐射通量可表示为：

(34)当水库表层无冰时，其进入到水体沿深度方向仍以指数函数衰减：

当水库表层结冰后，冰盖及其上雪层对太阳辐射具有一定的阻挡作用，其进入到水体沿深度方向仍以指数函数衰减：

其中，为到达地面的太阳辐射通量，w/m²；rice为冰(雪)表面对太阳辐射的反射率，取0.2；hs为雪层厚度，m；ηice为太阳短波辐射在冰盖内的衰减系数，取值1.0，m^-1；ηs为太阳短波辐射在雪层内的衰减系数，取值37.6，m^-1；hice为冰盖厚度，m；ηw为太阳短波辐射在水体内的衰减系数，与水体的透明度和水质相关，取值0.5，m^-1；z为从水面起沿水深方向的深度，m。

(4)水库库底边界条件为：库底采用无滑移边界条件，且为绝热边界。

边界条件中涉及的气象条件(包括气温、云量、相对湿度、风速和太阳辐射等)由气象站提供。

所涉及时间段内的水库调度(包括入库总流量、出库流量、水库水位和电站进水口高程等)采用水电站的实测数据。

(二)初始时刻的水库初始条件

本实施例中采用2010年9月2日观测的水库库区水温分布t(i,j)作为第(i,j)个控制单元体的初始水温场，初始水温场由各测量断面的垂向水温数据线性内插到由水库流向位移(距离水库入库位置的距离)-水深(距离水库库底的距离)坐标系中得到。

并给定如下初始流场、水内冰浓度和冰盖厚度条件为：w(i,j)＝0，u(i,j)＝qij/ai，p(i,j)＝0，k(i,j)＝0，ε(i,j)＝0，cv(i,j)＝0，hice(i,j)＝0。

本实施例所采用的热力耦合模型，如图1所示，包括水动力学模块、水温模块和冰模块；上述水温模块包括常规水温模型和水体过冷却模型；上述冰模块包括水内冰输移模型、静态冰盖生成判断模型和冰盖热力生消模型。本实施例提供的基于热力耦合模型的结冰水库水温-冰情反演方法核心思想是在对结冰水库水动力、水温计算的基础上，进一步开展水温与水体过冷却、水内冰的相互影响计算，并根据水动力学及表层单元温度分布情况进行不同类型的初始冰盖判断和生成，之后按“气-雪-冰-水”系统的热量交换和平衡过程实现冰盖的热力生消过程，进而实现水库冰-水水动力与热力的耦合计算。

上述水动力学模块为沿结冰水库宽度平均的水库立面二维水动力学模型，是将三维湍浮力流的时均连续方程、动量方程、k-ε方程沿河宽方向积分得到，其包括时均连续模型、动量模型和k-ε紊流模式模型组；

上述时均连续模型为：

上述动量模型为：

所述k-ε紊流模式模型组为：

其中：b为河宽，m；u和w分别为结冰水库水流x方向和结冰水库水深z方向上的速度分量，m/s；νe为分子粘性系数ν和紊动涡粘系数νt之和，m²/s；νt＝cμk²/ε；ρw为水体参考温度20℃下的水体密度，0.99822×10³kg/m³；p为压强，pa；β为水的热膨胀系数，℃^-1；δt为水温初始值与参考温度(20℃)的差异变化，℃；g为重力加速度，9.8m/s²；k为紊动动能；ε紊动动能耗散率；σk和σε分别为紊动动能和耗散率的普朗特数，分别取1.0和1.3；cμ、cε1和cε2为经验常数，分别取0.09、1.44和1.92；gk为紊动动能生成项，gb为浮力生成项，该项在水温稳定分层时可抑制紊动动能的生成，削弱热量向下的传递，是水库保持稳定分层的重要因素，

上述水温模块考虑了水气、冰水界面热交换以及冰水相互转化对水温和水内冰浓度的影响。其中，常规水温模块用以模拟结冰水库库区的水温场。水体过冷却模型单独运行，用以模拟水内冰的生消和水体的过冷却，其得到的水温变化值用来修正常规水温模型中的水温，得到的水内冰浓度变化值作为下一时刻水内冰输移模型的源项。该水温模块为冰模块提供了水温边界条件。

上述常规水温模型为：

当水体发生过冷(即tw,s＜0)时，所述水体过冷却模型为：

其中，tw为水温，℃；cp为水的比热容，4200，j/kg·℃；为透过水体表层，穿过z平面的太阳短波辐射通量，w/m²；σt为温度的普朗特数，取0.9；ρice为冰的密度，917，kg/m³；lice为冰的潜热，334000，j/kg；cv为单位体积内的水内冰浓度；为单位体积内由于热力生消的水内冰浓度。

进一步，将式(5)展开，并省略tw²的项，水体过冷却模型变换为：

上述水内冰输移模型模拟了水内冰的输移和上浮以及与并盖的转化。上述水内冰输移模型为：

其中，sice为单位体积内由于水内冰热力生消而产生的源项；kw为水的热传导系数，0.58，w/(m·℃)；为单个水内冰颗粒与水体的热交换nusselt数，取值4.0；nf＝cv/v0，为单位体积内水内冰的晶体个数；v0为水内冰颗粒的平均体积，m³；αo＝πdf·de，为正对α轴的单个水内冰颗粒表面积，m²，df为水内冰颗粒α轴的长度，2.0×10^-3，m，de为水内冰颗粒的厚度，1.0×10^-5，m；wb为水内冰在z方向的紊动上浮速度，0.001，m/s。

上述静态冰盖生成判断模型采用matousek经验判别，考虑了初始冰盖生成时的热力学和动力学条件。冰盖热力生消模型以冰盖的一维热传导方程为基础，考虑了水气、冰水界面热交换，并在冰盖上表面还考虑了雪层对冰气热交换的影响。

上述静态冰盖生成判断模型为：

当tws>0℃、tw,s≤tcr时，静态表层初始薄冰盖形成；其他情况均不考虑。

其中：tw,s为结冰水库表层计算水温，℃；tws为结冰水库表层平均水温，℃，由步骤s1计算得到的水库表层各控制单元体水温平均得到；为水气热交换通量，w/m²；uws为水库表层的纵向流速，m/s；b为与结冰水库表层水面宽度有关的系数，若b≤15.0m时，取b＝15.0；若b＞15.0m时，取b＝-0.9+5.7lnb；wwind为水面上空2m处的风速；tcr为设定的静态冰盖生成的结冰水库表层经验判别温度，-2.0，℃。

值得注意的是，式(9)中初始冰盖的形成需要满足水库表层单元的纵向流速uws<0.07m/s条件；对于初始时刻，按该方式确定的初始冰盖厚度，建议设定为0.001m。当流速大于0.07m/s，表层水体的紊流强度增加，水体发生过冷却现象并产生水内冰，当水内冰浓度大于0.8时，水库初始冰盖可由水内冰转化得到，水内冰转化成的冰盖初始厚度为表层网格厚度，但表层网格厚度不宜超过0.3m。

水库初始冰盖形成后，在“气-雪-冰-水”系统下的热量传递和作用下，冰盖会发生厚度变化。如图2所示，本实施例中“气-雪-冰-水”系统，是一个多介质、具有复杂的热交换发生和发展过程的系统。对于结冰期，寒冷大气通过雪和冰的阻冷效应后，进一步作用于冰水界面，在与冰-水热交换进行热平衡后实现冰的增长过程；此外，太阳辐射通过雪层、冰层的吸收和穿透后，仍可少部分进入水体并向水体提供热量。对于消融期，冰盖上、下表面分别受温暖大气和太阳辐射、冰水热交换而发生融化，此时冰盖由于上下界面同时处于零温状态而内部无热交换。本实施例的冰盖厚度热力生消模型正是根据上述机理建立起来的。

所述冰盖热力生消模型的建立基于以下假定：(1)冰盖的融化发生在冰盖的上、下表面，冰盖上表面的融冰水可通过冰盖裂隙或蒸发等作用及时得到排除；(2)冰盖始终浮于水库表面，并随水库水位的涨落而变动，忽略冰盖对其下水流的阻力作用；(3)冰水交界面处的水温始终处于冰点温度，即0℃；(4)冰盖及雪层内部温度呈线性分布。

所述冰盖热力生消模型为：

当气温ta≤0℃时，冰盖厚度增长期，冰盖厚度的变化主要发生在冰盖底部，其冰盖厚度的变化：

当气温ta>0℃时，冰盖开始融化，融化可在冰盖表面和底部同时发生，冰盖厚度变化：

表面：

底部：

其中，hice为冰盖厚度；为冰盖吸收的净太阳短波辐射通量；αice、βice为考虑了长波辐射、蒸发、热传导过程的经验系数，分别取100、15；ta为空气温度，℃；ts为冰盖表面温度，若有雪层覆盖则为雪层表面温度，℃；ks为雪层的热传导系数，0.3，w/(m·℃)；kice为冰盖的热传导系数，2.24，w/(m·℃)；hs为实时监测的雪层厚度，m；tice-f为水体冰点，0℃；hw-ice为冰水热交换系数，w/(m²·℃)，hw-ice＝ρwcwchuws，ch为无量纲的冰水热交换系数，取值为1.1×10^-3。

本实施例采用有限体积法对各模型进行离散，得到各模型的离散方程。

基于上述情形，如图3-4所示，本实施例按照以下步骤对东北fm水库在2010年9月2日～2011年5月31日的水温分布情况以及冰盖分布情况进行反演：

s1依据水库边界条件以及当前时刻的水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量值，利用水动力学模块、常规水温模型和水内冰输移模型获取当前时刻的结冰水库流场参数、水温分布、水内冰浓度。

依据当前时刻的水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量初始值(对于初始时刻，当前时刻水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量值对应的是初始条件)，然后按照以下分步骤得到当前时刻的结冰水库流场参数、水温分布、水内冰浓度：

s101通过水动力学模块获取结冰水库的流场参数。

本步骤，依据边界条件和当前时刻结冰水库流场各变量初始值、水温分布相应变量初始值，通过simple算法对均连续模型(式(1))、动量模型(式(2))和k-ε紊流模式模型组(式(3))的离散方程进行求解，得到流场参数，包括u(i,j)、w(i,j)、p(i,j)、k(i,j)、ε(i,j)。

s102通过常规水温模型获取结冰水库的水温分布。

本步骤，依据边界条件、当前时刻的水温初始值以及步骤s1得到的流场参数对常规水温模型(式(4))的离散方程进行求解，得到水温场参数，tw(i,j)。

s103通过水内冰输移模型获取结冰水库的水内冰浓度分布。

本步骤，依据边界条件、当前时刻水内冰浓度初始值以及步骤s101得到的流场参数和步骤s102得到的水温场参数对水内冰输移模型(式(6))的离散方程进行求解，得到当前时刻的水内冰浓度参数，cv(i,j)。

s104判断流场参数、水温和水内冰浓度是否收敛，若收敛进入步骤s2，否则返回步骤s101进行迭代计算。

本步骤，水库流场参数、水温和水内冰浓度的收敛条件是指相应变量的误差余量与结冰水库入口处相应变量值的比值小于或等于设定的容许误差(这里取2％)，误差余量指的是通过对相应变量方程离散处理后的方程两边差值。

为了提升上述获取水库流场参数、水温分布、水内冰浓度过程的准确度和反演效率，本实施例进一步设置循环变量n作为收敛参数，当水温、水内冰浓度和冰盖厚度误差余量不满足设计精度时，进入下一个迭代过程，同时循环变量n增加1，直至水温、水内冰浓度和冰盖厚度误差余量满足设计精度。在迭代计算过程中，流场与水温场进行耦合求解，利用上一次迭代常规水温模型获取的温度值作为下一次迭代步骤s101中输入值，从而对动量模型和k-ε紊流模式模型组进行修正。

s2遍历水库表层，识别遍历位置是否已存在冰盖，若已存在冰盖，进入步骤s7；若不存在冰盖，则进入步骤s3；

本步骤，遍历水库表层每个控制单元体，识别遍历位置(这里指的是每个控制单元体，以下同)是否已存在冰盖，若hice＞0m，表明遍历位置存在冰盖，进入步骤s7；若hice＝0m，表明遍历位置不存在冰盖，进入步骤s3。

s3遍历水库表层，判断遍历位置水温是否小于零，若遍历位置水温小于零，进入步骤s4；若遍历位置水温不小于零，进入步骤s6。

本步骤，遍历水库表层每个控制单元体，判断其水温是否小于零，当存在水温tw小于0的地方(根据每个控制单元体获取的水温进行判断)，表明水内冰生成，需要启动水体过冷却模型，获取水内冰生成情况，进而对水温进行修正。若不存在水温tw小于0的地方，表明不存在水体过冷却现象，需要进一步对冰盖情况进行进一步分析。

s4通过水体过冷却模型获取遍历位置的水内冰生成情况，进而对遍历位置水温进行修正。

本步骤，通过上述公式(5)或(8)给出的水体过冷却模型得到遍历位置的温度变化值，依据温度变化值可以获取水库遍历位置的水内冰生成情况，例如若温度变化值为正值，存在放热现象，表明水内冰生长；若温度变化值为负值，存在吸热热现象，表明水内冰融化。因此，将步骤s1得到的遍历位置的水温加上温度变化值即为遍历位置的实际温度，进而实现对表层遍历位置水温的修正。本步骤完成后，进入步骤s5。

s5根据遍历位置水内冰浓度大小，判断是否由水内冰转化为冰盖；若转化为冰盖，进入步骤s7；若没有转化为冰盖，进入步骤s8。

本步骤，继续根据步骤s1得到遍历位置的水内冰浓度大小，判断是否由水内冰转化为冰盖。本发明，当遍历位置的水内冰浓度大于0.8时，水内冰转化为冰盖，否则表明水内冰没有转化为冰盖。

不管任何时刻，只要存在水内冰转化为冰盖，均给定遍历位置的冰盖厚度为表层网格厚度(本实施例中为0.3m)，然后进入步骤s7。

s6通过遍历位置流速及冰盖生成判断模型判断静态冰盖是否生成，若生成静态冰盖进入步骤s7，若没有生成静态冰盖进入步骤s8。

静态冰盖的形成需要满足水库表层单元的纵向流速uws<0.07m/s，若满足，进一步通过静态冰盖生成判断模型(式(9))来判断是否有静态冰盖生成；若不满足直接进入步骤s8。

需要注意的是，对于初始形成静态冰盖的时刻，给定遍历位置的冰盖厚度为0.001m，然后进入步骤s7。

s7通过冰盖热力生消模型获取冰盖厚度，之后进入步骤s8。

不管是步骤s2、还是步骤s5、或步骤s6，都是为了计算冰盖的厚度。依据步骤s2遍历水库，得到的冰盖厚度是水库表层当前时刻之前冰盖分布情况。由于步骤s5是对水内冰情况进行分析，因此按照步骤s5给出的遍历位置得到的冰盖厚度是水内冰转化为冰盖的厚度。由于步骤s6是对无冰盖区域的静态冰盖生成情况进行分析，因此按照步骤s6给出的遍历位置得到的冰盖厚度是初始静态冰盖的厚度。

本步骤，依据边界条件、上一时刻的冰盖厚度以及步骤s1得到的水温场参数对冰盖热力生消模型(式(10)-(12))的离散方程进行求解，得到当前时刻的冰盖厚度参数，hice(i,j)。

s8在当前时刻基础上，增加时间步长；若增加后的时间不大于设定时间，进入步骤s9；若增加后的时间大于设定时间，进入步骤s10。

s9将当前时刻得到的水库流场参数、水温分布、水内冰浓度和冰盖厚度各变量值作为下一时刻的结冰水库流场参数、水温分布、水内冰浓度和冰盖厚度变量初始值，在当前时刻基础上增加时间步长后的时间作为下一时刻，返回步骤s1。

s10当前程序结束，得到水温分布情况以及冰盖分布情况。

在当前时刻基础上增加一个时间步长δt，并判断增加后的时间是否达到设定的时间段上限tend，若没有达到设定的时间段上限tend；将利用当前时刻得到的结冰水库流场参数、水温分布、水内冰浓度和冰盖厚度对各变量进行赋值，作为下一时刻的水库流场参数、水温分布、水内冰浓度和冰盖厚度变量初始值；增加后的时间作为下一时刻，然后返回步骤s1，重复上述步骤s1-s8的操作；若达到设定的时间段上限，则程序结束。对计算得到的水温数据、水内冰数据和冰盖厚度数据进行汇总统计，便可得到反演的结冰水库的水温分布情况以及冰盖分布情况。所述水温分布情况包括结冰水库的库区水温时空分布、下泄水温随时间变化过程；所述冰盖分布情况包括结冰类型(例如水内冰，静态冰盖)、结冰与解冻的顺序以及水库表层冰厚的时空分布。

对于任一给定的时刻，将所得水温数据添加到由水库流向位移(距离水库入库位置的距离)-高程坐标系中便可得到库区水温二维分布图。图6给出了对东北fm水库反演得到的结冰期2010年12月20日的库区水温二维分布图。该时段fm水库入库温度约为1.5℃，但在沿程-16℃气温作用下，水体沿程失热，并不断降温，而库区较深水域由于蓄热作用，仍具有6℃的高温区域，沿程温差达到4.5℃。此外，还可观察到，库尾的低温水行进到水深较大、水温较高的区域，与该段水体发生纵向和垂向混合。当混合区水温到达4℃后，库尾低温水的进一步行进和负气温的作用促成了水深较大区域的逆温分层现象，进而产生表层冷水浮力流动，这种表层冷水浮力流动逐渐向坝前水域推进。已形成逆温部分的表层水体，在负气温作用下，其水温进一步失热，逐渐降低至冰点附近，满足成冰条件后产生冰情。

对水库某一位置在给定时间段内的水温数据进行统计，便可得到该位置水温随时间的变化趋势。图7给出了对东北fm水库反演得到的水库坝下0km电站尾水下泄水温与坝下0.2km一期尾水下泄水温实测值的对比变化曲线图。整体来看，本发明提供的模型能够较好地模拟出fm水库冰期及前后下泄水温的变化趋势，即降温期下泄水温逐渐下降，冰期下泄水温基本保持平稳，开河后升温期下泄水温逐渐升高。降温期(9月～11月)模拟水温与实测值的误差在0.4℃以内；水库冰期下泄水温的模拟精度略差，表现在结冰期和稳封期模拟水温偏高，开河期模拟水温偏低，可能是由于模型未考虑风浪对结冰期冰的破坏以及库汊结冰(融冰)对主库水温的影响。

对水库某一位置在给定时间段内的冰盖厚度数据进行统计，变可得到该位置水温随时间的变化趋势。图8给出了对东北fm水库反演得到的坝前水库冰盖厚度随时间的变化曲线及与实测冰厚的对比。2010～2011年冰期，fm水库坝前封河日期为2011年1月2日，开河日期为2011年4月17日，整个冰期历时106天；反演得到的坝前封河日期为2010年12月31日，较实际日期提早2天，而反演得到的坝前开河日期为2011年4月17日，与实际日期一致。反演得到的冰厚均方根误差为4.58cm，最大相对误差为10.9％(2011年2月7日)；实测最大冰厚为70.0cm(2011年3月11日)，而反演值为70.8cm(2011年3月15日)，最大冰厚的相对误差为1.0％；特别是封冻初期(1月)和稳封末期(3月)模拟值与实测值基本一致。2月上旬至中旬冰厚模拟的误差较大，相对误差在4.9～10.9％。上述结果表明，本发明建立的模型可以很好地模拟结冰水库的冰厚演变过程。

本实施例基于东北fm水库边界条件及2010年9月2日初始条件，利用本发明提供的方法实现了在2010年9月2日～2011年5月31日的水温分布情况以及冰盖分布情况的反演。此外，也可以依据结冰水库边界条件及给定时间段内的初始时刻水温及冰情分布(即初始条件)，利用本发明提供的方法对给定时间段内的水库水温分布及冰盖分布情况进行预测。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。