土石坝溃决往往造成破坏性的洪水，尤其是滑坡、泥石流、终碛堤，以及火山碎屑等地质过程堵塞河道形成的天然土石坝及其溃决过程，由于事先未有任何征兆和防护措施，其造成的危害往往更为巨大。为了能够准确的对坝体溃决所造成危害进行评估、制定合理的安全预警方案以及实施有效减灾工程措施，模拟坝体的溃决过程以及溃口的洪水过程显得十分重要。
1溃口几何形态溃口的几何形状对溃口水力过程影响很大。溃口的形状通常概化为矩形、三角形或梯形，或将断面概化为抛物线形。根据对大量的人工土石坝溃决事件以及汶川地震灾区堰塞湖溃决的调查发现，溃决后的坝体溃口多呈梯形或两侧为均匀斜坡底部为抛物线的形态特征。在本计算模型中，将溃口的几何形态概化为梯形，溃口形状由:溃口底宽b、溃口深度z和溃口边坡参数s3个参数确定。溃口通道坡度较小，在计算时可假定为水平，u和d分别为坝体迎水面和背水面边坡。
2溃口堰流方程
溃口假定为梯形，溃口流态可以认为是急变流，通过溃口或坝顶水流的水力特征用宽顶堰流描述。溃口出流过程是溃口几何形状和溃口水头的函数，可以用下式表示
Qb=1.71b(H-Z)1.5+1.2s(H-Z)2.5(2)
=1H2-Z23(H-Z)1-2.783(H2-ZH-Z-0.67)3H2-Z23{(H-Z)(3)
式中H为堰塞湖水深，Z为溃口相对于坝顶的高程，s为溃口边坡，取决于坝体组成物质的内摩擦角()，可取s=tan(45/2)，为淹没系数，H2为坝下游正常流深，可以近似用曼宁公式表示
H2=(Qbnd1.49J0.5B0)(4)
式中n2为曼宁糙率系数，J为坝体下游河道坡度，Qb为前一个时段溃口出流量。
3库区水量平衡
堰塞湖或水库库容可以用水位-库容关系表示
W=aHn(5)
式中W为库容，H为水深，n为库容指数。对于特定的坝体，其水位库容曲线可以由地形图绘制。上式适用于堰塞湖表面为水平时，一旦溃决开始，湖水就会向溃口集中并沿溃口下泄，这时由于库区水面坡度较小，仍可认为上式适用。这时堰塞湖水量平衡可以表示为
dWdt=Qi-Qb-Qs-Qp(6)
式中Qi为入库流量，Qb通过溃口的流量，Qs为通过溢洪道的流量，Qp为管涌出流量。由于溃坝时间尺度较短，因为堰塞库区的蒸发作用可忽略不计;对于天然土石坝漫顶溃决，若无溢洪道且管涌量较小，可忽略通过溢洪道和管涌的流量。则通过式(5)和式(6)可得
anHn-1dHdt=I-Q(7)
式中I为入库总流量，Q为出库总流量。如果在坝体溃决的时间段内，入库流量可以认为是常数，则上式中只包含两个未知量，堰塞湖水面高程H和溃口的出流量Qb，通过水深和时间的对应关系，可以求解出库流量。
4溃口泥沙输移
当堰塞湖水位到达坝顶时，水流漫过坝顶从而对坝顶形成冲刷过程;被侵蚀的泥沙由水流带走，从而使溃口不断扩大。溃口的大小和溃口流量决定于溃口水流的冲刷能力，而水流的冲刷能力又与溃口大小和流量有关，二者是相互联系的。溃口泥沙侵蚀和输运能力的大小与溃口的形状、筑坝物质的性质、坝址下游的水位等水力因素有关。
目前溃口泥沙的侵蚀和输运机理，特别是高速水流下的泥沙侵蚀和输运机理仍不清楚，许多模型使用泥沙输移方程表示溃口的泥沙输移和溃口下切拓宽过程，例如广泛应用的Meyer-Peter-Muller公式(如Breach模型)和Einstein-Brown推移质公式(如beed模型)。
但是由于这些公式的基本假定为缓变流和水深远大于泥沙颗粒粒径的条件。而天然土石坝其组成物质小到粘粒，大至几十米的砾石，级配很宽，且溃口流态为急变流，因而这些公式是否适用于天然堰塞坝的溃决仍需继续研究。
5计算方法及程序化
模型可以用来模拟坝体溃口侵蚀过程，估计坝体溃决参数，包括溃口流量过程线，溃决持续时间以及溃口在溃决过程中和溃决结束后的几何形状。模型计算假定溃口形状为梯形，需要输入溃口边坡，坝体物质组成和几何特征，堰塞湖的水深-库容关系以及堰塞湖入流量等参数。
计算过程中主要变量为每一时间段的堰塞湖水位H和溃口底部高程Z(或溃口深度z)。设计算时间步长为t，则式(7)可以写差分形式
Hn(i+1)-Hn(i)t=1a[I-12(Qi+1+Qi)](20)
6天然土石坝溃决实例分析
天然堰塞坝的溃决过程鲜有较为详细的记录。汶川地震中形成的唐家山堰塞湖的安全泄流过程，为探讨模型模拟天然土石坝溃决过程的可行性，以及模型参数的选择等问题，提供了可靠的数据支持。
6.1唐家山堰塞坝概况
唐家山堰塞坝位于北川县城上游直线距离3.2km，苦竹坝水电站上游2km的通口河上，集水面积3550km2，坝址区通口河谷为不对称的V型谷，右岸较陡，坡度约为60左岸坡度较缓，坡度约为30。岸坡为残坡积碎石土层，最大厚度15m。
唐家山滑坡为顺层岩质滑坡，主要岩性为深变质砂板岩，滑坡体长600m、宽200m、厚度80～90m，滑坡堆积体积2.04107m3。堰塞体碎石土由粉质壤土、岩屑和块石组成，其中粉质壤土占50%～60%，岩屑30%～35%(粒径5～20cm)，块石占10%～15%;坝顶左侧最大高程785m，右侧为755m，中部最低为752.2m，坝趾高程为669.5m，最大坝高122.2m，有效坝高80.3m;坝体迎水面坡度为20，背水坡平均坡度为32。
6.2实测溃决洪水过程
2008-06-07T08:00，坝前水位涨至740.43m，水库开始通过人工溢洪道溢流。06-09T14:00，由于水流的不断冲刷，溢洪道开始初现下切和拓宽过程。06-10T01:00，堰塞湖达到最高水位743.1m，T07:00通过溢洪的流量已达497m3/s，堰塞湖开始泄洪，水流的强烈冲刷使得溢洪道迅速的拓宽和加深，通过流量迅速增大;T08:00泄洪流量已增至940m3/s，坝前水位742.8m;T12:30出现洪峰流量，坝前水位735.8m，流量为6500m3/s;随后，溢洪道的流量迅速的减小，T20:00，唐家山堰塞湖坝前水位回落至719.48m，蓄水量为1.12108m3，泄流槽泄流流量为290m3/s。至06-11T07:00，坝前水位715.2m，流量为71.8m3/s，基本与通口河5月中旬日平均流量相当。唐家山堰塞湖排险泄流过程基本结束，工程排险泄流取得成功，泄流过程线。唐家山堰塞坝泄流过程中，堰塞湖共排出水量1.7108m3，水位降低27.9m。
7参数敏感性分析
模型中需要输入参数包括溃口边坡，坝体物质组成和几何特征，堰塞湖的水深-库容关系以及堰塞湖入流量等4个方面的，溃口边坡可以通过坝体物质天然内摩擦角确定，坝体几何特征和堰塞湖的水位库容关系可以通过测量加以确定，这3个参数对于特定坝体而言都为定值。但是代表坝体物质组成的物料系数却为经验参数，其取值取决于对坝体特征的判断，文中在其他参数相同的条件下，通过不同的物料系数，对唐家山堰塞湖的溃决过程进行了演算。
从表5中可以看出，模型计算中物料系数对坝体溃决过程的相关影响较大，尤其是对峰值流量的影响最大，其计算结果相差将尽7倍，当和时计算相差也近3.5倍，因此在模型计算时，物料系数的选择应非常谨慎。从唐家山堰塞坝的物质组成来看，粉质壤土和较小岩屑占到80%以上，块石含量较少。模型计算过程中选择土类3(含土料较多的土石坝)，从计算结果与实测值得吻合程度来看，选择比较合理。
8结语
通过溃口水力过程，泥沙输移过程，坝体物质组成和几何特征，建立了土石坝漫顶溃决过程的计算模型。模型中，使用堰流方程来模拟溃口的出流过程，使用经验公式确定溃口最终的几何形状和溃口处得泥沙冲蚀量，并以此确定溃口在某一个瞬时的深度。通过堰塞湖的水量平衡确定其水深随时间的变化过程，从而确定溃口水头和溃口流量过程。由于天然堰塞坝体在几何形态和物质组成上的复杂性，对于溃口最终几何特征和溃口的泥沙输移过程，目前仍只能通过经验公式加以确定。
模型需要输入参数包括溃口边坡，坝体物质组成和几何特征，堰塞湖的水深-库容关系以及堰塞湖入流量等4个方面有，可以计算口流量过程线，溃决持续时间以及溃口在溃决过程中和溃决结束后的几何形状。
汶川地震中形成的众多堰塞湖中，唐家山堰塞坝体体积和蓄水量最大，其安全泄流过程引起了全世界的关注，也为天然堰塞坝溃决过程的模拟提供了有力的数据支持。但是从堰塞坝溃决风险评估角度来看，事先制定的1/3瞬溃方案显然有极大的不合理性，其计算所得峰值流量也与实测值有很大出入，其主要原因在于:1.天然土石坝的溃决过程不可能为瞬时溃决;2.其最终溃决的深度难以确定。模型中，考虑了堰塞湖的水量平衡，通过经验模型模拟溃口的冲蚀过程，通过宽顶堰流方程模拟溃口流量过程，与天然堰塞坝溃决物理过程比较吻合。通过唐家山堰塞坝溃决过程反演，可以看出模型计算结果与实测值比较吻合;但是由于在模型计算中未将溢洪道开挖考虑在内，致使计算结果都呈现出偏大的特征。
值得说明的是，溃坝模型对于坝体的物质组成参数物料系数反映比较敏感，类似于唐家山这种土体含量较多的堰塞坝体，其物料系数取值建议为3.65，并可根据土石比做出调整。而对于在物质组成和形成条件上比较类似的天然土石坝体，可以采用类比的方法，首现根据已知溃决资料推算出比较可靠的物料系数，然后可以通过模型预测其他类似的天然土石坝体得溃决过程。
