黄河内蒙古河段典型凌情过程模拟

王 涛，余弘婧，郭新蕾，刘吉峰，陈玉壮，佘云童

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

1 研究背景

黄河是中国凌汛出现最为频繁的河流之一，凌汛问题一直受到中央和水利部的高度重视，随着刘家峡、龙羊峡、万家寨和小浪底水库的运行，黄河干流凌汛灾害得到很大缓解，万家寨水库以下的中下游已经连续多年未出现一定规模的凌汛灾害，主要凌汛灾害集中在最北端的内蒙古河段。黄河内蒙古河段位居倒“U”字型河道大拐弯最底部（如图1）［1］，全长823 km，冬季最低气温可达-35 ℃，特殊的地理位置导致内蒙古河段封河自下游向上游、开河自上游向下游发展。无论从地理位置、河道走势还是水文气象条件看内蒙古河段都是黄河冰情最为严重的区域。最近十几年内蒙古河段凌汛又出现了新情况［2］：冬季平均气温偏高，气温变幅大；河道淤积严重；槽蓄水增量加大。这些新特点的出现使得内蒙古河段的冰情规律发生了明显变化，现有的冰情预报和模拟模型参数不能客观地反映冰情的新特点，需要改进现有冰情模拟系统以反映当前冰情特性，更准确地把握其变化规律。

从1950年代开始就涌现一批冰情研究学者开展黄河冰情观测［3-4］、局部河段冰情模拟［5-6］、冰情特性研究［7-9］等，为全面掌握黄河冰情机理奠定了坚实基础。随着科技手段的进步、科研水平的提高及新型观测仪器的研发，黄河冰情研究得到进一步发展。李志军［10］、郜国明［11］和邓宇等［12］通过在黄河典型河段取冰块样本进行实验室分析，研究了黄河不同结冰期冰盖结构和力学性质的变化特征。王军等［13-14］改进了冰塞下冰块起动的临界条件计算形式，并借鉴泥沙输运理论推导得到冰塞面变形方程，计算黄河河曲段稳封期河道中的冰塞堆积厚度。冀鸿兰等［15］根据实测资料，总结了万家寨水库运行后上游河段冰情的变化特点。王涛和郭新蕾等［16-18］采用神经网络和模糊推理系统预报黄河流凌、封河和开河日期。上述研究从不同方面推动着黄河冰情研究的进步。

较为系统的河冰发展演变过程模拟的数学模型可从Shen 的RICE 模型算起，该模型提出的模拟河冰过程的双层解析框架，考虑了水温分布和冰花的浓度分布以及冰盖热力增长和消退［19］，之后不断改进发展为完善的一维RICEN 模型［20］。类似的模拟模型还包括ICEJAM（Flato 和Gerard）［21］、HEC-RAS、ICEPRO、ICESIM（Carson 等）［22］、River 1D（She 等）［23］和RIVICE（Lindenschmidt）［24］等。Shen 等建立了具有世界领先水平的二维DynRICE 模型（CRISSP）来模拟冰凌发展过程，并考虑了河床变化和泥沙运动，该系列模型已成功应用于世界多条河流的冰情研究中［25-26］。郭新蕾和杨开林等［27-28］在Shen 研究基础上开发了一维树状数学模型模拟调水工程等人工渠道冰盖发展过程，并成功应用到南水北调中线工程和京密引水渠中。对黄河冬季冰水情发展过程的模拟可追溯到1990年代，Shen 和黄河水利委员会合作，采用RICE 模型模拟黄河下游冰盖下水流变化和河道冲淤过程［29-30］。黄河水文局和Delft 大学合作，采用黄河水利委员会的一维黄河河冰动态模型（YRIDM）模拟2007—2008年三湖河口冬季流量和水位的变化［31］。但所开发的黄河冰情模拟系统很快就不能适应新的河道断面、水文和气象条件。

因此，本研究分析了黄河冰情演变规律和河段断面变化规律，解决河段计算中断面资料、部分气象和水文资料缺乏的难题，率定出黄河冰情计算中所需的关键参数。在此基础上进一步完善了冰情模拟软件系统，开展了黄河内蒙古河段水温、冰盖厚度、冰盖前沿、流量等冰凌发展过程的模拟和实测资料的对比验证。

2 河段断面和边界条件

内蒙古河段现有水文站5 个（如图1）：石嘴山、巴彦高勒、三湖河口、包头和头道拐，其中包头水文站2014年1月1日才开始运行，计算所用的水文资料来自上述水文站。黄河内蒙古河段最新河道大断面数据为2012年10月黄河水利委员会实测资料，从巴彦高勒到头道拐下游共测量166 个河道大断面。

图1 黄河内蒙古河段

2.1 河床断面作为游荡性河道，黄河河道随水流条件的改变断面形状变化较频繁。图2—4 为巴彦高勒站、三湖河口站和头道拐站2012—2013年度凌汛前后大断面套绘图，分别比较了流凌前、稳定封冻期和开河后断面的变化。巴彦高勒和头道拐断面凌汛后形状变化大，巴彦高勒河床最大淤积深度可达4 m，三湖河口主河床淤积整体超过1 m。三湖河口位居倒“U”字形河道底部的中部，河床断面变化不大。总之，内蒙古河道冬季流量小，过流能力低，河道呈现不同程度的淤积。黄河河道断面冲淤变化频繁且断面形状改变较大，在冰情演变过程计算中，当年的河道断面同本年度水文信息、冰情变化特征相联系。本研究采用2012年10月测量的河道断面资料，为了计算的准确性参数率定、冰情模拟均采用2012—2013年实测的水文和气象资料。

图2 巴彦高勒站2012—2013年度凌汛前后大断面套绘图

2.2 冰桥位置冰盖发展过程的模拟从冰桥起（Ice Bridging locations）向上游发展，也就是局部河段的首封位置向上游推进，2012年内蒙古河段最先封冻时间为11月30日，位于东经109°54′12″和北纬40°31′53″，该位置距离上游巴彦高勒水文站286.7 km，作为计算的一个冰桥。其下游仍需设冰桥位置，但是下游局部河段封冻起始位置未有完整记录，在计算中参考历史首封位置设定易卡塞形成冰桥的位置。

图3 三湖河口站2012—2013年度凌汛前后大断面套绘图

图4 头道拐站2012—2013年度凌汛前后大断面套绘图

从2000年有记录的首封位置（见图5）统计：2000—2019年首封主要集中在三湖河口水文站上、下游（3 次），包头水文站上游50 km 范围（其中磴口3 次，包头市3 次，包神铁路桥附近5 次），头道拐水文站上游什四分子大拐弯处（5 次），有的年份在同一天多处现首封，这些首封地点都是容易形成冰桥的位置。根据历年多次发生首封的位置和实测水文资料，2012—2013年度冰情模拟在首封位置以下河段另增加6 个冰桥。

图5 内蒙河段首封位置

2.3 支流边界黄河内蒙古河段水系发达，支流众多，其中一级支流有10 条。内蒙古境内十大支流，位于黄河河套内，十大支流几乎等距离从南到北流入黄河［32］。但是当前未有10 大支流相关的水文、气象资料，十大支流全年干旱少雨，降水主要以暴雨形式出现，多集中于夏季汛期，在黄河冰情模拟中不予考虑。

黄河三盛公水利枢纽工程通过南岸干渠、北总干渠和沈乌干渠引水服务于工农业用水，巴彦高勒和三湖河口之间常用的规模较大的退水闸分别为二闸、三闸和四闸。位于三湖河口上游的乌梁素海一方面承担内蒙古河段分洪任务，另一方面也需定期进行生态补水，因此巴彦高勒到头道拐河段计算中模拟的取水和退水边界为二闸、三闸、四闸和乌梁素海，位置见图6所示。

图6 巴彦高勒-三湖河口支流边界位置

三湖河口和头道拐之间存在多个扬水站和水厂用于工农业用水和居民饮水，规模较大的取水口有：打不素扬水站、包钢水源地取水口、画匠营子岸边取水泵房、磴口净水厂、团结渠电力扬水站、民利渠扬水站6 个支流边界，冬季封冻期只有包钢水源地取水记录，取水口位置如图7所示。

图7 三湖河口-头道拐支流边界位置

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

参考作者《资料缺乏区河渠冬季冰情过程模拟研究》文章，根据文中对缺少资料条件下支流流量的分配方法推演出支流流量分配系数。具体流量分配系数的确定通过调整各支流流量分配值以计算下游边界流量模拟值，利用最小二乘法寻优达到相对准确辨识区间支流流量分配比例，寻求下游流量计算较优值，确保计算河段流量的平衡。巴彦高勒和三湖河口两水文站的流量差按照0.1∶0.1∶0.3∶0.5 的比例分配到二闸、三闸、四闸和乌梁素海分水口；三湖河口和头道拐两水文站的流量差按照0.1∶0.1∶0.1∶0.1∶0.3∶0.3 的比例分配到打不素扬水站、包钢水源地取水口、画匠营子岸边取水泵房、磴口净水厂、团结渠电力扬水站、民利渠扬水站分水口。区间支流流量的合理分配确保了计算河段流量的平衡。

3 参数率定

计算断面资料依据2012年10月黄河水利委员会实测的巴彦高勒至头道拐河段下游共166 个实测断面基础上增加桥梁断面18 个、浮桥断面3 个，总计187 个大断面。实测断面间距较大会影响模型计算的稳定性及模拟结果的精度，因此按照最大断面间距不超过1 km，对断面进行线性插值，最终得到河道断面580 个，计算得到的是两个断面之间的间距为直线距离，没有充分考虑河道弯曲情况，得到河道长度约为486 km，落差68 m。

3.1 糙率准确的河道糙率是正确模拟河道过流能力的关键。钱宁通过水文站实测资料分析表明对于黄河这样含沙量高的冲积河道，糙率变化不仅同水力要素有关，还同泥沙因素有关［33］。根据大量实测资料推求，受多泥沙影响黄河河道糙率系数n 值存在较低值情况［34］，张防修在黄河洪水演进中取花园口河床糙率为0.012［35］，Fu 在宁蒙河段冰情计算中，分段率定出河床糙率取值区间为0.011-0.040［31］。建立在黄河河道走向、断面形状和滩地特点分析基础上，通过无冰期河道数值模拟率定河段的糙率为：计算河道0～28 km 糙率为0.040，28～196 km 糙率为0.012，196～207 km 糙率为0.030，207～443 km 糙率为0.012，443 km-末端糙率为0.015。

3.2 热交换系数气温和水温变化是凌情变化的决定因素，凌汛期是“降温-流凌-封河-升温-开河”的过程，通过水温模拟率定出热交换系数。

在一维条件下，沿流向的热扩散方程［19］为：

式中：ρ 为水的密度，kg/m3；Cp=4.2 为水的比热，J/（kg·℃）；A 为渠道断面面积，m2；Tw为断面平均水温，℃；V 为断面平均水流速度，m/s；Ex为热扩散系数，m2/s；B 为水表面宽度，m；为水体与周围环境的单位面积热交换率，包括：明流水面与大气的热交换率、水面与飘浮冰块或冰盖的热交换率、以及大气与冰盖的热交换率等。

水面与大气的热交换率采用线性传热函数［23］：

空气和水的热交换系数hwa的值在不同区域河流中存在差异，Lal 和Shen 在美国北部河流中取值19.7W/（m2·℃）［19］，Andrishak 在加拿大the Peace River 河取值15W/（m2·℃）［36］，Blackburn 和She 在阿拉斯加Susitna River 河取值为20W/（m2·℃）［23］。黄河内蒙古河段热交换系数hwa、jwa、kwa、hia通过冬季水温模拟值和实测值对比率定得到。计算中当输入太阳辐射和率定出合适的热交换系数hwa气温的模拟值和实测值吻合较好时，线性热交换系数 jwa和kwa取0。 hia的取值参考Shen 同纬度地区参考值而定的，hia=15［37］。通过多组次水温发展过程的模拟优化，选择系数 jwa=0、kwa=0、hia=15，α =1［29］时，三湖河口2012年11月6日到2013年3月20日（全线开河时间）水温实测值和模拟值更为接近的4 种工况（Case 1—Case 4）比较分析，图8 为工况Case 1、Case 2、Case 3 和Case 4 的热交换系数hwa分别取14、15、16 和17 的模拟值和实测值的比较，表1 采用确定性系数、均方根误差和绝对误差均值作为模拟值和实测值的评定指标。其中确定性系数越接近1，模拟值和实测值吻合越好；均方根误差和绝对误差均值越小，模拟值和实测值吻合越好。观察图8，4 种工况模拟值相差很小，且模拟值和实测值较为接近，但整体看模拟值都比实测值略微偏大。从表1 明显看出：随着热交换系数hwa的变化，hwa=15 时确定性系数0.52，均方根误差0.90℃，绝对误差均值0.72 ℃，模拟结果比其它3 种工况差。 hwa=14 和hwa=17 的结果接近，好于hwa=15 的模拟值。hwa=16 时整个冬季模拟确定性系数0.79，均方根误差0.60 ℃，绝对误差均值0.19 ℃，无论是确定性系数还是均方根误差、绝对误差均值比较都优于其它3 种工况，因此选择热交换系数hwa=16 用于冰情计算，从模拟值和实测值的评定指标和模拟曲线看，模拟值和实测值吻合较好。

图8 三湖河口水温计算

表1 水温模拟值和实测值误差的评定

3.3 弗劳德数Fr冰盖发展有3 种模式：平铺上溯模式、水力加厚模式和力学加厚模式。存在一个临界弗劳德数Frc，当冰盖前沿弗劳德数Fr 小于临界弗劳德数Frc（Fr

式中：V 为冰盖前缘上游水流的平均流速，m/s；H 为冰盖前缘的水深，m；ti为冰块厚度，m；li为冰块长度，m；为流冰的形状系数，取值在0.66～1.3 之间变化；e 为冰块孔隙率；ρ 为水体的密度，kg/m3；ρi为冰体的密度，kg/m3；g 重力加速度，m/s2。

当弗劳德数超过Frc时，冰块将出现翻转、下潜，冰盖将以水力加厚模式推进，这时冰盖初始厚度hice的计算公式可用Michel 计算［19］：

黄河内蒙古河段冰情计算中，由于结冰过程中的某些观测数据（冰花浓度、面冰密度、锚冰等）缺乏，无法对过程相关参数进行率定，结冰过程主要率定临界弗劳德数Frc和最大弗劳德数Frm，计算中取Frc=0.06，Frm根据计算河段的糙率、地形和走向等综合考虑分段率定，计算河段1～50 km取Frm=0.11，50～174 km 取Frm=0.13，174 km-末端取Frm=0.097。

4 黄河内蒙河段冰情的模拟

冰情演变是一个非常复杂的物理过程，其发展过程模拟的数学模型包括水流的热扩散方程、冰花的扩散方程、冰盖下水流的输冰能力方程、水面浮冰的输运方程、冰盖和冰块厚度的发展方程等。本研究采用加拿大Alberta 大学开发的开放冰情计算软件RIVER 1D 模拟水温、冰花浓度、岸冰的形成、冰花的输移及冰盖等要素的发展过程［23］，并根据黄河冰情的特点优化和改进了该软件。

4.1 冰盖厚度的模拟2012年11月6日—2013年3月20日三湖河口冰盖厚度发展过程模拟如图9所示，冰厚测量数据显示，2012年12月23日除部分清沟外研究河段全线封冻，2013年2月1日冰厚已发展达到最大值。三湖河口模拟式和实测值的均方根误差和绝对误差均值见表2，其中整个冬季、冰盖增长期（2012年11月6日—2013年2月1日）和冰盖消融期（2013年2月2—20日）模拟值和实测值的均方根误差分别是0.063 m、0.047 m 和0.079 m，绝对误差均值分别为0.055 m、0.046 m 和0.070 m，可以看出冰盖增长过程模拟误差比融冰期小。

图9 三湖河口冰厚的发展过程

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

三湖河口冰盖厚度模拟值和实测值的误差分析如下：（1）在冰厚增长期，热力作用为冰盖厚度增加的主要因素，而在冰盖消融期，冰盖并不是按照静止不动通过热力作用就地融化消失，通常开河阶段是由热力和动力共同作用，随着冰盖厚度的减小，冰盖强度也随之降低，当冰盖消融到一定厚度时，就会在水动力作用下破裂，向下游流动，但在一维冰情模拟中冰盖发展过程数学模型未能考虑动力因素的影响，所以冰盖厚度增长期模拟精度高于冰盖消融期。（2）冰盖厚度采用冰面凿冰孔后量冰尺深入冰盖下测量，固定点冰厚测量位置是否具有代表性、凿冰过程中冰盖结构破坏程度等不可控因素对冰盖测量的精度影响较大。但整个冬季和冰盖厚度增长期实测值和模拟值的均方根误差分别为0.063 m 和0.047 m，且模拟值和实测值发展趋势一致。

4.2 冰盖前沿的发展计算中设置7 个冰桥位置，从上游到下游依次为：①昭君坟上游附近、②包头和昭君坟之间、③包头附近、④五犋牛浮桥、⑤五犋牛浮桥下游、⑥头道拐上游和⑦下游边界（如图10所示），冰桥①为内蒙古河段首封位置，该位置距离上游巴彦高勒水文站286.7 km，2012年12月23日封冻上首发展到上游巴彦高勒，图11 为整个冬季冰桥①前沿发展过程模拟值和实测值的比较。从图上看冰盖前沿发展过程实测值和模拟值吻合较好，整个冬季模拟值和实测值的均方根误差为18.76 km，绝对误差均值6.43 km；开河前（2012/11/6-2012/3/5）模拟值和实测值的均方根误差仅为6.27 km，绝对误差均值2.50 km，模拟的冰盖前缘发展过程与观测值基本吻合，仍然存在开河期模拟值和实测值误差大于冰盖增长期和稳定期的情况，原因同4.1 节冰盖厚度模拟误差分析（1）。

图10 计算中冰桥的位置

图11 冰盖前沿的发展过程

4.3 流量的计算封冻后河道主槽过流能力明显下降，滩地的围垦，增加了河道滞水量，为槽蓄量增加提供了天然条件。内蒙古河段通常下游先封冻，封河期间局部冰塞冰坝的出现，也会导致过流能力降低，加剧槽蓄量增加。2012—2013年凌汛期槽蓄量变化如图12所示，最大槽蓄水增量出现在2月中旬，其中，巴彦高勒-头道拐河段最大值可达10.7 亿m3，巴彦高勒-三湖河口河段最大值可达7.5 亿m3，本年度内蒙古河段最大槽蓄水量较往年均值偏大12%。槽蓄量特点为：封河开始快速增加，稳定封河期缓慢变化，开河前后迅速释放，巨大的槽蓄量变化影响凌汛期流量模拟的精度。

表3 冰盖前沿发展模拟值和实测值的误差

三湖河口封河时间和开河时间分别为2012年12月5日和2013年3月14日，内蒙古河段全线开河时间为3月20日。根据流量实测值和模拟值比较（见图12）计算过程可分为4 个阶段：封河前期、冰盖快速发展期、稳定封河期和开河期，4 个阶段流量的模拟值和实测值均方根误差和绝对误差均值见表4，每个阶段流量变化特点为：（1）封河前期（11/1—12/4）河道内无槽蓄水量，流量的模拟值和实测值吻合好，这个阶段三湖河口流量模拟值和实测的均方根误差为55 m3/s，绝对误差均值为42 m3/s；（2）冰盖快速发展期（12/5—12/27）因冰盖阻力增加，槽蓄量快速增加，流量下降，河道中大量的槽蓄水导致流量模拟值和实测值差别较大，模拟值和实测值均方根误差为281 m3/s，绝对误差均值265 m3/s；（3）稳定封河期（12/28—3/5）冰盖厚度缓慢变化，冰盖下糙率减少，槽蓄量趋于稳定，这个阶段三湖河口流量实测值和模拟值吻合好，其均方根误差为65 m3/s，绝对误差均值为55 m3/s；（4）开河期（3/6—3/20）河道内巨量的槽蓄量在短短几天时间释放到下游，从图12 可以看出三湖河口槽蓄量集中释放期最大洪峰流量为1100 m3/s，这个阶段流量的模拟值和实测值误差大，其均方根误差为191 m3/s，绝对误差均值146 m3/s。

表4 流量模拟值和实测值的误差

图12 三湖河口流量模拟和内蒙古河段槽蓄量变化对比

4 结论

黄河冰情模拟中面临的主要困难是资料不足：缺少河道冬季大断面的信息，河道断面资料为冬季之前测量数据，断面之间间距大；河道支流流量、退水和取水资料缺乏；太阳辐射和气温资料不完整；受到测量手段限制，测量数据精度有待提高等。本研究分析了黄河内蒙古河段冰情演变规律，解决了黄河冰情模拟中冰桥确定方法和支流流量分配方案等关键问题。开展了不同工况下水温和冰情发展过程模拟的比较，率定出符合内蒙古河段的河道糙率、热交换系数、临界弗劳德数Frc、最大弗劳德数Frm等关键参数。本研究应用一维冰情数学模型模拟2012—2013年凌汛期黄河内蒙古河段冰情发展过程，其中对比了水温、冰盖厚度、冰盖前沿发展和流量的模拟值和实测值。2013年3月6日开河前水温、冰盖前沿发展和流量的模拟值和实测值的绝对误差均值分别为0.07 ℃、2.50 km、55 m3/s，冰盖消融前冰厚模拟值和实测值的绝对误差均值为0.046 m，尽量实测资料不够完善，但模拟值和实测值吻合较好。

本研究为天然江河、人工渠道冰情全过程计算提供了研究思路和模拟方法，并为黄河防凌减灾和冰情预报提供科学依据。随着冰情观测仪器的改进和观测技术的提高，观测数据种类、数量和质量将不断增加，黄河冰情模拟精度将会进一步提高。