河流水冰沙耦合模型研究Ⅰ：原理和方法

潘佳佳，Hung Tao Shen，郭新蕾，王 涛

（1.流域水循环模拟与调控国家重点实验室，中国水利水电科学研究院，北京100038；2.Department of Civil and Environmental Engineering，Clarkson University，Potsdam，NY 13699-5710）

1 研究背景

我国北方河流如黄河、黑龙江和松花江等每年都有超过100天的冰期［1］，而冬季河冰运动对泥沙输运和河道演变的影响常常被忽略。一方面受地区和时间的限制，河冰影响下的水沙问题是季节性过程，不及明渠水沙研究更具代表性，常常被研究人员忽略［2-4］。另一方面冰期河流涉及冰体堆积释放、水位壅高、流量波动、河床冲淤变化和岸滩崩塌侵蚀等多种过程，存在复杂的水冰沙相互作用。北方河流水冰沙耦合作用机理问题是水力学、河冰动力学和河流动力学的交叉方向，涉及的物理因素多［4-6］，问题复杂，是河冰领域研究的前沿和难点。

冬季河冰过程对北方河流水沙运动的影响至关重要。全球气候变化和人类活动影响下，极端冰塞冰坝发生的可能性更大。冰塞冰坝能引起上游河道水位迅速抬高，流凌刮擦割蚀岸滩能导致堤岸崩塌破坏，由此引发的凌汛洪灾严重威胁北方河流冬季输水安全和河流管理。河冰不仅影响泥沙运动和河道变化，还显著影响水体温度和含氧量，例如锚冰和冰盖的形成会压缩水生物生长繁殖空间，进而影响水生态环境［7］。这些河冰过程吸引了众多学者的关注，并在水内冰、岸冰、锚冰、冰盖、冰塞和冰坝等方面取得长足进步［8-11］，但缺少耦合水沙运动和河冰动力过程的研究［3-4］。目前岸滩的崩塌侵蚀研究主要基于无冰条件［12-14］，不能满足北方河流季节性岸滩侵蚀的研究需求。

传统的河流动力学主要研究水流、泥沙与河道边界（包括河床和河岸）三者之间的复杂相互作用，而目前河冰动力学研究主要分析河冰运动与水流间的相互作用。本文首先简要总结了河冰水力学的研究进展，详细综述了有关河冰过程的理论成果和已有研究的不足。这些研究是水冰沙耦合分析的基础，为系统的河冰模拟奠定了坚实的理论框架。在综合各河冰数学模型优点的基础上，本文提出河流水冰沙耦合数学模型［4，15］，既能模拟无冰河流水沙运动和岸滩演变过程，也能计算冬季水冰变化，还能模拟北方河流全季节及河冰生消全过程的水冰沙耦合过程。该研究有助于分析流凌和冰塞冰坝引起的洪水过程及堤岸溃决问题，为北方河流防凌减灾提供技术支撑。

2 河冰水力学理论框架

Beltaos 从河冰热力生长和消融过程、水内冰、锚冰及岸冰的发展、冰塞冰坝形成和释放引起的水位流量波动等方面详细综述了相关理论、数学模型以及原型观测的显著进步［16］，但也强调关于气候变化和人类活动影响下的冰塞冰坝形成机理和预测预报仍有很多不足。Shen［17］详细综述了河冰全过程所涉及的热力、动力和水力过程，指出河冰研究是水力学、冰体力学、热力学和河流动力学等多学科的交叉领域，所包含的物理过程复杂，关于河冰的理论和数学模型在过去30年有了长足进步，能协助解决天然河流和实际工程中涉冰的防洪、发电、航道、生态及环境问题。河冰数学模型的建立与发展为北方河流冬季用水安全和管理提供了有力的技术工具。

受冬季低气温的影响，北方河流水体会释放热量，进而生成过冷却水和水内冰［18］。新生成的水内冰具有较强的吸附性，能黏附于河床、水工结构物及河岸，形成锚冰和岸冰；而大量聚集的水内冰上浮到水面后可以形成表面浮冰。随着气温的进一步降低，河流低流速区域的表面冰盖会出现聚集，直至整个断面被河冰覆盖，形成冬季河流的首封位置［19］。在封河过程中，由于上游径流的变化和气温的波动，在弯道、收缩断面以及阻水建筑物附近能形成冰体堆积。堆积体前大量上游来冰的聚集和堵塞会缩小断面过流面积、显著增加河道阻力，进而引起冰体卡塞点上游水位的壅高，形成封河冰塞［20］。在冬末春初气温回升时，河道冰盖发生热力消融和动力破坏，进而出现开河。在开河过程中，如果水流平稳，河面上的冰盖最终会完全消融，这种平稳的热力开河为“文开河”［21］。如果上游融雪融冰产流大或降雨较多，在过水断面受限的河道会出现大量上游来冰堆积，持续的冰体聚集和上游水位壅高会产生开河冰坝。严重的冰坝能引起洪水漫堤或堤岸溃决，进而诱发凌汛灾害。冰坝引起的河面冰盖破裂和释放常常伴随着急剧水位流量波动，这种水冰动力诱发的开河为“武开河”［11］。河冰的生长、运动和释放会影响水工建筑物结构稳定性，也能刮擦侵蚀河床岸滩，进而影响冲积河流的输沙过程和河道演变。详细的河冰过程和河冰水力学理论框架见图1，主要包括水体失热、河流产冰、封河、开河及河冰影响等五个方面。

图1 河冰水力学理论框架

2.1 水体失热北方河流水体的热交换包括径流和支流的能量汇入、下游出流的能量输出、空气与水体的热交换、水体动能和势能的转换及河床与水体的热交换。冬季水体失热主要是河流与空气间的热交换，占水体热循环的90%以上［22］。持续的水体失热是河冰形成的前提。河床一般在温暖的季节储存能量，在冬季向水体释放热量［23］。河床与水体间的热交换在水体热循环中占比较小，但显著影响河床上锚冰的生长和释放过程［24-25］。水体运动中动能和势能转换的热能一般较小，可以忽略［15］。因此，河流水体热循环主要考虑水面或冰面与空气间的热交换。

针对河流表面热交换，Shen 考虑太阳辐射、水体长波辐射、蒸发热损失、降雨降雪热损失和水体与空气间的热传导等因素，建立了冬季北方河流表面热交换的准确计算方法［26］。除了太阳辐射是水体吸收热量外，其它热传导均为河流向外界释放热量，进而促进水体产冰。该方法能提供详细的水体热交换过程及各影响因子的权重，但需要详细的太阳辐射角、降雪量、河道经纬度、空气透射程度、云层状况、气温、水温和风速风向等各项资料，不便于实际工程应用。为了简化计算，Shen等提出了近似的线性公式计算河流表面的热交换［27-28］

式中：φ为河流表面损失的净热通量，W/m2；φs为太阳辐射的热通量，W/m2；β为与气候有关的经验参数，W/m2；α为水面或冰面与空气间的热传导系数，W/(m2℃)；Ts为河流表面水面或冰面的温度，℃；Ta为空气的温度，℃。针对长期大尺度的河冰数值模拟，式（1）可简化为［27-28］

式中α′为考虑太阳辐射和其它热损失的综合热传导系数，W/(m2℃)。

2.2 产冰当河流水面温度降到0℃后，水体开始产冰。随着水面的持续失热，水体温度能出现一定的过冷却到达-0.01℃的量级，大量的水内冰产生并释放潜热以弥补水体热量损失，接着水体温度回升并维持在结冰温度［18，29-30］。湍流下水体过冷却和水内冰产生计算的理论公式相对完善，具体见文献［18，30-31］。然而，水内冰上浮、悬浮和下沉的运动机理尚不明确，有待研究［32］。

在河流的低流速区如河岸附近易形成静态的岸冰，进而促进断面冰盖的产生。Matousek［33］总结了不同类型流冰的原型观测资料，给出了静态冰盖和运动薄冰形成的水温和流速条件，冰盖的生长过程与湍流强度有关。水流速度过大或者水温不够低时冰盖均不能发展［34］。岸冰的横向发展速率与表面流冰密度成正比，与水流拖曳力和流冰与岸冰摩擦力的大小密切相关［17，35］。

锚冰是一种黏附在河床上的冰，相对浮冰的观测难度更大。锚冰的生长一般在水体达到结冰温度之后，在河面完全冰封之前。当锚冰受到的浮力大于冰体与河床间的吸附力时会从河床释放，进而上浮到水面。锚冰底部因吸收太阳辐射而融化时也会上浮释放。冷冻水槽试验显示锚冰的生长和释放能显著影响河床的综合糙率［36］。锚冰的释放能搬运所吸附的河床泥沙，甚至输运几十公斤的大型卵石或石块，并在融化后释放泥沙［37］。最近Pan等［25］的研究显示锚冰的形成与释放能引起河床高程的变化、断面流量的波动和河道整体糙率的急剧变化，进而造成水深和流量高达30%的增长波动。

2.3 封河当河道出现卡塞点和初始冰盖后，上游的来冰会在卡塞处堆积。如果水流速度较低，断面弗劳德数低于一定临界值时，水面浮冰会以平铺上溯的方式发展（Juxtaposition mode）；如果水流速度较大，浮冰所受的水流拖曳力能拖动冰块下潜或翻转到冰盖下，冰盖会以水力加厚的模式（Hydrau⁃lic thickening mode）向上游发展，又被称为窄河冰塞［16］；当下潜的浮冰过多，冰塞体的内摩擦力小于水流拖曳力和上游压力时，冰塞体会出现坍塌和力学加厚，进而以更厚的冰塞向上游发展（Mechani⁃cal thickening mode），这被称为宽河冰塞，大部分河流的冰塞属于该类型［38］；如果水流速度进一步增大，断面弗劳德数超过上限临界值时，上游来冰在水流拖曳作用下下潜并沿冰塞体底部滑移，最终冲蚀到下游河道，且不能停留堆积在冰塞体下，因此冰盖不会向上游发展［17］。这会在上游河段出现明流区域，并不断产生水内冰。水内冰在下游冰盖下的输移和堆积可进一步形成悬冰坝［29］。

经典的河冰水力学理论采用静力平衡原理分析了冰盖的形成和临界冰厚，将水力学理论应用到冰塞问题分析［39］。Michel 将能量守恒方程应用于冰塞头部计算，建立窄河冰塞下临界流速与冰厚的关系，进一步量化了冰塞的发展过程［40］。Uzuner and Kennedy 采用连续方程和动量方程分析了宽河冰塞下的水流特征，提出非恒定流下宽河冰塞的冰厚计算方法［41］。通过自然河流中大量冰塞的原型观测，Beltaos 与Fan 等［42-43］收集了平衡和非平衡冰塞资料，提出冰塞对水体的阻力与冰厚成正比，并指出冰塞释放引起的岸滩刮擦侵蚀和河床冲刷远大于夏季汛期的河道冲刷［44-45］。沈洪道团队针对冰塞过程中的非恒定水流运动，采用冰水双层流连续介质假定，建立一二维河冰动力学模型，能模拟河冰的产生、输运、堆积、消融及冰塞的发展释放过程，为河冰研究提供了有效的技术手段［17，46-50］。

2.4 开河相比于封河研究，冬末春初的开河研究较少。关于冰盖横缝和纵缝的开裂过程及碎冰盖的卡塞位置和时间仍是河冰研究的技术瓶颈。随着冬末气温的升高和太阳辐射的增强，冰盖消融，冰体强度下降。在上游洪水波的作用下，陡峭河段的冰盖破裂成块，并向下游输运。类似封河时的冰塞形成，开河堆积的冰块也能形成冰塞或冰坝。国内外关于冰塞和冰坝的定义存在较大差异。国外将开河和封河过程中浮冰或流动冰块堆积形成的冰体雍塞均定义为冰塞（Ice jam），将稳定冰盖下流冰花或流冰块的悬浮堆积称为悬冰坝（Hanging ice dam）［16-17，29］。在国内，一般将封河过程中产生冰体堆积和拥堵定义为封河冰塞，将开河时期冰块堆积和拥堵定义为开河冰坝。受流量和来冰量的影响，一般开河冰坝引起的水位和流量波动比封河冰塞大［21］。

Shulyakovskii［49］观测了大量河流开河过程，采用开河和封河时期的水位差作为冰盖破碎与开河的标准，建立冰盖强度与开河水位和封河水位的经验公式。Beltaos 通过大量加拿大河流的原型观测，进一步提出开河与封河期的水位比值与冰厚、冰封宽度、冰体强度和水流拖曳力有关，建立了开河冰盖破裂和移动的概化模型［50］。但这些经验公式具有较大的主观因素，不利于其他河流的应用。国内开河预报主要基于人工智能算法［10］。陈守煜等［51］基于BP神经网络模型，预报了黄河宁蒙河段开河日期。王涛等［52］基于GIS 地理信息系统和神经网络预报模型，建立黄河冰情预报专家系统，能有效预报宁蒙河段开河日期，预报期精度在10 d左右。王军等［53］采用BP神经网络模型模拟了实验条件下弯道冰塞的水位壅高，能较好地预测河冰影响下的水位变化。

受温升融水和降雨的影响，春汛引起的开河冰坝在北方河流最为常见，相应的凌汛灾害也更严重。基于河冰动力学理论的数学模型被成功用于日本渚滑川和加拿大圣约翰河等河流的冰坝过程模拟［54-55］，但开河相关的力学机制尚待完善。开河期碎冰块间歇性的运动和停滞易诱发链式冰坝的形成和释放，伴随的凌汛洪水风险更高。开河冰坝最大的难题是冰坝坍塌释放机理和碎冰块与冰盖间的相互作用机制，包括冰块与冰盖间的相互作用力及水对冰体的作用力等。此外，上游下泄的碎冰块可能撞击挤压下游稳定冰盖，进而导致更多冰盖破裂，相关力学机理还需进一步研究。

2.5 河冰影响冬季河冰运动和堆积可能影响水工建筑物运营及安全［56-58］。封河期生成的大量水内冰能吸附在拦冰栅和取水管口，造成取水口堵塞，进而影响冬季供水和输水安全［25］。流凌运动直接撞击或刮擦桥墩，严重的能引起桥墩混凝土表面脱落，影响桥梁结构的安全性。此外，桥墩结构所在的断面容易造成浮冰堆积，促进冰塞冰坝的形成，引起凌汛洪水［59］。严重的冰坝能在桥梁或其它阻水建筑物迎冰面形成较大推力，造成水工结构的破坏。

Ettema 和Kempema 详细综述了河冰对沙质河床地形、泥沙输运、河岸侵蚀的影响，强调封河期和开河期河冰运动对泥沙运动和河道演变的影响最为剧烈［60］。锚冰的释放能搬运大量泥沙，甚至输运非冰期无法启动的卵石［37］。冰塞冰坝的形成和释放会导致河冰刮擦割蚀河岸，严重的能导致堤岸崩塌破坏，诱发凌汛洪水。另一方面，岸冰的冻结能避免流凌直接刮蚀河岸，进而保护河岸［2，6］。受限于冬季恶劣的天气条件、仪器装备的不足和理论的缺失，目前仍缺少关于河冰对河床冲淤影响的研究，无法定量评估河冰对岸滩的刮擦侵蚀、冰盖下的泥沙输运和冰塞冰坝作用下的河道演变规律。

3 水冰沙耦合模型

3.1 河冰数学模型的基础自从1960年代Pariset和Hausser 将水力学基本方程引入到河冰冰塞分析，众多研究者完善发展了冰塞理论，并促进河冰数学模型的产生和发展［16-17，39，41，61］。基于静力平衡理论，Flato和Gerard开发了一维恒定流下的非平衡冰塞模型（ICEJAM）［62］，Beltaos 建立了考虑冰塞中水体渗流的宽河冰塞模型（RIVJAM）［63］。随后，静力平衡冰塞理论被包含在非恒定一维商业模型中，包括HECRAS和MIKE11-ICE等［64-65］。这些一维模型能预测冰塞冰坝等极限条件下的冰厚和水位上限，为冬季凌汛防治提供了依据，但它们忽略了河冰运动的动力过程，并不能判断冰塞和冰坝形成的条件和位置。

Shen 等将一维非恒定圣维南方程与河冰运动、质量守恒和面密度平衡方程相结合，奠定了河冰动力学模型的基础，能有效模拟河冰的运动以及冰塞的动态发展过程［46］。随后，Lal 和Shen 建立了一维非恒定流河冰数学模型（RICE），考虑水温变化、河冰热力生长过程及河冰运动与水流过程的相互影响，能模拟水内冰、表面浮冰、岸冰和冰塞过程［66］。Shen 等在RICE 的基础上进一步开发了RICEN 模型，扩展了复杂河网模块、锚冰模块和基于输冰率的水内冰输移模块，并成功应用于黄河下游和美国尼亚拉加河上游河道的河冰模拟［27］。随后，该一维河冰模型被升级为CRISSP1D和RICEE，能完整考虑冬季河冰全生命周期的产生、发展和演变过程［48，67-68］。Pan 等基于CRISSP1D 模型，开发了考虑锚冰生长和释放的锚冰洪水波模块，揭示了锚冰引起的河床地形抬高、断面流量变化和河床综合糙率变化，其中锚冰影响下的河床糙率变化是水位和流量波动的主要因素［25］。

为了考虑复杂地形和河岸的影响，Shen 等开发了二维河冰动力学模型（DynaRICE），采用基于欧拉场的有限元法计算二维浅水方程，采用基于拉格朗日式的无网格光滑粒子法计算河冰运动，并应用于密西西比河中游河段冰塞过程模拟，为寒区河流的防凌减灾提供了有效技术支撑［69］。随后，Liu等开发了河冰全生命周期的二维数学模型（Two-dimensional Comprehensive River Ice Simulation System，CRISSP2D），能模拟（1）复杂地形下的急缓流过程；（2）太阳辐射及风场影响下的水温升降规律；（3）水内冰的生成及锚冰、岸冰和浮冰过程；（4）冰盖下的浮冰输移和沉降过程；（5）冰塞冰坝的产生、发展和释放过程；（6）冰盖的热力增长和消融及武开河过程［70］。Knack和Shen耦合二维河冰模块与水沙数值模块，将CRISSP2D 发展为Hydro-Thermal-Ice-Sediment River Dynamics Model，能模拟北方河流河冰影响下的推移质和悬移质泥沙运动及河床变形［4］。此后，Pan 和Shen 将CRISSP2D 发展为RICES2D，考虑河冰影响下的泥沙输运、河床冲淤变化及岸滩侵蚀等，为北方河流水冰沙耦合机理研究奠定了基础［71］，也是本系列文章的前提。

3.2 数学模型的构建基于沈洪道团队长期开发的河冰动力学模型和河冰水力学理论［17］，本文建立了二维水冰沙耦合数学模型，以研究冰盖和流凌作用下的泥沙输移、河床演变及岸滩崩塌侵蚀问题。该模型耦合了水动力过程、热力学过程、河冰运动、泥沙输运、岸滩侵蚀和河床演变过程，模型框架及功能见图2。该模型包括二维水沙数值模块、河冰动力学模块和岸滩侵蚀模块三部分。其中，二维水沙数值模块采用基于非结构网格的有限元法计算非恒定水流运动、推移质输沙率和河道冲淤变形，河冰动力学模块利用无网格的光滑粒子法计算河冰的运动、分布、水力和动力加厚过程，而岸滩侵蚀模块采用双泥沙休止角法计算岸坡崩塌和土体再分布。通过给定耦合时间下的信息传递和反馈，调整变化河冰条件下的水流和泥沙过程，进而实现河流冰期水冰沙耦合过程的模拟。以下分节给出三个模块的控制方程和计算方法。

图2 二维水冰沙耦合数学模型

3.3 二维水沙数值模块二维水沙数值模块包括二维浅水运动、泥沙输运和河床变形三个部分。基于上层浮冰和下层水体的双层流体系，考虑河冰阻力和质量引起的水位流量变化，北方河流河冰影响下的二维浅水方程组为［69］：

式中：x、y 为平面直角坐标，m；t 为时间变量，s；qtx、qty为两个坐标方向上的单宽流量，m2/s；t′i为淹没在水下的冰厚，m；H 为总水深包括河冰的影响，m；Ht为单宽流量对应的有效水深，m；C为河冰的面密度或河冰面积占整个河面的比例；Txx、Tyx、Txy、Tyy分别为不同方向和作用面上的切应力沿水深的积分值，N/m；ρ为水体密度，kg/m3；τsx、τsy分别为两个方向上的表面冰对水体的切应力，Pa；τbx、τby分别为两个方向上的床面切应力，Pa；η为水面高程，m。

浅水方程组采用基于三角形非结构网格的有限元方法求解，能准确描述复杂地形下不规则河道边界。该有限元法基于具有迎风特性的皮德罗夫-盖尔金（Petrov-Galerkin）算法，能有效捕捉干湿边界并适应急缓流变化的条件，甚至能模拟溃坝引起的洪水波传播过程，在多个国家众多河流开展了应用［4，17］。

Knack 和Shen 在总结已有实验和原型资料的基础上，提出适应冰盖条件和流凌影响下的推移质输沙率公式［3］。基于该公式和非平衡非均匀沙方程，考虑河冰影响的泥沙质量守恒方程为［72］：

式中：qbk为k粒径泥沙单宽推移质输沙率，m2/s；k为粒径分组；ubk为k粒径推移质运动速度，m/s；zb为河床高程，m；p′m为河床泥沙的孔隙率；α′bx，α′by为推移质在水流和重力沿岸坡分量综合作用下的方向向量分量。根据不平衡输沙的经验公式，推移质的运动方程为：

式中：Lbk为k粒径泥沙不平衡输移长度，m；q*bk为k粒径泥沙平衡输沙下的单宽推移质输沙率，m2/s。岸坡上的泥沙运动受重力和水流拖曳力的综合影响，推移质泥沙运动方向向量的修正计算公式为［15，72］：

式中：αbx、αby为水流方向向量的分量；β为考虑岸坡和床面泥沙性质的经验性修正系数，可以由实测资料率定或经验公式计算［73-76］；φx、φy分别为床面沿两个方向的倾角。

3.4 河冰动力学模块河冰动力学模块的主要控制方程为河冰质量守恒方程、河冰面密度连续性方程和动量方程［64］。河冰厚度和面密度由河冰连续性方程计算，河冰运动速度由动量方程计算。河冰的动量方程、质量守恒方程、河冰质量及河冰面密度的连续性方程分别为：

式中：M为单位面积冰的质量，kg；为冰速向量，m/s；为冰内部的相互作用力，N；为风作用在冰上的力，N；为水流作用在冰上的力，N；为重力在水平方向的分力，N；Em为单位面积冰的增长速率，由热力生长消融或外部源汇控制，kg/s；ρi为冰的密度，kg/m3；ti为冰的厚度，m；Ra为机械作用引起的冰面积变化，s-1；Ea为冰面积的热力生长消融速率，s-1。为了计算河冰之间的内应力，河冰模块采用黏弹塑性应力应变方程计算冰块的内应力，采用摩尔-库伦屈服准则计算冰塞冰坝中冰块的屈服及冰与河岸的临界摩擦力［77］。

式（9）—式（12）假设河冰颗粒为连续介质，采用基于拉格朗日场的光滑粒子法求解方程组。光滑粒子法不受网格的限制，能准确计算冰盖前沿河冰下潜、冰塞体中冰块的坍塌变形及冰塞的冲蚀和释放等复杂过程，避免了复杂水冰作用下离散格式的不稳定。

3.5 岸滩侵蚀模块自然河流会出现周期性的岸坡冲刷变陡、河岸失稳破坏、坍塌土体落淤再平衡的演变过程。大部分研究采用泥沙休止角判断非黏性沙岸坡的稳定与否，通过原型观测或参数率定确定岸坡稳定的临界坡角，即泥沙休止角［78-80］。Zech 等针对淹没和非淹没岸坡不同含水量的岸坡稳定性，提出双泥沙休止角法：采用稍大的动泥沙休止角判断岸坡的失稳与否，采用静泥沙休止角决定崩岸坡面再平衡的稳定坡面，采用两组不同的动、静泥沙休止角分别计算水面上下的崩岸过程［81］。本文岸滩侵蚀模块采用双泥沙休止角法计算冰盖及流凌刮擦影响下的岸滩侵蚀、崩塌和再平衡过程［71］。计算方法如下。

（1）临界失稳点判断。如图3计算岸坡各个网格点的坡角，判断坡角是否大于相应的动泥沙休止角，并确定临界崩塌破坏点的位置。采用式（13）计算坡面上各点的坡角：

式中αi为横断面i点的坡角。通过对比各点坡角与动泥沙休止角的大小，确定临界失稳点位置，该点以上为不稳定坡面，以下为稳定坡面。图3中NT为横断面网格点总数，Nc为临界失稳点。

图3 横断面坡角及临界失稳点示意

（2）确定崩塌失稳土体质量。以临界失稳点为起点、静泥沙休止角为坡角，确定崩塌后稳定的坡面，计算崩岸土体的质量，具体示意见图4。图中，Ns为失稳坡面上端起点；Ne为结束河床调整的下端网格点；As为崩岸破坏土体面积，m2；Af为崩岸土体重新分布的面积，m2；Δzai为不同网格点崩塌岸滩调整的高程变化，m；Δzbi为不同网格点崩岸土体再分布后的岸滩调整高程变化，m。崩塌失稳土体的面积采用下式计算：

（3）失稳土体堤脚再分布：在保证土体质量守恒的基础上（As=Af），确定落淤土体分布位置，更新岸坡再平衡后各点的地形高程，见图4。通过试算法确定再平衡后堤脚的落淤土体分布。先假定落淤土体以静泥沙休止角为坡角平铺在临界失稳点Nc以下，计算相应的再分布土体面积A'f。如果试算的土体面积等于崩塌的土体面积As，则假定的稳定坡面为最终落淤坡面。如果As≠A'f，则采用式（15）平行调整临界失稳点以下的落淤高程变化。

图4 失稳坡面崩塌及再平衡示意

岸滩侵蚀模块通过坡面各点坡角与动泥沙休止角的对比，判断坡面的稳定与否；利用静泥沙休止角确定崩岸坡面，计算崩塌土体的质量；再通过土体质量守恒计算堤脚落淤土体分布，进而为水沙计算提供更新的河床地形。

3.6 各模块耦合步骤水冰沙耦合模型实施步骤和计算方法是：（1）在满足显格式稳定性的时间步长下运行二维水沙数值模块，直到给定的耦合时间，提供计算区域的水位、流速、水温、水流拖曳力、挟沙力和床面高程等水沙要素，并传递给河冰动力学模块；（2）运行河冰动力学模块获得耦合时间的冰厚、冰速、冰浓度及冰摩擦力等河冰要素；（3）将前两步计算的结果传递给岸滩侵蚀模块，并计算水冰沙影响下的岸滩侵蚀速率、岸坡稳定坡面及崩岸土体再分布位置，直至给定的耦合时间；（4）将第三步计算结果反馈给二维水沙数值模块，校正更新岸坡和地形下的水位、流量和输沙过程；（5）重复以上4 个步骤至计算结束，通过3 个模块的信息传递和反馈实现水冰沙耦合模拟。

4 结论

受季节性的河冰影响，北方河流的输水安全和河道演变需考虑水流、河冰、泥沙和边界的复杂相互作用。本文首先总结了近几十年来河冰研究的主要进展和基本理论，包括冬季水体失热、河流产冰、封河过程、开河机理和河冰影响等五个方面。基于上述理论基础，本文提出了二维河流水冰沙耦合数学模型，创新性地将复杂水流变化、河冰运动、泥沙输运、河床冲淤变形及岸滩崩塌侵蚀过程耦合起来。该模型分为二维水沙数值模块、河冰动力学模块和岸滩侵蚀模块。二维水沙数值模块采用非结构的有限元法计算非恒定水流运动、泥沙输移和河道冲淤变化，能适应复杂河道地形条件。河冰动力学模块采用无网格的光滑粒子法求解，能模拟冰塞冰坝等复杂河冰聚集和释放过程。岸滩侵蚀模块采用双泥沙休止角法计算，能模拟不同含水层岸坡的崩塌和再平衡过程。新模型将河冰理论和水沙理论创新性地融合起来，既满足河冰模拟需求，也适应水沙问题研究，能模拟北方河流全季节和河冰全过程的水冰沙相互作用。二维水冰沙耦合模型的初步建立为北方河流防凌减灾和河流管理提供了有力的技术支撑。

二维水冰沙耦合模型主要针对沙质河道和岸滩而开发，尚不包括黏性土质岸坡和河床冲淤计算模块，也不包括岸滩植被分析功能，不适于黏性岸坡稳定性分析和密集植被影响的河道模拟。针对非黏性土和黏性土混合的复杂河道条件，还需进一步拓展该模型的水沙数值模块和岸滩侵蚀模块。

关于河流水冰沙耦合模型研究的系列文章分为两篇，本文主要介绍原理和方法，模型具体验证与应用将在下篇给出。