赵高文 乔建平 姜元俊 孟华君 王萌 杨宗佶

摘要：基于离散元法（ DEM）对滑坡堰塞坝的几何特征进行了数值模拟，考虑沟道断面形态、滑坡速度以及沟床坡度三种因素的影响，共27种组合。结果表明：沟道断面形态、滑坡速度对堰塞坝横向、纵向几何形态都有影响，而沟床坡度只对堰塞坝纵向几何特征产生影响;滑坡速度增大，堰塞坝纵剖面顶宽和底宽的比值减小，堰塞坝的有效厚度减小;沟床坡度增大，堰塞坝下游坝坡内角减小，上游坝坡内角增大，两者的比值减小，堰塞坝坝坡倾角的变化规律与内角的变化规律相反。

关键词：滑坡;堰塞坝;泥石流;漫顶破坏;离散元

中图分类号：P694;P642.22

文献标志码：A

doi：10.3969/j.issn.1000- 1379.2019.05 .003

滑坡堰塞坝属于堰塞坝的一种，是滑坡造成沟道或河道堵塞而形成的天然坝，该类坝体稳定性很差。Costa等[1]统计分析了73座滑坡堰塞坝（堵江型）破坏案例，发现其中27%在形成后一天内破坏、50%在形成后10 d内破坏，破坏的主要原因是漫顶侵蚀。柴贺军等[2]在Costa等的研究基础上，增加24个破坏案例后进一步指出.21%的滑坡堰塞坝在形成后一天内破坏.48%在形成后10 d内破坏，78%在形成后半年内破坏，88%在形成后一年内破坏。滑坡堰塞坝若堵塞河道形成大型堰塞湖，可能会淹没上游农田与民宅，若溃决形成洪水或泥石流，往往给下游带来巨大损失。例如：1786年6月四川泸定摩岗岭滑坡堰塞湖溃决，溃决洪水不仅摧毁了下游一切基础设施，还造成了近10万人死亡[3]：1933年叠溪滑坡堰塞湖溃决，造成下游近2 500人死亡[4]。

目前，对滑坡堰塞坝几何特征的研究较少且多以经验性分析为主：Cui等[5]、Xu等[6]、Yin等[7]根据堰塞坝高度、库容及坝体材料特性，对汶川地震灾区堰塞坝的危险性进行评估，并将其分为高危、危险、一般、低危等类型：王珊珊等[8]和Zhou等[9]基于GIS技术提取河谷横断面特征参数，评价滑坡堵江的可能性;Casagli和Ermini采用地貌指数法，建立了以堆积指数I_b、蓄水指数I_i[10]和无量纲堆积指数DBI[11]为控制因素的指数理论，用来评价堰塞坝的稳定性;Komp[12]、Nash等[13]、Dong等[14]、Dal Sasso等[15]均采用地貌指数法进行了相应的研究。除经验性分析方法外，数值模拟是较为常见的研究滑坡堰塞坝几何形态的方法：Hun-gr[16]根据流变理论，建立了滑坡堰塞坝的动力学数值模型;Lo等[17]基于离散元法（ DEM），反演分析了台湾小林村滑坡堰塞坝可能的堆积形态：Zhao等[18]采用CFD-DEM耦合的方法模拟滑坡堵江形成堰塞坝的过程。但现有关于滑坡堰塞坝几何特征的研究仍然存在以下两个问题：

（1）现有经验分析方法主要是根据已形成堰塞坝的几何特征，分析堰塞坝未来的发展趋势，并不能对堰塞坝的几何特征作出预测。有的经验分析准确性不够，例如，KOILip[12]采用地貌指数法分析新西兰境内滑坡堰塞坝稳定性时，结果并不理想。

（2）数值模拟研究不成熟。Hungr的研究并不适合滑坡体含水量较低、未发生液化的情况;Lo等的研究采用了小林村所在区域的真三维地形，这增加了模拟结果的真实性和可靠性，但由于地形过于特殊化和具体化，这种研究往往只能用于具体案例的反演分析，研究结果并非普遍适用，而且真三维地形模拟需要准确的地形数据（包括滑坡前和滑坡后）;Dong等[19]指出，准确的数字地形数据在绝大多数情况下是无法短时间内获得的：Zhao等的研究主要关注流固耦合的作用，未考虑地形条件等外在因素的影响，且该研究设置的离散元颗粒生成方式与野外真实的滑坡堰塞坝形成过程有较大区别。

1 堰塞坝几何特征与溃口的关系

几何特征对滑坡堰塞坝的溃决破坏有重要影响，在地貌指数法中，滑坡堰塞坝的高度、长度、宽度和体积等变量常常作为堰塞坝稳定性分析的关键因素[19]。除了用于堰塞坝的稳定性分析外，滑坡堰塞坝的几何特征还是研究堰塞坝的溃决机理和过程以及开展模型试验的基础。由于滑坡堰塞坝的几何形态很不规则，因此不同滑坡堰塞坝漫顶时的溃口位置不同，如著名的叠溪和红石岩滑坡堰塞坝，如图1所示。

绝大多数滑坡堰塞坝的破坏方式为漫顶侵蚀，由图1可见，堰塞坝的几何特征决定了堰塞坝漫顶时的溃口位置，溃口位置又影响漫顶时坝体的侵蚀过程。当溃口位于坝体侧面时，坝体受到的侧向侵蚀为单侧侵蚀，当溃口位于坝体中间时，侧向侵蚀为双侧侵蚀。这种因堰塞坝几何特征引起的侵蚀差异对于滑坡堰塞坝的溃决过程和输沙过程有重要意义。笔者曾做过验证性试验，模拟侧向侵蚀差异对溃坝过程的影响，两个坝体除初始溃口外，其他条件均相同，结果如图2所示，可以看出，双侧侵蚀的效果明显强于单侧侵蚀，进而导致溃口的扩展速度和输沙浓度等参数也不相同。

综上所述，滑坡堰塞坝的几何特征是堰塞坝的稳定性分析、破坏侵蚀过程分析以及开展模型试验的基础，研究滑坡堰塞坝的几何特征具有重要意义。

滑坡堰塞坝的形成过程是一个复杂且有待深入了解的山坡一沟谷耦合过程[12.20].很多因素都可能对滑坡堰塞坝的几何特征产生影响，包括滑坡类型、滑坡的运动以及堆积过程等[19]。在滑坡堰塞坝的所有诱发因素中，地震是最主要的诱发因素之-[l]，汶川地震后，整个灾区地层出现松动，形成了上千个不稳定斜坡[21]，很多不穩定斜坡失稳后都有可能形成新的滑坡堰塞坝，这给次生地质灾害的发生提供了条件。如果能预测这些滑坡堰塞坝的几何特征，将为灾害识别和应急抢险工作提供有力参考。

2 研究方法

采用离散元法（ DEM）对滑坡堰塞坝的几何特征进行初步分析，颗粒接触采用Hertz-Mindlin with JKR模型。

离散元法通过力一位移法则计算颗粒间的接触力，将颗粒运动当作颗粒振动进行求解，并根据牛顿第二定律分析颗粒运动（考虑重力的影响），其控制方程为

2.2 Hertz-Mindlin with JKR接触模型

Hertz-Mindlin with JKR模型是将Hertz-Mindlin无滑动模型与JKR（Johnson-Kendall-Roberts）黏聚力模型[22]相结合的一种颗粒接触模型，该模型的优点是可以考虑黏聚力的影响，模型中颗粒切向接触力、法向阻尼力和切向阻尼力的计算方法与Hertz-Mindlin无滑动接触模型的计算方法相同[23]。但由于考虑了颗粒表面能（y）的影响，因此颗粒的接触半径（a）比Hertz-Mindlin无滑动模型中的接触半径（ao）要大。法向力（ F_JKH）和重叠量（δ，可以为负值，表示颗粒未发生物理接触）的计算公式为

当表面能γ=0时，Hertz-Mindlin with JKR模型即为Hertz-Mindlin无滑动接触模型。

在岩土工程数值计算领域，Hertz-Mindlin with JKR模型可以考虑土体含水量和黏性颗粒的影响，因此具有较高的适用性。

3 数值模拟

3.1 物理模型

一般地，滑坡堰塞坝的动力过程包括：坡体失稳一加速运动一减速堆积等过程，如图4（a）所示。由于数值模拟考虑所有过程的计算量过大，因此有必要进行物理模型概化，本研究不考虑坡体启动和加速运动过程，通过API（Application Programming Interface）编程直接赋予坡体相应的速度，模拟滑坡体以某个速度滑人沟道并堆积成坝，滑坡体的形态也设为较规则的块体，如图4（b）所示。

3.2 影响因素

滑坡堰塞坝的几何特征受多种因素影响，笔者在此建议将各类影响因素按能空间尺度和能量尺度进行分类，以便今后开展滑坡堰塞坝动力过程关键因素的提取。

从空间尺度和能量尺度而言，沟道断面形态、滑坡速度（ν）以及沟床坡度（θ）是决定滑坡堰塞坝几何特征的基本参数，本研究主要考虑这三种因素对滑坡堰塞坝几何特征的影响。

3.3 模拟条件

由于滑坡的规模变化非常大，从几万到几百万甚至上千万立方米，因此本研究采用单位体积分析，滑坡体总体积0.5 m³.考虑到今后开展室内模型试验（滑槽模拟）的需要，此处将滑槽与沟道统称为几何体，其材质按钢材建模。滑坡体位于一矩形滑槽内，尺寸为1.0mxl.0 mx0.6 m（长×宽×高），沟道断面形态包括矩形、梯形和V形3种（见图5），为保证滑坡物质能够完全堵断沟道，3种沟道的断面面积与滑槽断面面积相等（ 0.6 m²）。同时，为保证不同沟道内的计算结果具有可比性，沟道高度保持不变（0.6 m），即滑坡体的重力势能保持不变，滑槽与沟道的夹角保持为30°，沟道纵向长度为5.0 m，滑坡体滑人沟道的角度与滑槽的底部平行。

滑坡堰塞坝的动力过程实质上是能量的释放和耗散过程，当势能保持不变时，坡体动能是影响滑坡堰塞坝动力过程的根本性因素。本研究设置的滑坡体滑人沟道的速度为单位速度（1 m/s）的整数倍，分别为1、2、3 m/s，代表低速、中速、高速（相对于沟道尺寸而言）三种条件，动能比值为1：4：9。沟床坡度设置为00、50和100三种，其中：00针对底床坡度较小的沟道或河道（统称沟道），50和100针对底床坡度较大的普通沟道。

3.4 参数标定

岩土工程数值计算对输人参数较敏感，为保证模拟结果的可靠性，必须对输人参数进行标定。根据2.2节的介绍.Hertz-Mindlin with JKR模型中，影响颗粒接触力的主要因素有表面能、静摩擦因数、滚动摩擦因数、恢复系数等。本研究以验证性试验（图2）中的土体材料作为标定对象，该土体材料取自都江堰虹口乡银洞子沟滑坡堰塞坝，其力学参数已通过大型直接剪切试验测出，黏聚力（c）和内摩擦角（φ）分别为6.2kPa和32°。由于土体黏性较低，因此参数标定方法以数值计算颗粒的自然堆积角（p）等于土体内摩擦角时的参数组合作为数值计算的输人参数。参数标定结果见表1.所有输入参数见表2.

4.1 沟道断面形态的影响

因27种组合的几何参数多达243个，此处仅以高速条件（相对于沟道尺寸而言）下，平直沟道内滑坡堰塞坝的几何特征为例进行分析。该条件下不同沟道内滑坡堰塞坝的几何数据见表3。

由表3可见，高速条件下.V形沟道内的堰塞坝纵剖面底宽最大，梯形次之，矩形最小;与此规律相反，堰塞坝的坝坡倾角和内角以矩形最大，梯形次之.V形最小。上述规律说明：V形沟道内滑坡堰塞坝向沟道上、下游延伸堆积的范围最大。同时，由表3还可看出，高速条件下，矩形沟道内堰塞坝的坝顶宽度L_t=0，这说明堰塞坝的纵剖面为三角形，但梯形和V形沟道内堰塞坝纵剖面均为梯形。低速条件下，三种沟道内的堰塞坝纵剖面如图7所示，全为梯形。由此可见，沟道断面形态对堰塞坝的几何形态有直接影响。

4.2 滑坡速度的影响

此处仅以沟床坡度θ=0°时梯形沟道内堰塞坝的几何特征为例进行分析，不同滑坡速度条件下堰塞坝的几何数据见表4。

低速条件形成的堰塞坝，h_max與h_min的比值达到4：1，说明堰塞坝的横断面存在很大高差，这类堰塞坝比较常见，如Young河滑坡堰塞坝（图8（a））和都江堰虹口乡银洞子沟滑坡堰塞坝（图8（b））均属该类型，其中银洞子沟滑坡堰塞坝刚形成时的高差接近160m。横断面存在巨大高差的堰塞坝若发生漫顶溃决，则其溃口侧壁失稳滑动的规模会越来越大。

隨着滑坡速度的增大，堰塞坝横断面的高差h_max-h_min逐渐减小。由不同速度条件下梯形沟道内堰塞坝的横断面图（如图9所示）可知，随着滑坡速度增大，滑坡体可以向对岸爬高而形成马鞍状的断面形态，这与图1中的两个真实案例相吻合。由此可见，滑坡速度对堰塞坝漫顶溃决时的溃口位置有重要影响，低速条件下形成的堰塞坝，初始溃口位于坝体侧面，侧向侵蚀属单侧侵蚀（如图l（a）、图8所示），但当滑坡速度增大后，堰塞坝的初始溃口向中部移动，侧向侵蚀变为双侧侵蚀（如图l（b）所示）。

随着滑坡速度的增大，堰塞坝顶宽（L_t）逐渐减小，堰塞坝顶宽与底宽（L_b）的比值亦减小。若以（L_t+L_b）/2表示堰塞坝的有效厚度，则随着滑坡速度的增大，堰塞坝的有效厚度变小。

4.3 沟床坡度的影响

通过对试验结果的分析可知，沟床坡度主要对堰塞坝的纵向堆积产生影响。27种组合对应的坝坡倾角和内角见表5。

4.3.1 坝坡内角（φ'_u、φ'_d）

根据表5的统计结果，可得出以下认识：

（1）沟床坡度一定时，堰塞坝坝坡内角的大小关系是矩形沟道≥梯形沟道≥V形沟道。V形沟道内的堰塞坝坝坡内角大小不受滑坡速度的影响。

（2）整体而言，沟床坡度增大，堰塞坝下游坝坡内角（φ'_d）减小，上游坝坡内角（φ'_u）增大，φ'_u与φ'_d的比值减小，如图10所示。

5 结论

基于离散元法（ DEM）对滑坡堰塞坝的几何特征进行模拟，采用可以考虑颗粒黏性的Hertz-Mindlinwith JKR模型进行了27种参数组合的模拟分析，得到以下几点主要认识。

（1）滑坡速度和沟道形态对滑坡堰塞坝的横向和纵向几何形态都有影响，而沟床坡度主要影响滑坡堰塞坝的纵向几何特征。

（2）低速条件下形成的堰塞坝在横向存在较大高差，堰塞坝若发生漫顶溃决，初始溃口位于堰塞坝侧面，侧向侵蚀为单侧侵蚀。当滑坡速度达到高速条件后，滑坡体可以向对岸爬高而形成马鞍状的横断面，这类堰塞坝若出现漫顶溃决，初始溃口位于坝体中部，侧向侵蚀为双侧侵蚀。

（3）滑坡速度增大，堰塞坝纵剖面顶宽（L_t）和底宽（L_b）的比值减小，堰塞坝的有效厚度减小。

（4）沟床坡度增大，堰塞坝下游坝坡内角减小，上游坝坡内角增大，两者的比值减小。但是堰塞坝坝坡倾角的变化规律和内角的变化规律相反。

参考文献：

[1]COSTA J E，SCHUSTER R L.The Formation and Failure ofNatural Dams[J].Ceological Society of America Bulletin，1988， 100（7）： 1054-1068.

[2] 柴贺军，刘浩吾，刘汉超，等，天然土石坝稳定性初步研究[J].地质科技情报，2001， 20（1）：77-81.

[3]DAI F C，LEE C F，DENC J H，et al.The 1786 Earthquake-Triggered Landslide Dam and Subsequent Dam-Break Flood onthe Dadu River， Southwestem China [J]. Geomorphology，2005， 65（3）：205-221.

[4] 王兰生，杨立铮，李天斌，等，四川岷江叠溪较场地震滑坡及环境保护[J].地质灾害与环境保护，2000， 11（3）：195-199.

[5]CUI P，ZHUY Y，HAN Y S，et a1.The 12 May WenchuanEarthquake-Induced Landslide Lakes： Distribution and Prelim-inary Risk Evaluation[J]. Landslides， 2009，6（3）：209-223.

[6] XU Q，FAN X M， HUANC R Q，et al.Landslide DamsTriggered by the Wenchuan Earthquake， Sichuan Province，South West China[J].Bulletin of Engineering Ceology andthe Environment， 2009， 68（3）： 373-386.

[7]YIN Y P，WANC F W， SUN P.Landslide Hazards Triggeredby the 2008 Wenchuan Earthquake， Sichuan， China[J]. Land-slides， 2009， 6（2）： 139-152.

[8] 王珊珊，童立强，基于河谷横剖面形态特征的滑坡体堵江易发性评价研究[J].地理与地理信息科学，2016， 32（5）：97-102.

[9]ZHOU L，TAKASHI 0.Longitudinal and Transverse Profilesof Hilly and Mountainous Watersheds in Japan[J].Geomor-phology， 2009， 111（1-2）： 17-26.

[10]CASACLI N，ERMINI L.Ceomorphic Analysis of LandslideDams in the Northem Apennine[J].Japanese Ceomorpho-logical Union， 1999， 20（3）：219-249.

[11]ERMINL，CASACLIN.Prediction of the Behavior of Land-slide Dams Using A Ceomorphologcal Dimensionless Index[J]. Earth Surface Processes and Landforms， 2003， 28（1）：31-47.

[12]KORUP 0.Geomorphometric Characteristics of New ZealandLandslide Dams[J].Engineering Geology， 2004， 73（1）：13- 35.

[13]NASH T， BELL D， DAVIES T， et al. Analysis of the For-mation and Failure of Ram Creek Landslide Dam， South Is-land， New Zealand [ J]. New Zealand Joumal of Ceologyand Ceophysics ， 2008， 51（ 3） ： 187-193.

[14]Dong J J， Tung Y H， Chen C C， et al. Logistic RegressionModel for Predicting the Failure Probability of A LandslideDam[J]. Engineering Ceology， 2011， 117（ 1） ： 52-61.

[15]DAL SASSO S F， SOLE A， PASCALE S， et al. AssessmentMethodology for the Prediction of Landslide Dam Hazard [J] .Natural Hazards and Earth System Sciences， 2014， 14（3） ：557-567.

[16]HUNCR 0. A Model for the Runout Analysis of Rapid FlowSlides， Debris Flows， and Avalanches [ J]. CanadianCeotechnical Journal， 1995 ， 32（4） ： 610-623.

[17]10 C M， LIN M L， TANC C L， et al. A Kinematic Modelof the Hsiaolin Landslide Calibrated to the Morphology ofthe Landslide Deposit [ J]. Engineering Geology， 2011，123（1） ： 22-39.

[18]ZHAO T， DAI F， XU N. Coupled DEM-CFD Investigationon the Formation of Landslide Dams in Narrow Rivers [J] .Landslides， 2017， 14（1） ： 189-201.

[19]DONC J J， LAI P J， CHANC C P， et al. Deriving Land-slide Dam Geometry From Remote Sensing Images for theRapid Assessment of Critical Parameters Related to Dam-Breach Hazards [ J] . Landslides， 2014， 11（ 1） ： 93-105.

[20]COSTA J E， SCHUSIER R L. Documented Historical LandslideDams from Around the World[ R] . Virginia ： US Ceologcal Sur-vey， 1991：91-239.

[21]殷躍 .汶川八级震地质灾害研究 [J] .工程地质学报，2008， 16（4） ： 433-444.

[22] JOHNSON K L， KENDALL K， ROBERTS A D. Surface Energyand the Contact of Ilastic Solids [ C] //Proceedings of the RoyalSociety of London A： Mathematical， Physical and Engneering Sciences. London ： The Royal Society， 1971： 301-313.

[23] SANTOS K G， CAMPOS A V P， OLIVEIRA O S， et al.Dem Simulations of Dynamic An～e of Repose of AcerolaResidue： A Parametric Study Using a Response SurfaceTechnique [J] . Blucher Chemical Engineering Proceedings ，2015 ， 1（ 2） ： 11326-11333.

[责任编辑张帅]