崔弘毅编 译

（国家电力监管委员会大坝安全监察中心，浙江杭州310014）

1 专题简介

4家研究机构参与了IMPACT项目泥沙运动专题的研究，分别为：天主教鲁汶大学（UCL，比利时）、特伦托大学（UdT，意大利）、里斯本高级技术学院（IST，葡萄牙）和法国农业暨环境工程研究所（Cemagref，法国里昂）。

本文主要阐述研究团队在该项研究涉及到的不同区域中的主要发现：近坝址和远坝址溃坝洪水、对溃坝洪水的试验模拟和数字模拟以及通过与试验室数据和现场数据的对比进行基准试验来校验模型。

2 现场数据和试验室数据

2.1 目标和方法

2.1.1 试验室数据

试验室试验可使研究人员重点关注某些具体的过程，可以对其进行更详细的分析研究。通过充分的测量，可得到精确的数据，数据用途有二：一是可以帮助更好地理解正在发生的现象的物理原理，二是为数字模型提供校验基准。在天主教鲁汶大学（溃坝洪水）和特伦托大学（均一泥石流）的试验室中进行了一系列试验，考虑了两类性质：

（1）在近坝场，伴随着溃坝波的演进，发生了快速而强烈的冲蚀。洪流显示了强烈的自由表面的特性：波破碎发生在中间位置（在坝址附近），接近垂直的水和残渣墙在波浪前沿掀起了泥沙底床，造成瞬时强烈的泥石流（见图1）。但是，在溃坝波前沿，泥石流却意外地和均匀波相似。因此，该项研究的第一项研究内容是在更深入地调查研究溃坝流条件下的性态之前，研究均一条件下的泥石流的特性。

（2）在远坝场，泥沙运动依然很剧烈，但泥沙的动态作用减弱。因泥沙分散、两岸冲蚀和残渣沉积，山谷地貌发生了戏剧性的变化，见图2。第三项研究内容就是研究远坝场的性态。

图1 近坝场的泥石流Fig.1 Near-field geomorphic flow(UCL)

图2 断续的坝段溃决导致的两岸冲蚀Fig.2 Bank erosion resulting from intermittent block failure

2.1.2 现场数据

泥沙运动主题研究所利用的现场案例是Lake Ha！Ha！溃决事件。Lake Ha！Ha！位于加拿大魁北克沙格奈河的一条支流上，于1996年溃决（Brooks&Lawrence，1999年）。在该溃决事件中，发生了大量的泥沙运动，山谷地形地貌发生了巨大的改变：天然河道变位偏移，两岸冲蚀引起河道中等或大幅变宽，受不可冲蚀的基岩的影响，洪流进一步偏离发生变道（洪流被撕裂）。该项研究的目的是收集与高度瞬变流相符，包括强烈泥沙转移的现场数据。对相关文献进行核查，以便发现最合适的案例，并有足够可用的数据。

2.2 收集到的数据

2.2.1 试验室数据

可用的数据系列有：

（1）均一泥石流，用以调查研究作用力和速度分布（UdT）。

①均一材料（PVC颗粒）；

②级配材料（PVC颗粒和沙）。

（2）溃坝洪水，用以调查研究近坝场的影响：冲刷和泥石流流阻的形成（UCL）。

①大坝上下游河床高程一样（Spinewine&Zech，2002b）；

②河床高程有一个抬升，水库内的高程更高（Spinewine&Zech，2003）。

（3）梯形河谷中的溃坝洪流，用以调查研究远坝场的影响。两岸冲蚀和河道拓宽（UCL）（Le Grelle等，2003,2004）。

①均一材料（沙）；

②级配材料（沙和粗砾石）。

2.2.2 现场数据

关于Lake Ha!Ha!溃坝事件，收集到了洪水前后的大量数据。数据处理工作由天主教鲁汶大学（UCL）、台湾大学、魁北克大学和加拿大地质调查局联合进行，其中，加拿大地质调查局是数据所有者。整套数据包含完整的洪水前后的数字地面模型（30 km范围）以及重新制作的堤坝溃决的出流过程线。

2.3 价值（潜在应用）

2.3.1 试验室数据

试验取得了高质量的数据系列，用于校验数字模型。

2.3.2 现场数据

现场数据可帮助识别和理解由溃坝洪水引起的地形地貌变化的关键特征，以便涵括未来数字模型的必要方面。另外，通过现有数字模型得出的结果和观察结果之间的比较，可评价现有数字模型的性态。

3 模拟极端洪水条件下泥沙运动的新方法

3.1 目标和方法

研究的目的是对调查过程进行数学描述。在将现有的洪流描述延伸去解释试验认定的过程之前，需对文献进行复核。

3.2 分析和发现

3.2.1 均一泥石流

由于快速颗粒流和气体之间存在某些物理相似性，将动力理论用于颗粒材料需耗费大量工作。所有模型都假设微粒间的相互作用是由瞬时碰撞引起的，这就意味着只需考虑二元或双微粒碰撞。

Jenkins和Hanes（1998年）曾将动力理论用于薄层水流，其颗粒由碰撞互相作用支持，而不是由紊流的速度波动支持。颗粒压力的本构关系被视为与Chapman和Cowling于1970年提出的稠密分子气体的准弹性近似，其描述了颗粒间碰撞率集中度的变化。

假设颗粒浮重整个被颗粒碰撞接触所支撑，则从试验中可能推导出颗粒压力σs和剪切应力τs。该项目框架中关于明确本构关系的最主要进步为解释了附加质量效应，即通过公式（1）替换泥沙密度ρs：

式中，Cs为颗粒集中度。

3.2.2 近坝场溃坝洪水

3.2.2.1 2D-V水平集模型

考虑到溃坝洪水第一阶段速度的垂直分量不可忽略，最初是想开发一个2D-V模型，可以表示出大坝溃决后最初的瞬间垂直中正面发生的情况。最合适的模型貌似是水平集方法。该方法基于的假设为：洪流被细分为性态近似相同的层次，由尖锐的界面分开。多种介质（空气、水、泥沙）中行进的界面对应更高维函数Φ的零点水平集，被定义为到界面的符号距离（Sethian，1999年）。就涡度和流函数Ψ而言，可建立Navier-Stokes方程，从方程可得到速度场。根据速度场，则可得水平集方程，得到有用的定性结果，见图3。

图3 水平集方法原理Fig.3 Level-set method principle

3.2.2.2 双层浅水1D模型

最初的进展由Capart于2000年提出。洪流由三个层次代表：（1）上层水，由清水组成，深为hw；（2）运动的泥沙层，厚度为hs；（3）固定的基岩层，基岩高程zb，作为上限。在最初的模型中（Capart，2000年），假设泥沙的集中度不变（Cs=Cb），且水和泥沙（hs）混合体的上部都与清水层一样，以同样均一的速度运动（us=uw）。根据这些假设，剪应力应该在垂线上是连续的。模型得出了解析解（Fraccarollo和Capart，2002年），尽管很灵巧，但不能用于实际的几何构型。

关于模型的最主要进展之一（Spinewine，2003；Spinewine和Zech，2002a）是赋予集中度和三个层次间速度以新的自由度（Cs≠ Cb，us≠ uw），见图4。

该描述中得出的方程由二阶Godunov有限容积法解答，其中流量用LHLL Riemann解来计算（Fraccarollo等，2003年）。

图4 近坝场洪水流数学描述的假设Fig.4Assumptionformathematicaldescriptionofnear-fieldflow

3.2.3 远坝场溃坝洪水

3.2.3.1 二维模型

首先，开发了一个2D可拓模型来表示近坝场，包括岸坡冲蚀机理。这里主要总结其方法，详细信息参见Spinewine等（2002年）和Capart&Young（2002年）。关键点在于认定分离的水流和液状泥浆层各自独立流动，控制方程就完全能处理岸坡料坍塌陷入水流的滑塌事件。一旦溃决发生，则溃决后的洪水就像其它水和泥沙的运动模式。

因此，需要液化标准来明确什么时候及什么地点，岸坡从固态转换成液态介质。对此，假定以下基本机理：当局部边坡超过临界角度αc时，则发生坝段溃决；延伸的溃决表面被确定为锥形，中心在溃决处，以剩余角αr＜αc向外倾斜。最后，假设椎体以上的泥沙材料在溃决时立即液化。为解释观察到的淹没区域和未淹没区域的区别，如图5所示，定义了4个不同的安息角：αc,subm和φr,subm用于描述淹没区域，αc,em和αr,em用于描述未淹没区域。

图5 2D地压溃决算子稳定示意图Fig.5 Stability diagram for the 2D geostatic failure operator

3.2.3.2 全部岸坡溃决的一维模型

第二个选来耦合以上岸坡冲蚀机理的模型是一个一维方案，它包含了一个流体动力有限容积法和一个单独的泥沙运动路线。有限容积法，其开发目的是为应对复杂地形（Soares-Frazão和Zech，2002年），它解出了流体动力浅水方程，因在一个计算时间步长上的纵向泥沙运动（推移质）导致断面地形发生的部分改变可以从泥沙相的Exner连续性方程得出。除单元格上下游面的泥沙流量外，岸坡溃决导致的侧向泥沙流入量设为体积Vs，其在计算时间步长末端的断面上将重新分布。

水位升高Δh淹没岸坡，引起如图6所示的棱柱形部分材料失稳，最终引发溃决。试验中，岸坡原始角α小于水面上稳定角αs,em，但大于水面下的稳定角αs,subm。因此，当水位上升时，岸坡立即变得不稳定，当达到溃决角度时，则发生溃决。对应于淹没工况和出露工况，溃决角度分别为αf,subm和αf,em（在实际情况中，溃决角度αf比稳定角度αs小一点）。断面上的冲蚀体积Vs也重新分布。

图6 因淹没岸坡引起的岸坡溃决Fig.6 Bank failure triggered by the submergence of the bank

如图7所示，冲蚀材料沉积在河道中，淹没部分，对应水下的安息角，泥沙堆积角度为αr,subm，出露部分则稳定在角度αr,em（堆积过程发生后的水上安息角）。所有这些安息角都由试验中所用材料决定，测量都在静态和动态试验中完成。

图7 岸坡冲蚀材料堆积示意图Fig.7 Deposition of the material eroded from the banks

最终，数字一维模型关键在于以去耦合的方式解决过程中的三个不同关键步骤：（1）水流的流体动力路线；（2）纵向泥沙运动和由此导致的冲蚀和沉积；（3）岸坡溃决和由此导致的断面形态的变化。

3.2.3.3 局部岸坡溃决的一维模型

以上方案在理想工况中适应性良好，其中，理想工况是指断面明确，例如是长方形或梯形，在其棱角方面，也只有有限的几个凸峰。但对天然河道来说，其断面复杂得多，不能用这么简单的方式来描述。在描述天然河道时，推荐使用Schmautz和Aufleger于2002年提出的方法（本文不介绍）。

断面剖面被切分成小块，从山谷边开始，对每一小块的稳定性进行复核（如图8所示）。若岸坡断面AB局部比临界值（稳定角αs）更倾斜，则岸坡部分会发生转动，直至达到A’B’位置，对应安息角αr（淹没或出露）。结果，这一新位置会加剧周边区域的稳定性（例如BC已移至B’C位置）。需多次观察整个剖面，直到所有区域都稳定。

图8 局部岸坡溃决模型原理Fig.8 Principle of the local bank-failure model

对于纵向的泥沙运动，无论对全部或是局部岸坡溃决模型，以下规则都适用：（1）若发生冲蚀，根据 Meyer-Peter.Müller公式，假设运动与(τb-τb,c)3/2的局部值成比例，其中τb和τb,c分别为河床的实际剪切应力和临界剪切应力；（2）若发生沉积，则假设泥沙沿河床均匀沉积（不是水平地），这在后面被定义为断面因素，其坡度小于水下的安息角。

3.3 主要结论和对行业的建议

3.3.1 2D-V模型

事实很快表明，在现实案例中使用这么复杂的模型几乎是不可能的。提出整个现象沿特定垂直面的表达的约束条件是不现实的，因为对大多数真实河谷来说，其在大坝附近变得非常窄，而在横河向，水深变化极大。因此，2D-V方法只能作为通往完全3D方法途中的一个有趣步骤，而完全3D方法是建模者的远期目标。

3.3.2 1D和2D-H模型

描述溃坝洪水条件下泥沙运动的方程组还未完全建立起来。因此，现阶段不可能开发出用于商业用途的工具包来解决随溃坝事件而产生的与严重的泥沙转移有关的各种问题。现有模型，包括其延伸的各种可能（如之前提到过的2D-H模型），在实际应用中都太慢太缺乏效率。

3.4 主要结论和对未来科学研究的建议

3.4.1 2D-V模型

2D-V方法只能作为通往完全3D方法途中的一个有趣步骤，而完全3D方法才是建模者以后的重要攻克目标。

3.4.2 1D和2D-H模型

所有比较的模型中，至少在平坦河床基准试验中，都存在一个普遍的缺点，那就是它们推进前沿太快，这是因为没考虑垂直速度分量，因此造成了错失泥沙流动的初期阶段的事实。另一明确的结论是，要重现饱和残积物前沿的冲蚀性态很困难，特别是之后发生部分再堆积的冲蚀。

但是，与数年前可获得的结果比较，这样的建模过程是令人惊叹的。取得该进展的部分原因是基于新的数字成像测量技术，它使对液相和泥沙运动层的速度场进行实时观测成为可能。

二层一维方法及其在2D-H模型中的延伸看起来都是有趣且有前途的方法。未来的努力方向是对岸坡溃决算子的描述，包括岸坡溃决的引发，加入对土力学的考虑就可以实现。从数字模拟的角度看，还应努力使计算在合理的时间内完成。

4 现有模型的基准试验

4.1 目标和方法

共进行了3个基准阶段的测试，不同的建模者可以在IMPACT项目中进行盲模试验测试其模型。通过比较盲模结果及分析不同数学描述后的模拟结果，可以知晓每种建模方法的优点和局限性。对每条基准，要求建模者提供其数值模型的描述，从而可以深入分析结果。

4.2 收集到的结果

4.2.1 流过原本平坦河床的溃坝洪水

Spinewine和Zech（2002b）给出了基准描述，有4家研究机构得出了结果，即参与IMPACT项目泥沙运动专题的小组：法国农业与环境工程研究所（Cemagref，法国里昂）、特伦托大学（UdT，意大利）、里斯本高级技术学院（IST，葡萄牙）和天主教鲁汶大学（UCL，比利时）。图9显示了典型结果，其中可见计算所得的河床高程、运动的泥沙层高程和洪水高程与试验数据的对比。

4.2.2 流过原本阶梯状河床的溃坝洪水

图9 溃坝波流过原始平坦易蚀性河床的基准试验结果和数值结果Fig.9 Experimental and numerical results from the benchmark on dam-break wave over an initially flat erodible bed

Spinewine和Zech（2003）给出了基准描述，有4家研究机构得出了结果，即参与IMPACT项目泥沙运动专题的小组：法国农业与环境工程研究所（Cemagref，法国里昂）、特伦托大学（UdT，意大利）、里斯本高级技术学院（IST，葡萄牙）和天主教鲁汶大学（UCL，比利时）。阶梯状河床基准下的各模型结果与试验数据的比较见图10，图中显示了在给定时间的不同高程。

图10 溃坝波流过原始阶梯状河床的基准试验结果和数值结果Fig.10 Experimental and numerical results from the benchmark on dam-break wave over an initially stepped bed

4.2.3 原本棱柱形河谷中由溃坝洪水引发的岸坡冲蚀

Grelle等（2003年）给出了基准描述图11为其基准，4家研究机构得出了结果（图12），即参与IMPACT项目泥沙运动专题的小组：法国农业与环境工程研究所（Cemagref，法国里昂）、特伦托大学（UdT，意大利）、里斯本高级技术学院（IST，葡萄牙）和天主教鲁汶大学（UCL，比利时）。

图11 试验测量Fig.11 Experimental measurement

4.3 分析和结论

4.3.1 流过原本平坦河床和原本阶梯状河床的溃坝洪水

关于前沿波速，特伦托大学（UdT）的结果利用了校准过程，包括用这些波速数据作为校准参数。与之形成对比的是，运动泥沙层估值偏低，原因是假定运动泥沙层的集中度与河床材料一样，但天主教鲁汶大学（UCL）和里斯本高级技术学院（IST）的模型情况并非如此，事实是，运动泥沙层的集中度在降低，以使颗粒可以运动。由前沿移动造成的冲蚀仅在天主教鲁汶大学（UCL）和法国农业与环境工程研究所（Cemagref）模型中出现，尽管法国农业与环境工程研究所的简单模型不能提供任何关于运动泥沙层的结果，但仍得出了激波过后的关于水体表面的有用预估值。在这一点上，对冲蚀和沉积的不对称处理可解释为天主教鲁汶大学模型的成功之处。

4.3.2 原本棱柱形河谷中由溃坝洪水引发的岸坡冲蚀

法国农业与环境工程研究所（Cemagref）模型中仅有河床运动，不考虑断面的典型增宽，而其显然是河岸稳定标准的需要。在特伦托大学（UdT）模型中，河岸没有重塑形，而是所有河岸材料都被转移到了底部，明显导致河床高程的估值偏高。在描述岸坡失事时，里斯本高级技术学院（IST）的模型奇怪地失效了，而河床加深会引发该机理的发生。天主教鲁汶大学（UCL）所使用的模型中，对因岸坡垮塌而导致的材料沉积，明显利用了安息角的定义。

必须注意的是，试验中，原本岸坡角度比临界角大，强调了该现象。对于更平缓的岸坡，形态效应不那么重要。

5 评估模型不确定性

5.1 目标和方法

大坝或堤坝失事产生的洪水可引起多种形式的、严重的泥沙运动。某些案例中，洪水带走的材料体积可与溃坝释放的洪水达到同一数量级（达几百万立方米），因此，泥沙运动的风险是巨大的。

现有模型还不能对大多数不确定性相关因子的识别做出详细而敏感的分析。在模拟洪水事件中，重点在强调对泥沙运动核算的重要性上。

5.2 分析和发现

强烈的冲蚀和泥沙转移使河谷地形发生了巨大而快速的变化。反过来，地形变化也对波浪性态有强烈的影响，因此，到达时间和最大水位高程也受到影响，而这是风险评估和组织预警的两个关键参数。这就意味着影响泥沙运动预测的不确定性因素可能最终影响整个预测过程。为证实这些影响，对波浪在固定、光滑的河床和波浪在有移动的泥沙上的传播性态进行比较（见图13）。

图12 10 s后的断面图Fig.12 Cross section after 10 s

图13 固定河床（点线）和移动河床（实线）上的溃坝波在t=1.5 s时的比较Fig.13 Comparison between dam-break wave on fixed(dotted lines)and mobile bed(solid lines)at t=1.5 s

从图13中可见，泥沙的移动转移了部分有效势能，因此波的前进速度大大放缓，这对预警和下游人口的应急规划来说是一利好。但波前沿后面的水深大幅增加，至少在近坝场有大幅增加，增大了危险区域和溃坝受灾区人民的风险。

5.3 主要结论和对行业的建议

此类现象中的大多数过程是不确定的。获取建模所需的数据通常也非常困难。河床底部的组成材料不均一，厚度也不明确，材料会随时间严重离析，特别是细粒材料。大坝下游河谷的河床材料也不均一，其中包含的土和岩石以不可预知的组合方式排列着。对河床进行测量非常困难、成本太高而且乏味。

这就意味着，现在还没有标准的做法可推荐用于考虑泥沙对溃坝波的影响。同样意味着终端用户必须清楚知道其模型在这方面存在着局限性。

5.4 主要结论和对以后科研进程的建议

模型是在对溃坝的理想化状态基础上进行模拟的。问题被表示成为一个竖直平面，且大坝被假设为溃决后立即消失，没考虑横向效应。在近坝场泥沙运动中，仅考虑了河谷河床材料，忽略了溃坝本身释放出来的材料。模型发展现阶段，移动河床建模还没与溃坝建模耦合起来。模型对理想状况来说非常有应用前景，但要用于真实案例中的情况，还有很大差距。

对远坝场，关键是要重现河谷的演变过程，即因上游固体材料输移及坝堤溃决引起的持续冲蚀和泥沙沉积过程。对此，可以模拟部分地形变化，尤其是局部的地形变化，但仅限于几公里范围内。因存在许多随机现象，因此溃决发生后的一系列事件变得不可预测，从而形成类似“不确定树”现象，很难掌控。

6 现场特定案例研究

6.1 概览

经过数据整理，真实的案例研究包括：（1）数据解释；（2）将合理的假设应用到数字模型中：确定网格，一维还是二维？边界条件情况如何，有时候并不能获取所有边界条件的物理数据；（3）解决问题（运行计算）；（4）对结果进行评判性分析，评估价值和结果的“真实性”。

6.2 Lake Ha！Ha！案例的概况

1996年7月，Lake Ha！Ha！水库副坝溃决，产生了强烈的溃坝波，导致下游30 km长河道的地形发生了重大改变，在Ha！Ha！湾汇入Saguenay河（见图14）。

图14 Lake Ha！Ha！示意图（Brooks，2003年）Fig.14 Lake Ha!Ha!(Brooks,2003)

事实上，可以观察到由溃坝灾难引起的巨大地形变化的所有典型特征：大范围的沉积，以致于河流河道发生变化（见图15（a）），河道大幅加宽（见图15（b）），有时因河床基岩底坎和河岸的出现而发生阻塞（见图15（c）），改变了河床的地形地貌，改变了河流的路线等。

图15 Lake Ha！Ha！溃坝洪水波洗劫后的典型地形变化（Brooks，2003年）Fig.15 Typical morphological evolutions after the Lake Ha!Ha!dam-break wave(Brooks 2003)

加拿大地质调查局、魁北克大学、台湾大学和天主教鲁汶大学（Capart等，2003年）做了巨大的努力来对获得的数据进行判读，以获得有用的数据系列，这或许是现实案例中可用来校验模型的最好的数据系列之一。

图16为发生大规模地形变化的周边区域的河床剖面变化图。原河床剖面主要为不可冲蚀的岩石区域，对应图16中绿色的斜线段。因在河岸线的压低处发生了漫顶，河流河道发生改变。河床岩石区域被绕开，引起严重的溯源冲刷。

虽然法国农业与环境工程研究所（Cemagref）的一维模型没使用关于泥沙运动的复杂的描述（普通泥沙转移的Exner方程），但结果得到的河床剖面变化方向正确，尽管存在数值不稳定性。其中存在一些差距的原因与岩石位置有关，且与这一实际情况有关：计算在2天后就停止了，而此时冲蚀过程还没有完成。

图17为与图16相对应的相同河段的水位剖面图。从图中可见河床移动和地形变化的作用。水位高程数据来自两种不同的计算，即法国农业与环境工程研究所（Cemagref）的模型，计算工况为移动河床和与之形成对比的固定河床法（El Kadi和Paquier,2004年）。在有些地方，因河床的移动性，导致水位上升或下降达5 m（例如图17中21～22 km处）。

图16 Ha！Ha！河溃坝前后的河床剖面，对比数字模型Fig.16 Bed profile of the Ha!Ha!River before and after the dam break.Comparison with numerical models

图17 Ha！Ha！河接近洪峰时（7月20日19∶30）的水位高程剖面图（法国农业与环境工程研究所（Cemagref）的模型：固定河床和移动河床间的对比（2种方法））Fig.17 Water profile of the Ha!Ha!River near the flood peak(July 20,19:30).Cemagref numerical model:comparison between fixed and mobile bed(2 approaches)

台湾大学使用的扩散移流模型是二维的、各向异性的，但其基本假设却很简单，该模型对Lake Ha！Ha！这一典型案例的应用来说非常有效。根据当地具体坡降，用沿坡的扩散来表示地形变化。不要求明确的流体力学计算：在凹陷处，地下水面是水平的，而其它地方的深度为0。

此方法得出的结果令人印象深刻，至少在有些河段给人留下了深刻印象，比如上游河段（见图18）。二维模型貌似可以捕捉到刚好溃决的坝体下游处的冲蚀区域和主要沉积区域。一维模型也可以预测某些沉积特征，并且，在一定程度上，还能得到冲刷区域的时间精确的结果。对有些河段的加宽也有定性的表示。

在Lake Ha！Ha！这样的复杂案例中，仅有最简单的方法能成功得出某些结果，复杂的模型也无法应付海量的必需的数据，这是非常有意思的现象。

[1]Final Technical Report of Investigation of Extreme Flood Processes And Uncertainty[R].2005.

[2]H.CAPART,B.SPINEWINE,D.L.YOUNG and et al.The 1996 Lake Ha！Ha！breakout flood,Quebec:test data for geomorphic flood routing methods[R].Revision Submitted to the Journal of Hydraulic Research.