基于FLO-2D的都江堰市龙池镇黄央沟泥石流数值模拟

王纳纳，唐川

(成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都 610059)

基于FLO-2D的都江堰市龙池镇黄央沟泥石流数值模拟

王纳纳，唐川

(成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都 610059)

利用HEC-HMS软件结合ArcGIS技术模拟都江堰市龙池镇黄央沟2010年8月13日泥石流的清水流量过程线,然后生成的清水过程线输入到FLO-2D软件中模拟了黄央沟泥石流的运动和堆积过程，计算了堆积的深度和分布范围。文中分析比较了模拟结果与实际调查结果，误差较小。可见基于FLO-2D的泥石流规模预测具有一定的应用意义，预测结果可以为泥石流的防治和预警提供依据。

黄央沟；泥石流；HEC-HMS；FLO-2D

泥石流是一种固、液两相的连续流体，由固相物质(土、砂和石块)和液相物质(水)组成[1]。它具有暴发时间的突发性、短暂性、周期性和季节性，运动过程的直进性、阵流性、非定常性、不均匀性和侵蚀性等特点。泥石流的运动特征(流速、流量、流深等)随着发展过程而变化，其动力学过程相当复杂。

在泥石流灾害的研究方面，国内外学者已经做了大量工作。泥石流的研究方法主要有野外调查、实验研究、数值模拟。数值模拟因成本低，重复性强，近年来被广泛应用到泥石流研究中。堆积区是泥石流致灾的主要地区，研究泥石流堆积的数值模拟方法，定量分析预测泥石流灾害的发生过程及结果，具有重要的实际意义。对堆积区研究最多的是日本学者。高桥堡提出以质量守恒在动量守恒为基础，建立二维非恒定流方程，采用有限差分思想，并应用于水工模拟[2]；石川芳治等通过水工模型实验用数值模拟研究了泥石流的堆积泛滥范围[3]。国内学者以唐川、刘希林做的工作最多，他们应用剖开算子有限差分方法进行求解，并借此对云南怒江州芭蕉河泥石流等实际的泥石流进行了计算模拟，得到了泥深及流速的分布，对泥石流的危险度进行评价[4-7]。罗元华考虑了侵蚀因素的控制方程，进行泥石流堆积数值模拟[8]。余斌用SIMPLER法求解离散方程，计算了河床底部有凸台情况下的泥流速度分布[9]。台湾学者詹钱登等也做了较深入的研究。他们都是以高桥堡的数学模型为基础，从质量守恒原理和动量守恒方程推导出来的，所不同的是阻力项的确定和运动形式的简化各有不同，计算机数值模拟方法各有特色。

本文重点使用FLO-2D软件对都江堰市龙池镇黄央沟泥石流运动堆积过程进行模拟。FLO-2D是由O’Brien提出的二维洪水与泥石流数值模拟软件，可用于模拟洪水灾害和风险评估，能够推算流体运动的流速及堆积深度，合理推估淹没面积。

1 研究区概况

黄央沟位于都江堰市龙池镇南岳村1组，属于龙溪河流域。沟口地理坐标为东经103°30′40″，北纬31°15′20″，沟长为2.1 km，平均坡降为420 ‰，汇水面积1.0 km2。黄央沟所在地地质构造复杂，断裂发育。最高海拔为1 900 m，最低为945 m，是典型的中山地貌，沟道两侧均为陡峭斜坡[10]。

2010年8月13日，都江堰市龙池镇遭遇特大暴雨，16时降雨量突然增大，16∶00～17∶00的小时降雨量达75 mm(图1)。达到20年一遇洪水量，龙溪河流域遭遇山洪泥石流袭击，造成龙溪河河床整体抬高近5 m，对龙溪河流域“5．12”地震灾后重建带来严重影响，经济损失巨大。

图1 都江堰市龙池镇“8.13”降雨数据图[11]Fig.1 Dujiangyan Long Chi Zhen“ 8 ． 13” rainfall data in figure[11]

2 基于FLO-2D的泥石流数值模拟

2.1 FLO-2D计算模型假设、限制条件

2.1.1 假设条件[12]

(1) 假设泥石流运动属于浅水波模式。

(2) 假设差分时间间隔内为定常流。

(3) 假设泥石流沟道断面及粗糙度比较规则。

(4) 假设泥石流液体压力分布为静水压力分布。

(5) 假设每个网格只能由单一的高程值及曼宁系数值。

2.1.2 限制条件[13]

(1) 因计算过程中假定为定沟床模式，故无法模拟沟道的侵蚀现象。

(2) 无法模拟激波与水跃现象。

2.2 FLO-2D运动控制方程

泥石流是一种多相体，包含了石块、砾石、泥沙、水等物质，在运动过程中这些物质混杂为一体，具有流体的力学性质，其复杂的流动过程同样遵循流体运动的连续性原理和动量守恒原理。

FLO-2D软件利用牛顿体模式与有限差分法求解二维模式的连续方程和动量方程，可以得到泥石流体的流速和泥深。其连续性方程和动量方程如下：

连续方程[12]：

(1)

动量方程[12]：

(2)

(3)

式中,h是泥深；Sf是摩擦坡降；S0是沟床坡降；g是重力加速度；t表示时间；u表示x方向的平均流速；v表示y方向的平均流速；i表示有效降雨强度。

2.3 FLO-2D模拟所需的参数

2.3.1 地形参数

地形数据为四川省龙池地区1∶50 000等高线和高精度遥感影像图。因FLO-2D需要ASCII格式的数字高程模型(DEM)或高程数据( Elevation data )，可以通过ArcGIS软件获得。具体过程如下：

DEM获取：利用ArcGIS里的3D Analyst扩展模块里的Create TIN From Features 功能，可以用等高线创建TIN(不规则三角网)。再将TIN文件栅格化，利用Tin to Raster功能得到研究区的DEM。

ASCII格式的转化：利用ArcGIS里的Arc Toolbox工具箱中的Raster to AscII功能将DEM转化成ASCII格式的高程文件。

将高程文件导入FLO-2D后，建立10 m×10 m的计算网格，按流域范围确定计算范围，并对网格进行高程插值，完成地形数据的处理。

2.3.2 重度

泥石流的固相物质一部分来自形成区，另一部分来自泥石流沟道两岸滑坡堆积物及冲刷沟床形成的床面物质。泥石流在运动过程中会随着地形的起伏而产生淤积变化，其物质组成也随之发生变化。在沟道中不同位置取得的沉积物样品的颗粒组成也可能不同。根据野外调查可初步判断属于粘性泥石流沟，我们在粘性泥石流的堆积区和流通区分别取样，计算黄央沟的沉积物重度。黄央沟沉积物的重度为1.67 g/cm3。

2.3.3 流量计算

本文中清水流量过程线选用美国工程兵水文工程中心(Hydrologic Engineering Center，HEC)所研发的HEC-HMS软件及其系列软件GeoHMS进行模拟。HEC-HMS模型包括流域模型、气象模型、控制模型等三大模块，如图2所示。HEC-HMS水文模型内嵌有多种产流、汇流、基流以及河道洪水演进的计算模型，能够模拟多种环境条件下流域的次降雨径流过程。GeoHMS是一基于GIS的ArcMap的软件包。模块作为HEC-HMS的数据接口，对输入的数字高程模型DEM进行地形处理，通过分析数字地形信息，生成HEC-HMS接受的格式文件后，把GeoHMS的水文结果导入到HEC-HMS中，以进一步进行清水过程线的模拟[14]。

图2 HEC-HMS模型基本构成Fig.2 HEC - HMS model

首先，利用GeoHMS生成可导入HEC-HMS的.hms文件。具体过程如下：根据黄央沟泥石流的实测地形建立黄央沟泥石流的数字高程(DEM)，并选取两条小支沟的交汇处作为泥石流数值模拟的起始点(即图3中的W140)。为了避免数字高程模型DEM上存在的洼地和平台影响沟道流域的提取，对龙池地区的DEM要进行填洼处理，然后确定水流流向，计算汇流累积流量，提取流域河网即沟道，定义流域出口位置，获取整个流域的地形特征和沟道特征参数等等，即可生成数字小流域模型.hms文件(图3)。

图3 Geo-HMS构建的basin模型Fig.3 The basin model of Geo - HMS construction

其次，将生成的流域模型导入到HEC-HMS软件。根据龙池地区的降雨径流条件、地形条件、气象条件、土地利用情况，产流计算模型选定SCS曲线法，汇流计算模型为Clark单位线法；河道洪水演算模型为Muskingum Cunge法。由于洪水的历时相对较短，加之流域面积小，因此，本次模拟忽略基流对径流的影响。利用ArcGIS技术，根据研究区DEM、土地利用情况和土壤类型设置CN值，由GEO-HMS预处理获得的沟道的基本特征数据计算集中流动时间、蓄水常数。参数设置完成便可得到如图4所示，FLO-2D设置的泥石流启动点处HEC-HMS所演绎出的“8.13”整个降雨过程的降雨损失、降雨深度和清水流量过程线图。

根据黄央沟泥石流(2010年8月13日16∶30～18∶00)的清水流量过程线，结合泥石流的的流量计算公式，可以得出泥石流的流量值，获得泥石流的运动过程(图5)。由图5可以看出，16∶20左右水量开始增大，16∶56左右水量达到最大值，随后开始减弱，持续到18∶00左右，最大水流量为6.5 m3/s。

泥石流流量计算：粘性泥石流在运动的过程中会对沟道两侧产生侵蚀，泥石流的流量会产生放大效应。FLO-2D二维泥石流模拟中，泥石流流量是以清水流量乘以体积膨胀因子BF，计算公式为BF= 1/ (1-Cv)，Cv为体积浓度。泥石流的流量过程线见图5。

2.3.4 体积浓度参数

泥石流体在流动过程中，随着地形坡度、地表粗糙程度的不同，体积浓度不同。由于现场取样往往很难取到或者很难配到暴发时的泥石流体，因此在本研究的模拟中，用FLO-2D使用手册[12]建议值，不同部位大致在0.58～0.62。

2.3.5 流变参数

流变学所研究的就是流动、变形与应力间的关系。为了定量地描述这些关系，即为了表达泥石流体流变性能，需要定义一些参数。我们称这些参数为流变参数泥石流体有很多流变模式，而FLO-2D基于的是O'Brien-Julien模式即

(4)

式中,τc为粘性屈服应力;τmc为莫尔-库仑剪应力;τv为粘性剪应力;τt为紊流剪应力;τd分散剪应力。

按Meyer-Peter、Müller (1948) 和Einstein (1950) 提出的思路重新定义(4)式中的5个剪应力项，将其写成无量纲的坡降形式，如下：

(5)

式中,Sf是摩擦坡降;Sy是屈服坡降；Sv是粘性坡降；Std是紊流-分散坡降。

图4 启动点处“8.13”降雨损失及深度图、清水流量过程线图Fig.4 Starting point “8.13”rainfall loss and depth map, water flow process chart

图5 “8.13”黄央沟的清水流量、泥石流流量过程线Fig.5 “8.13” Huang Yang ditch water flow, debris flow hydrograph

这样，将三项代入(5)，最终可以得到：

(6)

式中,η=α1eβ1·Cv和τy=α2eβ2·Cv。

其中α1，β1和α2，β2的取值可以按王裕宜、詹钱登相关文献[15]确定，本次模拟的流变参数取值α1= 0.000 247，β1=15.48，α2=0.03，β2=14.42；层流阻滞系数K按Woolhiser(1975)建议值，取K=228 5。

2.3.6 曼宁系数

曼宁系数需要为每一计算网格赋值，默认系数为0.04。具体值可以根据表1取值，也可以取当地的经验值。本文按王裕宜、詹钱登等人提出的不同阻力介质状态下的统一公式(7)计算[15]，公式如下：

(7)

式中,Rns为体积浓度;h为泥深。

根据现场调查黄央沟沟的堆积扇厚度为5 m左右。由体积浓度和泥深利用公式(7)可以得到黄央沟堆积扇的曼宁系数为0.08～0.12。

2.3.7 模拟总时间的确定

根据现场走访调查，居民反映黄央沟泥石流运动的总时间大致为90 min，与当天的降雨强度峰值出现的时间大致相当，见图5。

2.4 模拟结果及误差分析

通过对黄央沟泥石流运动堆积过程的模拟，可以获取泥石流整个流动过程中各个时刻，每个格网的最大流深分布。模拟结果见图6、图7。

图6 黄央沟的堆积范围模拟结果图Fig.6 Accumulation of yellow central ditch simulation results map

图7 黄央沟堆积扇厚度模拟结果图Fig.7 Huang Yang Gou accumulation simulation results Fig. fan thickness

图8 黄央沟泥石流沟口堆积扇Fig.8 Huang Yang gully Mizoguchi fan

模拟结果显示，黄央沟“8.13”泥石流堆积扇厚度最大为4.8 m，位于堆积扇中部，最远冲出距离达到202 m，最大宽度120 m。从泥石流数值模拟结果(图6、图7)中可以看出，泥石流冲出山口向两边扩散，为典型的扇形堆积扇，主流沿堆积扇中央冲入龙溪河。泥石流流深除了在沟口处较大之外，在堆积扇最前缘泥石流入汇主河处也可达到2.0～3.0 m。这些结果都与野外调查相较吻合。从流深的大小分布和堆积扇的空间堆积状态可以验证数值模拟结果的合理性。遗憾的是野外调查距离发生泥石流时间较长，没有精确的现场资料，无法将数值模拟与野外情况进行更为详细的比较验证。

误差分析：用相对误差验证计算结果的准确性，计算方法是用残差值除以调查值的百分比。残差值为实际调查值与模拟计算数值的差，计算结果见表1。

(8)

表1 “8.13”黄央沟泥石流模拟结果与实际比较

Table 1 “8.13” yellow central debris flow simulation results and the actual comparison

最大泥深最大堆积长度最大堆积宽度实际情况模拟结果误差5m4.8m-4%207.4m202m-2.6%180m120m-33.3%

对于表1中的误差值均为负数，表示两条沟的模拟结果都小于与实际调查结果。原因是模拟计算没有考虑运动过程中泥石流体对沟道两侧的侵蚀和对沟道的下切影响，也没有考虑到河水对堆积体前缘的冲刷影响。侵蚀会使新的固体物质加入泥石流运动，加大运动规模；河水冲刷对使堆积扇前缘产生扩散，加大堆积长度和宽度。

从结果来看，黄央沟沟堆积扇的最大长度模拟精度高，但最大宽度模拟精度为-33.3%。原因是实际过程中，泥石流冲出物在进入龙溪河后，河水的冲刷改变了堆积扇前缘的堆积形态，持续而来的泥石流体补充在被冲刷的位置，这一过程持续进行直至堆积完结。而模拟过程没有考虑这一影响，这种情况下的计算出来的堆积宽度必然小于河水冲刷条件下堆积的宽度。

3 结论

黄央沟泥石流是龙溪河流域一条较大的泥石流沟，本文是对“8.13”发生的泥石流使用HEC-HMS和FLO-2D模拟黄央沟泥石流，得出如下结论：

(1) 本文选用HEC-HMS进行清水过程线的模拟，能够比较真实地反映黄央沟流域的产流、汇流以及河道洪水的演算过程，避免了概化过程线带来的误差。

(2) FLO-2D数值模拟的黄央沟泥石流堆积扇与实际调查，堆积扇的最大长度模拟精度高，但最大宽度模拟精度较低。说明堆积长度受河水冲刷的影响不大，敏感度较差；堆积宽度受河水冲刷的影响较为明显，敏感度较高。

(3) FLO-2D数值模拟的流深的分布、大小都在合理的范围之内。从模拟结果可以看出，泥石流主要在沟口堆积，沿流向堆积厚度逐渐变小，向两侧也逐渐变小。

(4) 模拟结果比实际野外调查的结果偏小，除了与泥石流冲出物被龙溪河冲刷带走外，还因为在采用FLO-2D在进行模拟的过程中，并没有将泥石流运动过程中对沟道的下切和冲刷作用考虑进去。

(5) 泥石流运动数值模拟不仅可以模拟泥石流堆积扇的堆积范围，获得泥石流流深等重要参数在泥石流堆积扇区的时空分布，确定泥石流的危害范围和危害程度，同时还可以检验泥石流减灾工程的减灾效果。数值模拟将是未来学习研究的一个重要的有效的、不可替代的方法之一。

[1] 康志成.中国泥石流研究[M].北京：科学出版社，1989：23-24.

[2] 高桥堡，中川一，佐藤宏章.扇状地における图砂泛滥灾害危险度の评价，京都大学防灾研究所年报，31(B-2).1988:655-676.

[3] 石川芳治，水山高久，井户清尾.堆积扇上泥石流堆积泛滥机理.泥石流及洪水灾害防御国际学术讨论会文集(成都)，A(泥石流).成都，1991: 27-31.

[4] 唐川，刘希林.泥石流动力堆积模拟和危险范围预测模型[J].水土保持学报，1993，5(3)：37-40.

[5] 唐川.泥石流堆积泛滥过程的数值模拟及其危险范围预测模型的研究[J].水土保持学报，1994，8(1)：45-50.

[6] 唐川.泥石流堆积扇危险度分区评价的数值模拟研究[J].灾害学，1994，9(4)：7-13.

[7] 唐川.平面二维泥石流数值模拟方法的探讨[J].水文地质工程地质，1994，21(5)：9-12.

[8] 罗元华.泥石流堆积数值模拟及泥石流灾害风险评估方法研究[D].武汉：中国地质大学，1998.

[9] 余斌.二维定常泥流的模拟.自然灾害学报，1995，4(4)：96-99.

[10]杜霏.都江堰市龙池镇黄央沟泥石流基本特征及成因分析[J].科协论坛，2012，9：113-114.

[11]常鸣，唐川，付荣，等.水打沟泥石流形成条件及其静动力学特征[J].水电能源科学，2012，30(8)：104-106.

[12]O’Brien,J.S. FLO-2D Reference Manual.2009.

[13]Jan C D，Shen H A. Review of Debris flows Flow Analysis Proceedag XXV cong ress,I AHR[C].1993，3：25-32.

[14]蔡新明，董志勇，张永华.HEC系列水利软件的应用[J].浙江水利科技，2005，11:23-26.

[15]王裕宜，詹钱登，韩文亮，等.粘性泥石流体的应力应变特性和流速参数的确定[J].中国地质灾害与防治学报，2003，14(1)：09-13.

NUMERICAL SIMULATION OF HUANGYANG GULLY DEBRIS FLOW LONGCHI TOWN DUJIANGYAN CITY BASED ON FLO-2D

WANG Na-na，TANG Chuan

(State Key Laboratory of Geo-hazard Prevention and Geo-environment Protection ，Chengdu University of Technology，Chengdu 610059，China)

In this paper, the Clearwater flow hydrograph of Huangyang gully debris flow on August 13, 2010 in Long chi town Dujiangyan city was studied by using the HEC-HMS software combined with ArcGIS method, and then input the result to FLO-2D software to simulate the Huangyang gully debris flow and calculate the accumulation of depth and distribution range. The error between the simulated and actual result according to field investigation result is small. The prediction of debris flow scale based on HEC-HMS and FLO-2D software can be used in application significance, the prediction results can provide the basis for the debris flow prevention and prediction.

Huangyang Gully；Debris Flow；HEC-HMS；FLO-2D

1006-4362(2014)01-0107-06

2013-09-17 改回日期： 2013-10-17

科技基础性工作专项项目(编号：2011FY110100)，地质灾害防治国家重点实验室自主项目(SKLGP2001Z008)；四川省教育厅自然科学重点项目(11ZA047)

P642.23

A

王纳纳(1987- )，女，河南洛阳人，硕士研究生，研究方向为泥石流灾害评价与防治工程。