摘要
开放渠道中的冰过程本质上是多物理现象，其特征是流体动力学、热传递和相变动力学的强烈耦合。尽管已经提出了许多模型来研究冰过程，但以往的研究往往受到流动简化表示、对渠道几何形状的假设以及对热和相变相互作用处理不当的限制。本文提出了一种双层平均数学模型，用于解决具有复杂横截面几何形状的开放渠道中耦合的水-冰-热动力学问题。该模型明确考虑了上层冰水混合物层与下层清水层之间的质量、动量和热交换，从而能够模拟相变过程和温度演变。该模型通过与实验室和现场数据的对比测试，证明了其与观察到的冰生长、衰减、温度分布和阶段流量图的满意一致性。应用于引水渠道的案例表明了该模型预测耦合水-热冰演变及其对不同环境和操作条件响应的能力。这种建模框架为在变化的水力和热条件下模拟开放渠道中的冰过程提供了一个有前景的工具。
通俗语言总结
河流和渠道中的冰形成可能会阻塞水流、损坏基础设施并扰乱供水。迄今为止，预测这些过程仍然是一个挑战，因为水的运动、热传递和相变以复杂的方式相互作用，而许多现有模型对此进行了简化。在这里，我们提出了一种新的双层平均模型，该模型包括一个上层冰水混合物和一个下层清水层，模拟质量、动量和热交换，以捕捉任何形状渠道中的冰生长和融化过程。通过与实验室、池塘、自然河流和大规模渠道观测数据的对比测试，证明了在水位、冰厚度和水温方面的满意一致性。该模型揭示了流量、界面粗糙度和渠道几何形状如何影响冰演变，为减轻冰塞和洪水风险提供了指导。
1 引言
开放渠道中的冰过程是多物理现象的典型例子，它们源于流体动力学、热传递和相变动力学的复杂耦合。这些相互作用的物理过程控制着河流、水库和水输送渠道中冰层的形成、生长和衰减，尤其是在中高纬度和高海拔地区。世界上超过一半的主要河流都会经历季节性的冰覆盖（X. Yang等人，2020年）。虽然冰覆盖可以支持某些经济活动，例如北极地区的冰路（Rojas等人，2011年），但其不利影响也是巨大的：冰塞会干扰航行并引发洪水，冰覆盖的演变会危及基础设施和河岸，而冰的存在会改变水生栖息地，并使水力发电和引水操作复杂化（Brown & Mackay，1995年；Carson等人，2011年；J. Wang等人，2013年）。这些现象的广泛影响凸显了深入理解和建模开放渠道中的多物理冰过程的重要性，以改善风险缓解、确保基础设施安全并支持有效的水资源管理。近几十年来，通过实验室实验（Daly & Axelson，1990年；Dow Ambtman等人，2011年；Healy & Hicks，2007年；Wazney等人，2019年）、现场观测（Beltaos & Burrell，2006年，2010年；Deng等人，2019年；Shen & Wang，1995年）和数值模拟（Beltaos，1993年；Shen，2010年；Shen等人，2000年；Zhai等人，2022年，2024年）对开放渠道中的冰过程进行了广泛研究。实验室研究通常集中在冰塞形成（Healy & Hicks，2006年；J. Wang等人，2011年，2019年）、流-冰相互作用（Healy等人，2002年；Robert & Tran，2012年；Tang & Davar，1985年）或热生长和衰减等特定方面，而不是全面考虑多物理冰过程。现场调查提供了宝贵的实际洞察，但受到恶劣天气条件下数据收集的挑战（Shen，2010年）。因此，能够整合流体动力学、热和相变过程的数值模拟已成为关键工具，特别是在数据稀缺或操作至关重要的情况下。早期的冰过程数学模型主要是单层模型。值得注意的模型包括River Ice Model（Lal & Shen，1991年；Shen等人，1995年）、MIKE-ICE（Timalsina等人，2013年）、RIVICE（Lindenschmidt，2017年）和River1D（Blackburn & She，2019年）。其中一些早期模型将冰覆盖表面视为纯水，从而忽略了冰的物理效应（Beltaos & Krishnappan，1982年；Blackburn & Hicks，2003年）。其他模型通过在流动方程中加入额外的阻力项来简化冰动量的表示（She & Hicks，2006年）。然而，观测显示开放渠道中的冰水流动通常是垂直分层的（Shen等人，2000年），包括一个上层冰水混合物层和一个下层清水层。为了解决冰水流动的垂直分层问题，开发了如DynaRICE模型（Shen等人，2000年）这样的双层模型。这些模型区分了上层冰水混合物层和下层清水层，提供了更符合物理实际的分层流动条件表示。值得注意的是，DynaRICE框架中引入了一个额外的面积冰质量守恒方程，类似于海冰和湖冰模型中使用的方法（Hibler，1979年；Wake & Rumer，1983年）。然而，这个额外的方程基于经验假设，并不是严格从基本的质量守恒原理推导出来的（Thorndike等人，1975年），限制了其理论一致性。Zhu等人（2025年）提出的双层平均冰模型是一个较新的发展，它建立在水-沉积物-形态动力学建模中建立的质量和动量守恒原理基础上（Cao等人，2015年；Li等人，2013年）。在这个框架中，上层被视为冰水混合物，而下层由清水组成。该模型基于深度平均的质量和动量守恒原理推导出，明确包含了层间质量和动量交换，并动态求解冰浓度的时空变化，提供了比假设恒定或经验相关冰浓度的模型更好的物理一致性和灵活性（She & Hicks，2006年；Zufelt & Ettema，2000年）。尽管Zhu的模型得到了有希望的验证，但它尚未完全解决热和相变动力学问题，并且仅限于矩形渠道。虽然对于短期事件（如冰塞释放）可以忽略热效应（She & Hicks，2006年；Zhu等人，2025年），但对冰覆盖形成和衰减的全面理解需要包括热传递和相变动力学。此外，现实世界的河流和渠道很少呈现矩形横截面，冰覆盖的演变受到横截面形状和纵向渠道地形的影响。本研究旨在开发一个改进的双层平均冰模型，该模型整合了具有复杂横截面的开放渠道中的流体动力学、冰、热和相变过程。在Zhu等人（2025年）的基础上，所提出的模型引入了温度演变，并考虑了不规则开放渠道中的多物理冰水相互作用。与基于Lagrangian SPH的DynaRICE模型不同，当前的Eulerian框架直接从深度平均守恒定律推导出冰演变，避免了需要一个关于冰面积的经验方程。此外，与使用经验冻结标准的传统1D模型（例如River1D）不同，该模型通过明确解决层间质量和动量交换来再现冰前沿的前进和增厚。该模型首先通过实验室实验、池塘中的现场观测和自然河流冰塞事件进行了验证，以证明其在捕捉冰覆盖生长、衰减、热演变和高度动态的水-冰涌动方面的能力。然后，它被应用于模拟引水项目中的大规模冰演变，进行了系列情景分析，以研究关键因素对冰动力学的影响，包括界面粗糙度、流量排放和渠道横截面形状。
2 数学模型
为了解决以往研究的局限性——特别是最近的双层平均模型（Zhu等人，2025年）中排除了热效应和仅限于矩形渠道几何形状的问题，本研究开发了一个改进的双层平均模型，用于模拟非矩形开放渠道中的耦合水-冰-热动力学。该模型明确表示了上层冰-水混合物层（包括固体冰、浆状冰和雾凇）与下层清水层之间的质量、动量和热交换。这两层由一个理论界面分隔，该界面定义在冰浓度接近零的深度。分别对每一层求解横截面平均温度方程。此外，该模型考虑了由于热交换导致的上层冰浓度变化，从而能够动态模拟相变过程。当下层完全冻结时，它被视为冰；融化会将水转移到下层，而冻结则从下层吸取水在上层形成新的冰。这种表述使得能够动态模拟冰形成-融化周期和相关的上层冰浓度变化。图1显示了模型变量和物理过程的示意图。
2.1 控制方程
在Zhu等人（2025年）提出的双层平均模型的基础上，本研究扩展了该模型，以纳入热效应并考虑非矩形渠道几何形状。控制方程是使用Xie（1990年）提出的含泥沙流动的广义浅水方程的方法推导出来的（详细推导见支持信息S1中的Text S1.1）。它们包括上层和下层的质量、动量和能量（温度）守恒定律，以及上层冰浓度质量守恒方程。这些方程捕捉了上层冰水混合物层和下层清水层的耦合动力学，并如下所示：
对于上层冰-水混合物流动层：
(1)
(2)
(3)
(4)
对于下层清水流动层：
(5)
(6)
(7)
其中t是时间；g是重力加速度；x是顺流方向坐标；Am和Aw分别是上层冰-水混合物流动层和下层清水流动层的横截面积；Qm和Qw分别是上层冰-水混合物流动层和下层清水流动层的流量；Tm是上层冰-水混合物的平均温度（当上层只有冰时，它代表冰盖的温度）；Tw是水温；Bw是界面宽度；B是表面宽度；Ew是两层之间界面的水吞蚀通量；ρw = 1,000 kg/m³和ρi = 920 kg/m³分别是水和冰的密度；ci是体积冰浓度；ρm = ρw(1 − ci) + ρici是冰水混合物的密度；Um = Qm/Am，Uw = Qw/Aw分别是上层冰-水混合物流动层和下层清水流动层的平均速度；η = ηin + hm是表面高程；ηin是下层与上层界面高程；hcm是上层横截面的质心高度；τin是上层和下层界面处的剪应力；τs是作用在上层的边界剪应力（由于渠道岸坡或河床）；τb是作用在下层的床面剪应力；Cp是水的比热；Cpi是冰的比热；Φ代表冰水相变率；Φ* = Φ仅在上层完全是冰时适用（ci = 1）：下层水的冻结有助于上层冰盖的生长，而冰盖的融化则为下层提供水；否则，这个项是不活跃的，即Φ* = 0；Cf代表在下层生成的雾凇冰的通量，仅在过冷条件下激活，代表向上层的水质量交换；ϕwa、ϕia、ϕiw分别是上层水与空气、下层水与空气以及上下层之间的单位表面积净热交换率。方程1表示了冰水混合物流动层的质量守恒方程。其右侧的第二项和第三项量化了上层冰水相变引起的体积变化，最后一项表示从下层交换的雾凇冰的体积通量。方程2表示了冰水混合物流动层的动量守恒。右侧第一项包括作用在该层上的重力和静水压力；第二项代表界面摩擦力和河岸阻力；第三项和第四项描述了两个运动层之间清水交换引起的动量传递；第五项表示冰浓度空间变化的影响；第六项和第七项代表冰-水相变的作用；最后两项表示由下层生成的浮冰向上运动引起的动量交换。方程3控制上层中冰的质量守恒，同时考虑了由热交换和浮冰流入驱动的相变。方程4和7描述了两层中的温度变化，考虑了热对流和层间热交换，其中方程4右侧的前两项与上层中的清水有关。方程5和6控制清水层中的质量和动量守恒。在方程6中，第一项表示重力和静水压力；第二项描述界面摩擦力和河床阻力；第三项考虑了冰-水混合物层对下层施加的静水压力；第四项表示通过界面进行的清水交换引起的动量传递；最后一项描述了由下层生成的浮冰向上运动引起的动量损失。虽然双层平均框架提供了现实的宏观描述，但它假设每层内部的条件是均匀的。目前，该模型没有规定明确的定量阈值（例如，就冰浓度、弗劳德数或湍流强度而言），超出这些阈值时这种双层描述就不再适用。因此，其适用性取决于冰层是否可以合理地被视为一个垂直平均的连续体。从这个意义上说，双层模型旨在表示冰-水系统的整体行为，而不是解决冰层内部的详细垂直结构。因此，它最适合于冰层表现为相对均匀的连续体的情况，如冰的传输、积累和早期冻结阶段。当更复杂的垂直结构发展起来时，例如在底层泥浆层上的固化冰盖或强机械分层时，现有的双层公式无法明确解决这些次层特征，尽管它们的整体效应仍然可以被表示出来。
2.2 模型闭合
为了闭合当前双层平均模型的控制方程，必须引入一组关系来确定τb、τin、τs、Ew、Cf、Φ、ϕwa、ϕia和ϕiw。尽管Manning方程最初是为稳定均匀流动推导出来的，但在数值建模中，将其应用于非稳定分层系统通常被视为标准的准稳态近似，前提是计算时间步长足够小，能够分辨出流动变化。这种方法已经在之前的双层平均水沙流模型中成功采用并得到了验证（例如，Cao等人，2015年；Li等人，2013年）。根据这些研究，界面和河床剪切应力使用Manning粗糙度方程来估计：
其中nin是上下层之间的界面粗糙度；nb是河床粗糙度；Rm和Rw分别是上层冰-水混合物流层和下层清水流层的水力半径；χb = Am/Rm − Bw是上层受河床影响的湿润周长。使用整体水力半径不可避免地会平均河床剪切 stress 的横向变化。虽然在本研究中采用均匀的粗糙度nb是对相对规则通道的一个合理近似，但可以使用横向积分方法（Cao等人，2006年）来扩展该模型，以适应具有明显横向变异性的自然河流。此外，界面粗糙度nin代表由两层之间的剪切相互作用和形态阻力引起的有效阻力。尽管在流体-流体界面使用Manning型方程是一种宏观近似，但在之前的分层流和双层平均模型中也采用了类似的处理方法（例如，Cao等人，2015年）。实际上，nin 选择一个物理上合理的范围（通常为0.01–0.05 m−1/3 s），并通过与现有观测数据的比较进一步细化。根据Zhu等人（2025年）的方法，水夹带质量通量Ew由（Parker等人，1986年）计算：
其中Richardson数。虽然这种关系最初是为含沙流动开发的，但将这种基于沙输送的公式应用于描述冰动力学是河冰建模中的成熟做法（例如，Blackburn & She，2019年；Shen & Wang，1995年）。由于缺乏针对冰-水混合物的成熟闭合公式，当前的公式为表示界面质量交换提供了一个实用的近似。为了处理将这种基于沙的闭合应用于冰-水混合物可能带来的不确定性，我们评估了模型结果对夹带通量的敏感性（详见支持信息S1中的文本S3.2），评估了其对冰输送行为的影响。当下层水体过冷（Tw < 0°C）时，将会生成浮冰并将其交换到上层。水体中浮冰的生成率由（Lal & Shen，1991年）表示：
其中λt表示从过冷水到浮冰的相变特征时间尺度。在当前模型中，采用了瞬时平衡假设，意味着任何过冷都在一个计算时间步长内有效消除（Beltaos，2013年；Liu等人，2025年）。因此，假设生成的浮冰会由于浮力立即转移到上层，同时将下层水的温度恢复到冰点。使用热交换率计算冰-水相变率由界面处的热平衡决定：
其中Li = 334 kJ/kg是冰的潜热。在这个模型中，热通量使用Andrishak和Hicks（2008年）提出的线性热传递模型来计算，所有热通量以W/m2为单位。水和空气之间的热交换ϕwa由以下公式给出：
其中αwa是水面的线性热传递系数，Ta是空气温度，ϕs是入射的太阳辐射。上层和下层之间的热交换ϕiw按（Blackburn & She，2019年；J. Yang等人，2023年）计算：
其中是界面的线性热传递系数（Ashton，1973年），采用浓度加权方法应用；Tf是界面温度，可以设为0°C。由于之前的开放式渠道双层平均模型没有明确包括界面热交换，我们采用了这种已建立且经过验证的公式来模拟两层之间的热通量。αiw也直接与水的热导率成正比，与水深成反比，在静态水条件下适用。它适用于静止覆盖层和移动的冰层，确保在冻结过程中热力过渡的平滑性。水和空气之间通过上层的热交换ϕia表示为：
2.3 数值算法
两层的控制方程系分别同步求解。每个系统都写成矩阵形式（Li等人，2013年；Zhu等人，2025年），并使用适用于一维开放式渠道流动的拟黎曼求解器进行求解，并采用湿-干处理。根据Jiang等人（2024年）的方法，控制方程组可以写为：
其中U、T是变量向量；F、G是通量的对流通量向量；Sb和Rb、Sf和Rf分别是与坡度和摩擦相关的源项；S0、R0是其他源项：
显式有限体积离散化如下：
其中Δt是时间步长，Δx是空间步长，i是空间节点索引，j是时间节点索引，Fi+1/2、Gi+1/2是在x = xi+1/2处的单元间数值通量。为了数值稳定性，时间步长满足Courant–Friedrichs–Lewy条件：
其中λmax是从雅可比矩阵计算出的最大速度，Cr是Courant数。为了确保数值稳定性，本文使用Cr = 0.5。本研究采用了显式方案。这一选择与之前的双层平均模型一致，并为解决两层之间的质量、动量和热交换提供了直接和稳健的框架。在当前工作中，主要目标 là 建立水-冰-热双层平均框架，而不是优化离散化方案。值得注意的是，半隐式方案在计算效率方面可能有优势，将来可能会考虑将其纳入当前的双层平均框架中。整体数值过程在图2中示意性地概述。采用Harten-Lax-van Leer（HLL）近似黎曼求解器来求解流动方程（Ying & Wang，2008年）。对于其他方程，包括温度和浓度方程，使用了Slope Limiter Centered（SLIC）近似黎曼求解器以及MINMOD坡度限制器（Jiang等人，2024年）。
3 基准测试
3.1 热过程：冰盖的生长和衰减
当前模型的热部分包括上层和下层的横截面平均温度方程，使得能够模拟冰覆盖渠道中的热过程。重新审视了两个测试案例来评估模型的热性能：(a) 在受控热强迫下跟踪冰-水温度演变的实验室实验；(b) 对自然河道中季节性冰盖动态的现场观察。这两个测试案例都涉及静态水条件和完全固态冰覆盖的时期。在这些条件下，第2节描述的全部控制方程可以通过将流速设为零并将冰浓度设为1来简化为一组常微分方程，从而得到此处使用的简化模型。这些简化方程保留了冰层和水层之间的关键热交换过程，使模型能够捕捉冰的生长和衰减，并确保充分代表了这些场景中的主导物理过程。简化模型在支持信息S1的文本S1.2中给出。
3.1.1 实验室实验
进行了一系列实验室实验来研究冰盖生长和衰减的热过程（Hao等人，2009年；Teng等人，2011年）。在Hao等人（2009年）进行的实验中，使用了一个可调节的温度冷却箱作为冷却装置，而一个具有非常低热导率系数的泡沫箱作为水的容器（见支持信息S1中的图S1）。冷却箱由泡沫板绝缘，泡沫箱之间的缝隙填充了具有高比热容和低热导率的土壤。这种设计有效地减少了水向周围环境的热量损失，因此可以假设热交换仅发生在水和空气之间。在受控的冷却条件下，冰层厚度在水面上基本保持均匀；因此，在建模中忽略了其空间变化。分析了两个具有不同初始条件的实验案例（表1）。在这两种情况下，泡沫箱的初始水深均为22.4厘米。在冰盖生长期间，空气温度在案例1中从-14.0°C降至-16.1°C，在案例2中从-14.5°C降至-17.0°C。由于两种情况下的平均室温分别为19.8°C和20.1°C，因此发生了冰盖的衰减。假设冰层（上层）中的温度分布是线性的，模型计算仅考虑了上层和下层的平均温度和厚度。表1. 实验的初始和边界条件
3.1.1 实验室实验
进行了一系列实验室实验来研究冰盖的生长和衰减的热过程（Hao等人，2009年；Teng等人，2011年）。在Hao等人（2009年）进行的实验中，使用了一个可调的温度冷却箱作为冷却装置，而一个具有非常低热导率系数的泡沫箱作为水的容器（见支持信息S1中的图S1）。冷却箱由泡沫板绝缘，泡沫箱之间的缝隙填充了具有高比热容和低热导率的土壤。这种设计有效地减少了水向周围环境的热量损失，因此可以假设热交换仅发生在水和空气之间。在受控的冷却条件下，冰层厚度在水面上基本保持均匀；因此，在建模中忽略了其空间变化。分析了两个具有不同初始条件的实验案例（表1）。在两种情况下，泡沫箱的初始水深均为22.4厘米。在冰盖生长期间，空气温度在案例1中范围为-14.0°C至-16.1°C，在案例2中范围为-14.5°C至-17.0°C。由于两种情况下的平均室温分别为19.8°C和20.1°C，因此发生了冰盖的衰减。假设冰层（上层）中的温度分布是线性的，模型计算仅考虑了上层和下层的平均温度和厚度。
图3a显示了两种情况下计算出的和测量到的冰盖厚度随时间变化的比较。模型成功捕捉到了冰盖生长和衰减的过程，这与测量数据的一致性得到了证明。图3b显示了两种情况下计算得出的水和冰的温度，以及案例1中测量得到的冰盖温度。需要注意的是，大多数温度探头都设置在靠近水面的位置（详见支持信息S1中的图S1），因此对于较深区域的测量数据较少。此外，水温度的非线性垂直分布使得难以对温度数据进行可靠的深度平均处理。与早期的冰模型（如RICE）不同，当前模型能够准确解析冰层（上层）的平均温度，并且与测量结果吻合得很好。研究结果表明，在冰层消融阶段，冰的温度首先升高到熔点，然后保持恒定，直到冰完全融化。
3.1.2 在池塘中的现场观测
与第3.1节中的实验室实验不同，该模型还使用了从中国内蒙古的一个池塘收集的现场观测数据进行了进一步测试（Ji等人，2016年）。与实验室中恒定的空气温度条件不同，现场环境中的冰盖会暴露在不同的温度变化下，提供了更为真实的验证环境。这个池塘的周长为306米，面积为5,673平方米，平均深度为2.46米。选定了五个测量点，水深分别为20厘米、30厘米、50厘米、70厘米和120厘米，详见支持信息S1中的图S2。从2014年11月23日到2015年3月23日，每天在上午8点和下午2点进行两次测量，涵盖了冰层的生长和消融阶段。相应的空气温度数据见支持信息S1中的图S3。图4展示了两个阶段中各位置的计算值和测量值。该模型能够以合理的精度捕捉到观察到的冰层生长和消融模式，尤其是在水深超过20厘米的区域。然而，在靠近岸边的地方，由于水深非常浅，偏差较大。Ji等人（2016年）也使用了一个能够解析深度的模型来模拟这些过程。表2总结了两种模型的平均误差。值得注意的是，这两种模型都没有考虑冰层厚度的空间变化性，也没有考虑风和太阳辐射的影响。
3.2 动态过程：阿萨巴斯卡河的冰坝释放
2002年加拿大阿萨巴斯卡河发生的冰坝释放事件被用来验证该模型捕捉自然河道中高度动态机械过程的能力（详见支持信息S1中的文本S2.3）。这一由机械因素驱动的事件涉及一个快速传播的涌浪，其影响范围约为100公里，提供了一个专注于水冰动力学的验证机会，同时排除了热过程的影响。由于缺乏详细的原始数据，本研究和Zhu等人（2025年）中的河道几何形状都是基于She和Hicks（2006年）的描述来估算的。然而，一个关键的区别在于几何表示方式：Zhu等人（2025年）将河道近似为一个宽度恒定为500米的棱柱形矩形通道，而当前模型则将河道表示为一个非棱柱形的矩形截面，其宽度从上游边界（x=380公里）的500米逐渐增加到下游边界（x=280公里）的700米。这种配置更准确地反映了She和Hicks（2006年）描述的河道加宽现象。为了确保比较的一致性，当前模型中的所有物理参数都与Zhu等人（2025年）使用的参数完全相同。此外，由于数据限制，两个模型中的模拟都在冰层到达下游边界时终止。为了量化数值解与测量数据之间的差异，定义了一个无量纲误差指标，具体公式如下：
图4展示了两种模型计算的冰层厚度与实际测量值的对比。第一组数据（左侧面板）展示了冰层的增长过程，第二组数据（右侧面板）展示了冰层的消融过程。表2总结了两种模型计算和测量得到的冰层厚度值的平均误差。
4 实地案例研究
4.1 研究案例
所提出的模型被应用于中国北方一个大型引水项目中的人工渠道。这项研究聚焦于一个位于水力控制结构（涵洞）上游的传统梯形混凝土衬砌渠道段（见图6）。该渠道的底部宽度为8米，边坡比为1:2.5，因此在典型的冬季流动条件下，水面宽度约为27米。作为引水网络中最北端的开放式渠道部分，这个渠道经常受到严酷冬季条件的影响，容易形成大量的冰层。在极端寒冷事件中，浮动的冰块会在上游积聚，形成冰坝，从而减小有效流通面积。这会导致上游水位升高，对结构安全造成威胁，并干扰向下游地区的计划供水。为了监测和分析冰情，在研究现场建立了一个冰情监测站。长期观测表明，冬季的冰层覆盖可以持续长达50天，记录到的最大冰层厚度达到30厘米。在本研究中，使用所提出的模型模拟了两次典型的冰情事件，以评估其在实际现场条件下的表现。
4.2 事件1：2019–2020年冬季的冰层生长与消融
在2019–2020年的冬季，由于下游泵站及相关管道正在进行维护，渠道必须在静态水条件下运行。水深保持在3.76米，初始温度为4°C。冰期从2019年12月3日持续到2020年2月29日，共持续了89天。冰层在2019年12月7日开始形成，并在2020年2月26日完全融化。在这方面，冰层的演变与第3.1节中的案例类似，因此也使用了简化版本的模型。该模型利用现场监测站测量的空气温度来模拟冰层的生长和消融过程。图7显示，计算结果与实际测量的冰层厚度演变非常吻合。冰层的发育可以分为三个阶段：（a）2019年12月11日至2020年1月22日的初期生长阶段，（b）2020年1月23日至2月10日的稳定阶段（仅有小幅波动），以及（c）2020年2月11日开始的融化阶段。
4.3 事件2：2021年1月的水冰-热动力过程
2021年1月6日，水温的急剧下降触发了冰层的快速形成。到下午6点30分，监测站观察到大面积的表面冰层，覆盖了超过90%的渠道宽度。这些冰块的厚度为2–3毫米，移动速度与水流速度相当，并迅速在冰坝前积聚，形成了初始的冰层覆盖，并逐渐向上游扩展。冰层一直持续到2021年1月17日融化。为了模拟这一事件，使用了完整的模型（即方程1-7）。这一时期的流量Qw约为24立方米/秒，因此被假设为恒定值。同时假设下游水位保持不变。由于存在倒置虹吸效应，没有考虑上游的冰流入。边界温度是根据每日测量值设定的。模型中的关键水力和热参数，以及网格设置详见表4。模型参数包括河床坡度为1/24,000，混凝土河床的曼宁粗糙度系数nb为0.015米^-1/3秒，界面粗糙度nin为0.010米^-1/3秒。热传递系数αwa根据第4.2节的测量结果校定为17瓦特/平方米/摄氏度。这个值处于类似气候条件下的典型范围（Andrishak & Hicks, 2008; Blackburn & She, 2019）内。流体变率系数βiw直接参考现有文献（Andrishak & Hicks, 2008; Ashton, 1973; Blackburn & She, 2019）设定为1,187Ws0.8/米^2.6/摄氏度。需要指出的是，αwa是一个经验系数，用于表征空气-水的热交换，因此会引入不确定性。特别是，用单一的线性经验系数来表示复杂的热通量过程不可避免地会导致不确定性，因为未考虑气象因素（例如风速和云层覆盖）。从定量角度看，当前模型的平均误差为18.23%，而Zhu等人（2025年）的模型为22.17%（表3）。这种改进在物理上是由于采用了线性变化的河道宽度，使得模型能够考虑河道存储和水力阻力的纵向变化，而这些在棱柱形近似中被忽略了。这种对流场面积和水力半径的精确表示提高了波浪速度和扩散衰减的预测精度，从而减少了对水位峰值和波浪速度的过高估计。此外，当前模型的表现与She和Hicks（2006年）的模型相当，甚至略好（18.92%），表明所提出的非棱柱形矩形表示方法能够捕捉自然河流的基本水动力响应，同时保持了冰-水相互作用过程的简化几何描述。
图5展示了当前模型、She和Hicks（2006年）模型以及Zhu等人（2025年）模型预测的阿萨巴斯卡河冰坝释放事件的水位，以及实际测量数据：（a）G140站点，（b）G135站点，（c）G130站点。表3列出了2002年阿萨巴斯卡河事件的各模型预测值（单位：%）。
4.1 研究案例
所提出的模型被应用于中国北方一个大型引水项目中的人工渠道。该研究关注的是一个位于水力控制结构（涵洞）上游的传统梯形混凝土衬砌渠道段（见图6）。在典型的冬季流动条件下，渠道的底部宽度为8米，边坡比为1:2.5，因此水面宽度约为27米。作为引水网络中最北端的开放式渠道部分，这个渠道经常受到严寒天气的影响，容易形成大量的冰层。在极端寒冷事件中，漂浮的冰块会在上游积聚，形成冰坝，减少有效流通面积。这会导致上游水位升高，对结构安全构成威胁，并干扰对下游地区的计划供水。这种中断可能严重影响供水效率，对下游的生产和日常生活产生潜在影响。为了监测和分析冰情，在研究现场建立了一个冰情监测站。长期观测表明，冬季的冰层覆盖可以持续长达50天，记录到的最大冰层厚度达到30厘米。在本研究中，使用所提出的模型模拟了两次典型的冰情事件，以评估其在实际现场条件下的表现。事件2的相关条件及主要参数设置
符号 | 值（单位） |
------------------------|-------------------------------------------|
Qw | 某时间段内的排水量 | 24（m³/s） |
Sb | 床面坡度 | 1/24,000 （–） |
nb | 混凝土床面的曼宁粗糙系数 | 0.015 （m⁻¹/³s） |
nin | 界面处的曼宁粗糙系数 | 0.010 （m⁻¹/³s） |
αwa | 水面的热传递系数 | 17 （W/m²/°C） |
βiw | 界面的热交换系数 | 1,187 （Ws⁰.⁸/m².⁶/°C） |
ki | 冰的热导率 | 2.22 （W/m/°C） |
Δx | 网格长度 | 2 （m） |
支持信息S1中的图S5展示了在初始冰层形成后的四个时间点，三种不同网格（Δx = 1、2和4 m）下的水面和冰层剖面。细网格（1 m）和中等网格（2 m）的计算结果相似，而粗网格（4 m）的结果则有所不同。两种较细网格之间的差异几乎无法区分；因此，采用了2 m的网格以确保解决方案与网格尺寸无关。为了进行比较，还应用了阿尔伯塔大学开发的River1D模型来模拟这一事件。为了定量评估两种模型的性能，基于L1范数定义了冰盖厚度和水温的差异度量指标：
（26）
（27）
其中， 和 分别是冰盖厚度和水温的L1范数； 和 分别是实测的冰盖厚度和水温，而hm、Tw是由数值模型计算出的冰盖厚度和水温。水温的绝对L1范数值可能显得较大，主要是因为大多数观测温度接近零，只有少数非零值被记录下来。为了在这种情况下提供更有意义的评估，引入了相对平均偏差（RMD）作为额外的指标：
（28）
其中， 是实测的平均冰盖厚度，n代表测量数据的数量。
4.3.1 冰盖厚度和水温
图8a展示了2021年1月7日至1月17日期间监测站计算和实测的冰盖厚度。现场观察显示，1月7日开始形成冰盖，初始厚度为4.5厘米，在持续寒冷的条件下，到1月11日冰盖迅速增长到大约15厘米。在最寒冷的时期，冰盖相对稳定，1月13日之后随着气温升高，冰盖开始融化，并最终在1月17日消失。当前模型成功捕捉了这些过程，与观测到的冰盖厚度变化情况吻合良好。
图8a显示了2021年1月冰盖形成后的模型结果：(a) 计算和实测的冰盖厚度比较；(b) 计算和实测的水温比较。模型合理地再现了水温达到0°C时的冻结过程，以及随后的冰层冷却过程。值得注意的是，现场数据或模型计算中均未观察到显著的过冷现象。随着气温上升，水温逐渐升高，冰温也升至0°C，随后冰盖开始融化。模型有效地再现了冰盖的热增长和融化阶段。值得注意的是，在冰盖融化阶段，模型相比增长阶段的误差较大。在模拟情景中，冰层演变主要由水力传输和热增长驱动，内部应力是次要机制。因此，模型准确捕捉了热过程驱动冰盖厚度增加的生长阶段。然而，在自然渠道中，冰盖的消融阶段通常由热融化和机械破碎机制共同作用。总体而言，冰盖的恶化很少仅通过原位融化发生；相反，如破碎和传输等动态过程起着重要作用。目前的一维模型在捕捉这些复杂的机械动态方面存在局限性，因此仅考虑了冰盖的热融化。由于缺乏这些关键的机械破碎过程，模型在消融阶段的计算准确性相对较低。为了定量评估模型性能，表5使用了之前提供的定义计算了差异度量指标。当前模型在冰盖厚度上的L1范数较低（3.37%），而River1D为6.78%，表明其在预测冰盖生长和融化方面的性能更好。对于水温，两种模型的L1范数相似（当前模型为31.92%，River1D为31.63%）。然而，如前所述，这些值较大是因为测量温度的幅度较小，大多数值接近零。为了更好地解释温度预测性能，还报告了相对平均偏差（RMD），显示当前模型的RMD较低（16.28%，而River1D为20.51%），反映了与实测温度数据的一致性改进。
图9进一步对比了两种模型在计算稳定冰盖形成时的性能。现有的单层模型（如River1D）可以再现几个宏观特征，包括平衡冰盖厚度和水面剖面。然而，像积累和停止这样的冰过程通常通过预设的经验准则（例如，局部弗劳德数阈值）来表示。相比之下，当前的双层平均框架直接求解了冰-水混合物层和下方流动的质量和动量方程。因此，冰的积累、停止和前沿推进是直接从控制方程中得出的，而不是基于预设的阈值。这提供了一种更基于物理的冰流相互作用表示，并能够明确模拟冰趾形成的动态，如图9所示。
4.3.2 冰盖形成动态
2021年1月6日的现场观察到，表面冰层迅速形成，并在涵洞入口处的冰坝前快速积累，从而开始形成冰盖。冰前沿稳定地向上游推进，最终形成了稳定的冰盖。图10显示了这一时期冰前沿位置和相应弗劳德数的计算演变。时间t = 0小时对应于1月6日的18:30。蓝色虚线表示通常引用的经验阈值（Fr = 0.06），用于区分冰层并置和机械增厚。需要注意的是，当前模型并不使用这个经验阈值作为输入参数；相反，计算出的Fr自然稳定在这个值以下。洋红色点划线标记了实测的最大Fr值（0.063）。模型结果显示，到1月7日，冰前沿已经向上游延伸了大约320米，与现场观察报告的稳定冰盖长度约300米相符。冰前沿的弗劳德数在冰层开始积累时起初波动，达到了大约0.063的最大值。这一结果与观测估计相符，证明了模型在解决与冰盖推进相关的流动动态方面的能力。
4.3.3 界面粗糙度对冰动态的影响
在该模型中，界面粗糙度表征了表面冰层与下方水体之间的相互作用，影响了速度剪切和冰积累动态。这一特性使当前模型与早期依赖简化冰传输假设的方法区分开来。图11展示了不同界面曼宁粗糙度值下的计算水面和冰层剖面。结果表明，较小的粗糙度值减少了冰-水混合物层与清水层之间的速度剪切，从而加快了冰层的传输和在冰坝前的积累速度。这导致冰前沿的推进更快，上游水位上升得较少。相比之下，较大的粗糙度值增加了界面处的速度剪切，降低了表面冰层的传输速度，减缓了冰前沿的推进速度。这导致冰坝附近的冰层积累更厚，上游水位上升得更明显，如图12所示。这些结果突显了界面粗糙度在调节冰盖动态和水位变化中的重要作用。理解这种关系对于制定操作策略很有帮助，特别是在寒冷季节控制水流和管理与冰相关的风险方面。
4.3.4 流量对动态的影响
流量是调节渠道中冰盖动态的主要控制参数。为了评估这一影响，进行了在不同流量场景下的一系列模拟。图13展示了这些场景下的计算水面和冰层剖面，而图14展示了相应的冰前沿推进情况。结果表明，较高的流量增强了冰坝上游表面冰层的传输和积累，导致冰盖更厚。在高流量条件下，增加的流速还带动了更多冰层进入下方水体，减缓了冰前沿的推进速度。这种水力增厚效应导致上游水位显著上升，并放大了回水效应。相比之下，较低的流量促进了通过并置过程形成较薄但更广阔的冰盖。在这些条件下，冰前沿更快地向上游推进，相关的水位变化显著减少，如图13和14所示。这些发现表明，低流量操作有利于稳定冰盖的形成，同时最小化水位波动，为冰季期间的渠道管理提供了实际指导。
4.3.5 渠道横截面形状对动态的影响
开放式渠道的几何配置，特别是横截面形状，在控制冰盖演变和相关水力过程中起着关键作用。为了研究这些效应，使用之前描述的事件2场景（第4.3-4.3.4节）进行了额外的数值实验，但将原始的梯形横截面替换为复合横截面。这两种配置都在相同的初始水深下构建，以便直接比较不同横截面几何形状对冰盖动态的影响（详见支持信息S1中的文本S3.3）。图15展示了冻结后的两种横截面配置的计算水面高度和冰盖剖面。在相同的初始水位和流量条件下，复合横截面的水面宽度（29.3米）比梯形截面（27.1米）更大。这种更宽的表面有助于更有效地容纳冰盖的生长，从而在冰积累期间减少了水位的上升。此外，复合横截面中的冰前沿推进速度更快，形成了沿渠道分布更薄但更广阔的冰盖。这些发现与低流量场景下的结果一致（第4.3.4节），其中降低的流量和增加的渠道宽度促进了稳定、薄冰盖的发展，并减少了上游水位的波动。结果表明，独立于流量的横截面几何形状可以显著影响冰盖动态和渠道系统在冰形成事件中的水力响应。从操作角度来看，这些结果表明，具有更大表面宽度的渠道设计可能在缓解不利的水位上升和促进冻结期间更可预测的冰盖推进方面具有优势。然而，需要指出的是，建议进行进一步的现场调查和敏感性分析，以评估在不同河道几何形状和水文条件下这些发现的稳健性。图15展示了冻结后不同横截面的计算水面和冰层剖面。
4.3.6 热过程对动态的影响
本研究通过引入热过程，扩展了双层平均冰水混合物流动模型，这是模拟开阔河道中冰盖完整演变的关键改进。图16比较了有无热耦合情况下的计算水面和冰层剖面。当热耦合被关闭时，上层冰水混合物中的冰仅来源于初始的冰块。上层保持在积累状态，在模拟过程中不会发展成稳定的冰盖。在这种情况下，冰坝上游的水面冰层会迅速堆积，在2小时内形成冰塞，从而抬高下游水位。在稳定流入流量且没有热耦合的情况下，一旦形成冰塞，上层就会停止演变——这种行为与观察到的物理过程不符。相比之下，具有活跃热耦合的模拟能够捕捉到驱动冰生长和消融的冰-水-热相互作用，从而更真实地再现冰盖的发展过程。这一比较强调了在开放河道中准确模拟冰盖演变时包含热过程的必要性。
5 结论
本研究提出了一个用于开阔河道中冰-水-热动态的双层平均模型，明确表示了上层冰水混合物层与下层清水层之间的质量、动量和热量交换。该模型通过引入热效应并将其适用范围扩展到非矩形和纵向变化的河道几何形状，改进了现有方法。通过与实验室和野外池塘数据的验证，证明了该模型能够解析冰盖的水-冰-热演变过程；而应用于冰塞释放过程的应用则证实了其再现动态机械过程的能力。考虑河道宽度的纵向变化通过更好地表示河道储存能力和水力阻力，提高了水文图预测的准确性。将该模型应用于导水渠进一步证明了其捕捉冰盖形成和消融的时空演变的能力。结果强调了热过程在控制冰盖演变中的重要性，以及流量、界面粗糙度和河道几何形状对冰积累和传播的联合影响。较低的流量和较宽的河道往往会产生较薄但更稳定的冰盖，而增加的速度剪切会影响冰盖发展的速度和范围。与River1D模型相比，所提出的模型在解析热动态和预测关键特征（如冰趾位置）方面表现得更好，这是由于其明确处理了层间动量交换。当前模型主要关注热主导的消融过程，没有包括冰块破碎等机械破坏过程，这可能会在消融后期增加不确定性。界面质量和热量交换的参数化也存在不确定性。未来的工作将集中在精细化这些过程，并结合额外的现场测量数据，以在更广泛的条件下进一步评估模型性能。
致谢
本文报告的研究工作得到了国家自然科学基金（项目编号52239007）和国家重点研发计划（项目编号2024YFF1306904）的财政支持。
作者声明
作者声明与本研究无关的任何利益冲突。
数据可用性声明
用于生成正文中图3-5和7-16以及补充信息S1中的图S3、S5和S6的所有数据和代码均可在以下链接获取：https://doi.org/10.5281/zenodo.16885126 (H. Wang, 2026)。
