火山泥石流是伴随着火山喷发，由火山湖泊水的突然溢出或火山喷发造成火山顶冰雪融化所形成的突增水流，夹杂着火山碎屑物质、岩浆等，在重力作用下，沿着火山的沟谷或山坡一起快速向下流动的一种特殊洪流^[1]。根据Withan在2005年的统计，20世纪内已经有超过30 000人因为火山泥石流而丧生。随着火山地区的经济发展和人口增加，减轻火山泥石流灾害的研究已经越来越紧迫。长白山是目前公认的中国最具喷发危险的活火山之一，其火山口积聚着20.4×10⁸m³湖水，一旦天池火山喷发将可能会形成大规模的火山泥石流，对其下游27×10⁴多人口以及10余座水电站将造成巨大的破坏，严重威胁到该地区人民生命财产安全。火山地区由于自然条件、科技水平等的限制，对火山泥石流的直接测量和观测有时候是很难实现的，在这种情况下要对火山泥石流进行研究就需要借助于数值模拟。另外，通过数值方法成功的模拟出火山泥石流发生时的情况可以为政府部门提供防灾减灾的参考资料，有效的减轻火山泥石流灾害。因此火山泥石流的数值模拟被众多学者所重视，但火山泥石流的形成、流通和堆积是一个相当复杂的过程，目前对火山泥石流模拟的研究还处于起步阶段。近年来，Aguilera、Emmanuel、Alison等学者分别通过不同的方法对厄瓜多尔Cotopaxi火山、菲律宾Pinatubo火山、新西兰Ruapehu火山等进行了数值模拟，对这些地区的火山泥石流危险性进行了预测^[2-4]，这些学者的工作为火山泥石流模拟研究奠定了一定的基础，但是这些研究主要是对火山泥石流复杂的形成机理等进行研究，虽然具有很高的理论意义，但是在实际应用中不能被简单有效的利用本次研究并非模拟火山泥石流复杂的形成、流通和堆积过程，而是对长白山天池火山泥石流的影响范围、流通速度及流通时间等进行研究，因为以上几方面内容是确定灾害预警时间和灾害影响范围的依据，是制定火山泥石流防治方案的基础。研究主要是通过流体动力学模型，结合假设的火山喷发时形成火山泥石流的水源量，利用支持曲面格式的一种3D水动力软件计算系统对长白天池火山泥石流进行模拟，通过模拟结果预测天池火山一旦喷发，所形成的火山泥石流的影响范围、流通速度及流通时间等特性。这些研究成果可以为制定火山泥石流灾害防治方案提供参考资料，使减轻火山泥石流灾害从被动救灾转变为主动防灾。

长白山天池火山位于吉林省东部中朝边境，主峰白云峰海拔2 700 m，是中国东北地区最高的山峰。长白山天池位于白云峰顶部，是火山喷发形成的火山口湖，湖水面高程为2 190 m左右，湖水平均深度为204 m，湖面面积约为9.4 km²，总蓄水量约20.4×10⁸ m³。长白山天池是松花江、图们江和鸭绿江三江之源，天池北部的缺口形成一个落差68 m的瀑布，天池水从瀑布处跌水直泻流入二道白河河谷。本次研究的范围主要是长白山地区二道白河流域的上游，海拔高度在500~2 700 m之间，为典型的火山构造地貌，自下而上主要由熔岩台地、山麓倾斜熔岩高原和火山锥体3大地貌单位构成，围绕天池呈同心环状分布^[5]。区内主要分布新生代火山岩，只有很小面积出露的华力西期侵入岩、侏罗系和中元古界地层。从大地构造来看，该区域形成于中元古代至早古生代，基底为太古代和早元古代的变质岩系、混合花岗岩、混合岩，属于胶辽断块隆起带。新构造时期，研究区域山体整体发生间歇性抬升，从中更新世至今，区域以断块差异上升运动为主，形成二级夷平面，沿河流两侧形成多级阶地堆积，伴有多期玄武岩喷发。从新构造单元分区看，区域属于小兴安岭—长白山—千山隆起带断块差异运动上升区。

长白山天池火山泥石流主要沿天池下方的二道白河流域上游沟谷流通，在流通过程中沟谷两边由于崩塌、风化等形成的松散物质不断对火山泥石流进行物源补充，最后在二道白河流域两岸的平坦区域形成大量火山泥石流堆积。天池火山泥石流随着距离瀑布口的距离增大而显示出了明显的分带性，在距离天池瀑布口大约0~6 km范围内为火山泥石流的形成及流通区，在这个区段主要是泥石流的形成、流通和物源补充阶段。在距离瀑布口6~50 km范围内火山泥石流形成堆积，该区域的火山泥石流堆积属于岩屑流堆积^[6]。堆积物中的松散物质的粒径和厚度由天池火山口往下逐渐变小，在火山口下游10 km左右的堆积物中有粒径4~5 m的粗面岩和碱流岩，厚度约为15~20 m左右。到距离火山口50 km左右的区域，堆积物中岩石的最大粒径递减到2~3 m，岩性以粗面岩和玄武岩为主，堆积厚度在10~12 m左右。距离火山口大于50 km处所形成的堆积物属于消退流堆积^[6]，堆积物主要是粗—细粒的浮石、火山碎屑、晶屑等，堆积厚度为8 m左右。

本次研究拟采用的Delft 3D是支持曲面格式的一种3D水动力模型系统，能非常精确地对水动力、波浪、泥沙等进行计算。该系统目前已经实现了有限单元法计算，计算过程快速、稳定，并能保证质量、动量和能量守恒。

利用Delft 3D进行分析时，将研究区域离散为若干个正交的网格(如图 1)，对每个网格分别进行计算，最后得到流速等信息。

图 1(Figure 1)

图 1 Delft 3D计算网格划分

火山泥石流是一种典型的多相体，在计算过程中将其简化成理想化的Newtonian流体，其特性符合Navier-Stokes方程及Boissinsq假设。将计算网格和计算步长划分得足够小，可以认为在每一个计算网格中，火山泥石流的流动是连续的。基于上述假设，得到火山泥石流水平面上的连续方程

式中：h是高于坝体的水位深度(m); d是坝以下水的深度(m); H=h+d是总的水位深度(m); t是时间(s), ${{G_{xx}}}$和${{G_{yy}}}$是Delft 3D中曲线转换矩形的调整系数，x和y是曲线水平坐标，U和V分别是x和y方向的平均速度(m/s)，Q是每个网格的流量(m³)。

式中：q_in和q_out分别为单位时间内每个网格中流体的流入和流出量(m³)，是垂直方向上水位变化幅度$\sigma = \frac{{z - h}}{H}$。

考虑到速度在垂直方向的变化，将式(2)改写为

式中：u、v和ω分别为水平和垂直方向上的速度(m/s)(如图 1)。

火山泥石流的形成机制与火山活动密不可分，火山爆发使积聚的热能瞬间转变成巨大的动能，冲破地表和火山山体，将其附近的巨大岩石、土体粉碎。岩浆混合着火山湖泊中的水以及由于火山喷发的高温所融化的火山顶的冰雪夹杂着这些岩石、土体碎块一同向火山口处的缺口及低处快速流动，形成了火山泥石流。同时，山体或附近地表被火山地震所振动，有大量崩塌和滑坡发生，为火山泥石流的形成提供了大量固体物源^[7-13]。火山泥石流实际上就是一种特殊的泥石流，也可以分为形成区、流通区和堆积区3个区段。长白山天池火山泥石流的形成区就是山顶的火山湖泊—天池。根据以往历史记载和考察资料等，长白山天池火山喷发后，泥石流主要是沿着二道白河流域发育，因此选取天池瀑布口以下12 km左右的山谷区作为本次研究的流通区，选取流通区下游的平坦地带为堆积区(见图 2)。火山泥石流的发展在本次模拟中是通过流体的波状传播将初始的“干”地区转变到“湿”地区。这个过程的描述需要每个网格中心点的高程、在速度点上的水位以及区分“干”、“湿”地区的水位或速度界线这样3个因子。对于二道白河这样的沟谷，地形起伏十分大，用平均水位去计算流体速度会因为沟底地形而造成很大的误差。为了避免这些误差，需要将计算网格划分的足够小，因此计算的步长也相应的要求划分到足够小，至少每1个计算步长的时间内不能跨越一个计算网格。因此，本次研究选取了较小的网格大小(见表 1)及计算步长(计算步长设为0.001 min)。需要指出的是，本次研究所建立的模型是假设z方向上均匀的2D模型。

图 2(Figure 2)

图 2 研究区域及其网格划分示意图

表 1(Table 1)

表 1 模型参数 参量 最大高程/m 最小高程/m 长/km 宽/km 网格精度/m 曼宁系数 库朗数 流通区 2 600 1 250 12 5(最宽处) (0.5, 5)×(0, 0.5) 0.15 4.2 堆积区 1 250 900 16 19(最宽处) (5, 10)×(2, 10) 0.05 3.2

表 1 模型参数

在火山泥石流模拟中，流体性质是最难处理的一个部分，因为火山泥石流爆发时，火山喷发涌出的火山口湖的水一般都混杂着岩浆、火山碎屑等物质，火山顶上冰、雪融化所形成的水也都夹杂着粘土等物质，这使得火山泥石流的成分和性质都相当复杂。盛俭等(2012)对天池火山泥石流堆积物的6个样品进行粒度维度分析，其平均维数为2.564，推断长白山天池火山泥石流主要为稀性泥石流^[14]。本次研究将火山泥石流简化成Newton体，不考虑流体中物质的沉积以及其非Newton体的性质^[15-17]。模拟火山泥石流流量是很难确定的，天池火山口积聚着20.4×10⁸ m³的湖水以及大量的冰雪，这些在火山喷发时都可能成为引发火山泥石流的水源。本次研究分别假设一次中等规模和一次较小规模的喷发造成瀑布口溃决。中等规模的喷发造成的溃决引起火山湖水25%泄流，喷发的高温造成火山顶冰雪融化，共有8×10⁸ m³的水作为火山泥石流的水源，峰值流量达到4 000 m³/s。较小规模的喷发引发火山湖水10%泄流，汇集冰雪融化共有3×10⁸ m³的水源，峰值流量达到2 500 m³/s。利用Delft 3D软件，对上述建立的模型进行数值模拟，得到了这2次不同规模喷发所引发的火山泥石流在距离火山口不同距离处的平均流速、峰值流速(见图 3)。

图 3(Figure 3)

图 3 距离火山口不同距离处火山泥石流平均速度和峰值速度变化曲线

从图 3可以看出，中等规模和较小规模喷发所引起的火山泥石流的流通平均速度和峰值速度随着距离火山口的远近变化较为明显，在火山口水流流出的地方速度较大，然后由于重力等的影响，速度逐渐加大，到距离火山口5~8 km左右平均速度和峰值速度都达到了最大值，到15 km左右速度开始骤降，然后是一个相对平稳的缓慢下降过程。从模拟的结果看，喷发规模的大小只是在幅值上对流通速度起到了影响，但对速度变化规律的影响不大。

中等喷发规模和较小喷发规模所引发的火山泥石流的平均速度和峰值速度所模拟出来的变化规律是相似的，只是在幅值上略有变化。因此本文只对中等喷发规模所引发的火山泥石流进行进一步计算，得到火山泥石流不同时间内的影响范围(见图 4)以及火山泥石流流通途径中不同时间段各剖面上泥石流厚度(见图 5)等。

图 4(Figure 4)

图 4 不同时间火山泥石流的影响范围

图 5(Figure 5)

图 5 各剖面处不同时间火山泥石流的深度

从图 4可以看出火山泥石流的流通是需要经过一段相对较长的时间，从火山口流通到下游28 km左右的地方需要经历近140 min，合理的利用好这段时间对于预防火山泥石流灾害具有重要意义。从图 5可以看出在距离火山口不同距离的剖面处(剖面位置见图 4)火山泥石流的深度随时间变化的规律，从该图也可以看出火山泥石流不同时段的影响范围。

根据本次研究主要得到以下结论

1) 喷发规模的大小对火山泥石流流通的平均速度和峰值速度的变化规律影响不大，只是在幅值上产生变化。

2) 中等喷发规模引发的火山泥石流从火山口流通到下游28 km左右的区域，也就是二道白河镇附近，需要140 min左右，在制定应急预案时该时间可以作为一个参考数值。

3) 模拟给出了距离火山口不同距离处剖面位置上的火山泥石流随时间的深度变化规律，从该图件可以了解到火山泥石流的影响范围。

本次研究只对研究区内的流通区和堆积区进行模拟计算，实际的流通区和堆积区可能都要大于本次工作所划分的区域，如刘祥等(2000)在研究1 000年前长白山天池火山泥石流堆积时就认为该次喷发所形成的火山泥石流一直堆积到吉林市以远^[18]。虽然本次研究的假设和实际情况存在一定的差异，但对长白山地区火山泥石流的流速、堆积情况等在数量及程度方面上进行的分析计算结果对火山泥石流的防治具有参考价值。