少资料条件河渠冬季冰情过程模拟研究

王 涛，郭新蕾，陈玉壮，刘吉峰，佘云童，潘佳佳

（1.中国水利水电科学研究院 流域水循环与调控国家重点实验室 北京 100038；2.阿尔伯塔大学，阿尔伯塔省 加拿大； 3.黄河水利委员会水文局，河南 郑州 450004）

1 研究背景

北半球较高纬度的地区约有60%以上河流在冬季会经历冰凌过程，我国北纬30°以北占四分之三的国土面积上都有冰凌现象发生。河冰发展演变包括流凌、冰盖形成、冰盖热增长和封河、冰盖的热衰减和消融开河等过程，是复杂的水文、气象、水力学、热力学和动力学交互作用的结果。为了准确模拟和预报河渠冰凌发展全过程，国内外学者在该领域做了大量工作。较为系统的模型可从Shen的RICE模型算起，提出的模拟河冰过程的双层解析框架，考虑了水温分布和冰的浓度分布以及冰盖热力增长和消退。Beltaos［1］研制了RIVJAM模型计算宽河型冰塞所引起的水位升高，该模型可以较好的模拟非平衡的冰塞或接地冰塞。类似的模拟模型还包括ICEJAM［2］，HEC-RAS，ICEPRO，ICESIM［3］，River 1D［4］和RIVICE［5］等。Shen等［6-8］建立了具有世界水准的一维、二维DynRICE模型（CRISSP）来模拟冰凌发展过程，并考虑了河床变化和泥沙运动，该系列模型已成功应用于世界多条河流的冰情研究中。最近Wazney等［9］在上述模拟冰盖形成演化时提出了新的公式以尝试建立热力学和动力学之间的联系并得到验证。国内杨开林［10］、王军［11］、茅泽育［12］、穆祥鹏［13］、李润玲［14］等也建立了冰塞形成及演变发展动态数学模型、河道冰塞堆积厚度的数值模型及其冰厚增长度日法等，计算了典型河道中的冰塞堆积厚度，郭新蕾等［15］针对长距离明渠-闸门-泵站系统冬季正、反向输水可能出现的冰问题，开发了调水工程冬季输水冰情过程模拟平台并在南水北调、密云水库调蓄工程中取得较好模拟效果。此外，Morse［16］、郭新蕾［17］、王涛［18］、Sun等［19-20］先后将统计模型、人工神经网络、模糊数学方法等与冰情预测预报结合，较为成功的预测了封开河时间、冰厚、冰坝发生日期等冰情参数，也应用到黄河、黑龙江等河流，取得了较好效果。

上述进展表明冬季河冰过程的模拟研究已取得较大进展，但上述数学模型在一些天然河道冰情模拟应用时常常会遇到一个难题，即缺乏细致系统的河道边界条件、水文、气象和水力学基础资料。具体表现一是国内水文站主要设立在较大流域上，而且站与站之间断面间距长，用于计算的水文站数据空间步长大；二是测量的水文和气象资料通常为日均值或者每日最大、最小值，在时间分布上现有数据的时间步长太大；三是天然江河支流纵横，但支流通常未设水文观测站，支流流入和流出的水量未知；四是水文站观测的水温数据通常精度不够或数据不完整；五是太阳辐射和气温等相关气象资料采用气象局发布的信息，气象站同水文观测站的距离较远，数据同步性差。比如我国黑龙江上游900 km也只有8个基本水位站，主要观测水位信息［21］。黄河是国内水文站系统较为健全、水文观测数据相对系统的河流，但黄河内蒙古河段823 km也只有4个水文站［22］，2014年以后也只增加了包头水文站。1990年代，黄河水利委员会先后同芬兰和美国合作开发了黄河下游冰情预报数学模型，但黄河上的观测资料不能满足模型计算中对实测资料精度和种类的需求，再加之黄河局部河段为游荡性河道，河床地形变化频繁，导致模型建成后很快也不能适应实际运行应用的需要［23-25］。目前，北方天然河道冰情模拟中普遍存在着实测资料不全或缺失、河道断面资料难以测量等条件制约。

鉴于此，本研究提出少资料条件河渠冬季冰情发展过程模拟的实用方法和处理技术，重点研究了支流流量未知时的流量分配方案、每时气温和太阳辐射过程的计算方法等，并将上述模型和方法应用到黄河内蒙河段冰情模拟计算中，目的是为这种少资料条件河渠的冬季冰情过程模拟提供方法支撑。

2 少资料条件下支流流量分配

天然江河进行水动力学计算时如果沿途支流流量信息缺乏，有必要根据实测流量和水位进行模拟反演，通过辨识沿途支流流量分配比例以达到计算中区间流量的平衡。

明渠非恒定流的连续性方程和运动方程［26］为

式中：x为距离，m；t为时间，s；z为水位，m；Q为流量，m3/s；A为过流面积；B为水面宽；ql为单位流程上的侧向出流量，负值表示流入。

在进行数值计算时，需要给出研究河段的初始条件、上下游边界、支流的流量边界条件。边界条件通常有3种类型，分别为水位边界条件Z=Z（t），流量边界条件Q=Q（t），水位流量边界条件Q=Q（Z）。计算中通常采用上游边界和支流边界给定初始流量、下游边界给定初始水位的方法，前后计算断面应满足流量平衡和能量守恒。

计算河段支流边界示意图如图1所示，其中边界条件包括6个：上游边界①，下游边界②，支流边界③-⑥，在冰情模拟时常会遇到只有上游边界和下游边界的流量信息。因支流③-⑥流量未监测，计算中为了满足流量平衡，需要对区间汇入和流出的流量进行合理分配。如果仅仅简单采用上、下游边界流量差进行比例分配的话，忽略了实际流量自上游向下游传递的时间差很可能导致分配方案的不合理。本研究提出下述流量分配方法进行支流流量分配，可确保计算河段流量的平衡。

图1 计算河系支流边界示意图

在计算中只输入上游边界条件①的流量过程和下游边界②的水位条件，上游边界①流量传播到各支流边界③-④-⑤-⑥-下游边界②的模拟值为Qi（i=1，2，3，4），以黄河内蒙古河段巴彦高勒至三湖河口4个分退水口作为支流边界计算的流量变化曲线如图2所示，其中上游边界①为巴彦高勒站，下游边界为②三湖河口水文站，图2中支流边界③-⑥模拟值和下游边界2的模拟值分别为支流未分配流量条件下计算结果。图2中可以看出，下游边界②模拟流量和实测流量存在明显差异，因为模拟中区间支流流量未计入，导致流量存在不平衡。为此，需给支流边界③-⑥分配适当流量。由于支流边界③-⑥并无流量观测数据，采用实测下游边界②的流量Qdown分别减去流量传递到③-⑥支流断面处的模拟流量，差值按照比例ai（i=1，2，3，4）分配到各支流边界③-⑥作为计算中支流分配的流量qi（i=1，2，3，4），确保计算河段流量平衡，各支流分配的流量qi为：

流量分配系数αi的确定可通过调整各支流边界αi值以计算下游边界②流量模拟值，再与实测值对比，以使关键时间点上的流量值最接近率定得到，或可以尝试利用最小二乘法寻优达到相对准确辨识区间支流流量分配比例的目的，确保计算河段流量的平衡。通过支流边界③-⑥流量的合理分配，计算得到支流分配流量后下游边界②的模拟值如图2所示，由图2可知，下游边界②三湖河口水文站模拟值和实测值吻合良好。

图2 未分配支流流量传递过程模拟曲线（实测值暂按±5%的误差考虑）

3 少资料条件下每时太阳辐射的计算

太阳辐射资料来自气象局，是以日为单位的净太阳辐射，而模型计算需要每小时的太阳辐射值，本研究采用计算和参数率定相结合解决缺少每时太阳辐射资料的问题。

太阳辐射计算，参考文献［27］无云状态下的短波太阳辐射φcl，可以由下式计算：

式中：φcl为短波净太阳辐射，W/m2；φso为每单位面积总的外来太阳辐射，J/m2；Iso为太阳辐射常数，冬季约为1380 W/m2；ω是时角，午时为0°，每小时变化15°，早上为正值，下午为负值；δ为太阳倾角，弧度；φ为纬度，度，北纬为正，南纬为负；dn为一年中的天数，从1月1日算起，如1月20日=20 d；m为当地气压 pa时海拔z米的光学气团，m； p0为水平面气压，Pa；m0为水平面处的光学气团；太阳纬度α=90-θz；θz=arccos( )sinδsinφ+cosδcosφcosωi，弧度； ωi为前一个小时与当前小时时角的均值；E0为地球轨道的偏心校正系数。

φcl为无云状态下的太阳辐射，在有云条件下，太阳辐射将减小，用下式估算［28］：

式中：φri为有云状态下的太阳辐射，W/m2；C为云的覆盖度，分为10个等级，C=0表示晴空无云，C=10表示乌云全部覆盖。如果云的状态没有数据，通常采用估计数据。

当太阳辐射到达水面或冰面，部分辐射将反射回大气，太阳辐射公式再次修正为φR：

式中Rt为反射率，同纬度和照射物体表面性质有关。

通常在冰情演变过程计算中，需将河段分为尽可能多的区域，确保计算断面使用的气象数据的准确度，每个区域计算一个站点的太阳辐射和气温。通过式（4）—式（10）可计算出晴空无云状态下冰情太阳辐射，式（11）计算实际有云时状态太阳辐射，但是云的覆盖程度没有资料记录，本研究采用每天日照数同统计时段2012年11月—2013年3月内最大日照数比值衡量云层覆盖状态。于是，式（11）可修正为：

式中：SSH为计算日日照数；SSHmax为统计时段内最大日照数。

冰情演变过程中太阳对冰面和水面的反射率不仅跟冰的表面状况有关，还和河流中水的含沙量、水的浑浊度等相关。实际计算中不同的河流不同冰状态反射率都不相同，一维模型计算中考虑太阳辐射时通常忽略太阳反射率的影响，导致热扩散计算中存在误差，本研究针对不同时段、不同河道冰情状况，采用不同的系数反应太阳反射率的变化，通过模型率定反演得到表达太阳反射率的系数。在冰情模拟中该系数分为流凌前Rb1、流凌期Rb2、封河期Rb3、开河前期Rb4、开河期Rb5、冰盖表面有积雪覆盖Rb66种反射率，不同时段冰水情状态如表1所示，则式（12）可写为：

表1 不同冰情状态表达太阳反射率的系数

式中Rbi为表达太阳辐射变化的系数，i=1，2，3，4，5，6。通过冰情实际情况率定得到。Rbi同冰面或水面太阳反射率Rti、河道情况Rii、河道水质Wai、冰情情况Ici和积雪覆盖情况Sni等有关。

黄河内蒙河段有磴口、临河、五原、乌拉特前旗、包头、土右旗气象局、托克托县7个典型气象站（如图3），黄河内蒙古河段冰情计算中，分别计算上述7个气象站太阳辐射值，7个气象站太阳辐射值分布规律和走势一致，在这里只列举出临河的太阳辐射计算结果，临河气象站位于东经40.73°，北纬107.37°，高程1041.1 m，图4为2012年11月1日—2013年3月31日临河每时无云状态下净太阳辐射φcl（W/m2）计算值，图5列出了2012年11月临河每时无云状态下净太阳辐射。为了更清楚了解太阳辐射局部变化规律，根据计算太阳辐射值对应的年份和日期，可率定出不同结冰时段和冰盖状态下太阳反射率的影响，并根据日照数估算出云层遮挡度的影响，代入式（14）即可求出符合实际情况的太阳辐射值。

图3 气象站分布和气象数据分区示意图

图4 2012年11月—2013年3月临河每时太阳辐射计算值

图5 2012年11月临河每时太阳辐射计算值

4 少资料条件下每时气温的计算

水文站和气象站记录的气温值为日最高值、日最低值或者日均值，且站与站之间距离远，冰情模拟需要每小时的气温值，现有资料不能满足模拟模型的需求。本研究中利用Parton and Logan（1981）［29-30］的气温计算公式进行推演。

式中：Ti为白天或者夜晚的第i个小时的温度；Y为白天时长，h；Z为夜晚时长，h；Tmax和Tmin为日最高和最低气温；TS为日落时的气温；m为最低温度出现后到日落的小时数；n为日落后到最低温度时间的小时数；a为最高气温的滞后系数；b为夜间气温系数；c为从日出时开始的最低温度的滞后系数；HR为计算温度对应的小时数，取值为1～24 h。系数a、b、c是由土壤和气温决定的，根据当地实测气温进行率定得到，本研究率定得到的a=0.15，b=2.28，c=-0.2。

仍以图3所示气象站资料来说明，根据气象站提供的日最大气温和最小气温值率定所需站点的时均气温，图6为2012年11月1日—2013年3月31日临河每时气温计算值，图7将2012年11月每时气温单独列出，便于更清楚观察气温的局部变化规律。

图6 2012年11月—2013年3月临河每时气温计算值

图7 2012年11月临河每时气温计算值

5 黄河内蒙古河段冰情模拟的应用

将上述少资料条件河渠冬季冰情发展过程模拟的实用方法和处理技术应用到黄河内蒙古河段的冰情模拟中。

该段最新河道断面数据为2012年10月的大断面资料，从巴彦高勒到头道拐共测量断面166个，因此本研究冰情计算、验证、校核及冰情模拟以2012—2013凌汛期作为研究对象。冰情演变过程模拟中必需的物理量气温、太阳辐射和支流流量分配数据均不能满足计算需要，需要采用上述方法对所需资料进行推演和计算。

5.1 少资料条件支流流量分配 模拟河段汇入支流众多，沿途分布多处用于工业、农业、畜牧业和居民饮水的取水和退水口，支流汇入没有水文观测资料，沿途取水资料模糊，因此需要通过实测流量和水位模拟，对沿途流量进行分配。内蒙古河段只有巴彦高勒、三湖河口和头道拐3个水文站流量资料。为了确保巴彦高勒和三湖河口、三湖河口和头道拐之间流量平衡，需按“少资料条件支流流量分配方法”将流量进行分配。计算中采用的支流边界如表2所示，巴彦高勒-三湖河口之间分配支流边界4个，分别为支流边界1、支流边界2、支流边界3和支流边界4，三湖河口和头道拐之间分配支流边界6个，分别为支流边界5、支流边界6、支流边界7、支流边界8、支流边界9、支流边界10和头道拐，根据研究中区间流量分配方法，经多次反算率定，建议巴彦高勒至三湖河口区间支流流量分配系数比例αi分别为0.1∶0.1∶0.3∶0.5，三湖河口至头道拐区间支流流量分配系数比例αi分别为0.1∶0.1∶0.1∶0.1∶0.3∶0.3。

5.2 少资料条件太阳辐射和气温的模拟 由于研究河段较长，流经区域气候条件有明显差异，因此将计算流域划分为8个气象数据计算区（如图3所示），区域1—区域8所采用的气温和太阳辐射资料分别对应气象站为：磴口、临河、五原、乌拉特前旗、包头、包头、土右旗气象局、托克托县。根据有限的气象资料，采用本研究少资料条件每时太阳辐射和气温的计算方法，可计算出上述8个水文站每时气温和每时太阳辐射作为8各区域计算输入值。其中临河气温和净太阳辐射计算如图4—图7所示。

5.3 黄河内蒙河段冰情发展的模拟 冰情过程模拟采用河冰过程冰水基本数学模型，包括：明渠非恒定流模型、水流的热扩散模型、冰花扩散模型、冰盖下的水流输冰能力模型、水面浮冰的输运模型、冰盖和冰块厚度发展模型、冰塞下冰花含量和冰塞厚度模型等。通过相关参数率定选择出适合黄河冰情演变规律的较为优选的方案，计算中取水和空气热交换系数hwa=16（W/m2/°C），平铺上溯模式发展的最大费劳德数Fr_jux=0.06和水力加厚模式发展的最大费劳德数，该值分段给出，即：1～50 km河段Fr_max=0.11，50～174 km河段Fr_max=0.13，174 km至下游段河段Fr_max=0.097。

巴彦高勒封河时间为2012年12月23日。图8为巴彦高勒固定冰盖厚度发展过程模拟值和实测值的比较，表2为冰盖厚度模拟值和实测值的误差。2012—2013年度冰情信息显示，截止到2012年12月27日研究河段全线封冻，根据冰厚测量数据显示，到2013年2月1日，冰厚已充分发展达到最大值。巴彦高勒整个冬季冰盖模拟值和实测值均方根误差为0.065 m，冰盖增长期（封河—2013年2月1日）均方根误差为0.056 m；巴彦高勒整个冬季实测量值和模拟值绝对误差均值是0.059 m，冰盖增长期绝对误差均值分别为0.041 m。模拟结果显示冰盖增长期模拟误差小于冰盖消融期，原因是：在冰厚增长期，热力作用为冰盖厚度增加的主要因素，而在冰盖消融期，冰盖并不是按照静止不动通过热力作用就地融化消失，通常冰盖消融期是由热力和动力共同作用，但在一维冰情模拟中，冰盖发展的数值模拟尚未能考虑到动力作用，所以冰盖模拟中消融期计算误差比增长期大。

图8 巴彦高勒冰厚模拟值和实测值比较

表2 冰盖厚度模拟值和实测值的误差

内蒙古河段最先封冻时间为2012年11月30日，位于昭君坟上游附近（109°54′12″E和40°31′53″N），该位置距离上游巴彦高勒水文站286.7 km，从该位置冰盖向上游发展过程如图9所示，2012年12月23日封冻上首发展到上游巴彦高勒。从图9上看出冬季冰盖前沿发展过程模拟值和实测值吻合较好，整个冬季模拟值和实测值的均方根误差为18.76 km，绝对误差均值6.43 km；2013年3月5日前模拟值和实测值的均方根误差仅为6.27 km，绝对误差均值2.50 km，模拟的冰盖前沿发展过程与观测值基本对应，但仍然存在开河期模拟值和实测值误差大于冰盖增长期和稳定期的情况。由上面的对比来看，提出的时均气温和时均太阳辐射计算方法、支流流量的动态分配方法能够较好的解决冰情计算中相关少资料的难题。

图9 冰盖前沿的发展过程

6 结论

针对计算河段区间退水和取水资料缺少问题，本文提出了流量动态分配方法，并应用到黄河巴彦高勒到头道拐区间10个分水口支流流量分配中，通过合理的比例分配支流流量，确保了计算区间流量的平衡。为了得到计算中时均太阳辐射资料，在计算净太阳辐射基础上，提出了云层遮挡计算方法和不同结冰期太阳反射率影响的推演方法，通过计算水温和实测水温的比较，率定出冬季不同冰情时段时均太阳辐射值。时均气温采用Parton和Logan的气温计算公式通过实测日最低、最高气温率定得到。将上述方法得到的支流流量、气象站的时均气温和太阳辐射应用到黄河内蒙古河段冰情模拟中，通过对比巴彦高勒冰厚和冰盖前沿发展过程模拟值和实测值的均方根误差和绝对误差均值对比，表明总体上模拟值和实测值吻合较好，证实所提出的时均气温和时均太阳辐射计算方法、支流流量的动态分配方法能够解决冰情计算中相关资料缺少的难题。

本文提出的不同河段支流流量分配比例的优化方法尚待进一步研究，不同冰盖结构和封河期太阳辐射率的变化将通过原型观测和理论研究相结合进一步验证和率定，并探索其变化规律。随着水文和气象观测仪器的改进，观测数据更加详实，支流流量数据、时均气温和太阳辐射数据的逐步完善，冰情模拟的精度会进一步得到提高。