引言

长白山天池火山是我国境内保存最为完整的新生代多成因复合火山,位于吉林东部的中国-朝鲜边境上(北纬41°20′-42°40′,东经127°00′-129°00′),处于长白山山脉的最高峰,属于大型的复式火山。全新世以来,天池火山的喷发在时间上初步可分为四期:5000 年前、公元946 年(千年大喷发)、公元1668年和1702年。喷发时均以爆破式火山喷发为主,产生了大面积的空降浮岩堆积、火山碎屑流、火山泥石流和熔岩流(于红梅等,2007;2011;赵波等,2010;赵波,2011;Xu等,2013;Wei等,2013;Yu等,2013;Zhao等,2013)。天池火山的每个活动阶段几乎都有不同规模的火山碎屑流产生,其中以千年大喷发形成的火山碎屑流规模最大,影响范围达2160km2。火山碎屑流具有以下灾害特征:①在成分上主要为火山碎屑物(浮岩碎屑、岩屑和晶屑等)与气体的混合;②流速快(高达200m/s),可攀爬高达1000m的障碍物;③温度高,最高可大于900℃;④搬运距离远,能量大,大规模的可达50—100km。鉴于上述特点使得火山碎屑流具有极强的致灾性。

在国外早期的火山碎屑流灾害研究中,Malin等(1982)开发了Energycone程序用于模拟1980年圣海伦斯火山爆发。为了弥补Energy Cone(energy line)模型中线性流体翻越地形障碍物的不足,Sheridan等(1992;1995)开发了Flow2D;随后Flow3D等模拟模型相继发展起来(Kover,1995)。2002年美国纽约州立大学布法罗分校、国家地理信息和分析中心、纽约州工程设计和工业创新中心参与了一项为期3年的山地物质流研究项目,该项目的研究领域涉及数学模型、地质模拟和地理信息科学,并研制出了Titan2D模型(Sheridan等,2004;余斌,2005)。

实际上在国外的研究中主要分为两步:早期是单纯的火山碎屑流灾害分布研究,如基于GIS的Energycone模型,研究用于人员分散和应急救援;随着研究的不断深入,逐渐变成为基于物理模型的数学研究即科学问题,如Titan2D和PDAC(Todesco等,2002)需要采用快速的数学运算来研究流体复杂的运动过程。而未来的趋势是将上述两方面的研究综合起来,将大量的数据转化为三维可视化图形,使普通人都可以参与火山碎屑流灾害的应急响应之中。

国内火山碎屑流研究起步较晚,研究程度相对薄弱。“九五”期间吉林省地震局在长白山地区灾害区划图的编制与减灾对策的子专题的实施过程中,把长白山地区火山碎屑流分为重灾区和中等灾害区(杨清福,2005)。整个重灾区以天池火山口为中心,半径在40—60km范围。在这个区域内将会发生冲撞、掩埋、高温灼烧等破坏,造成建筑物、桥梁等基础设施毁灭性破坏,在此地区的人员存活几率较小。火山碎屑流中等灾区是指火山碎屑流派生的灰云浪部分,以湍流搬运为主,堆积物具有交错层理,影响范围为75—90km。以上灾害划分均是针对超大规模的火山碎屑流,尽管不够详细,但的确是长白山天池火山碎屑流灾害区划的雏形和基础。

近年来的相关研究成果表明,半经验的Flow3D模型模拟的结果与实际分布吻合较好(Todesco等,2002;Sheridan等,2004;Saucedo等,2005),因此,本文采用了基于Flow3D模型原理的分段网格计算方法,对长白山天池火山未来大喷发可能产生的火山碎屑流进行了灾害区划研究。

1 碎屑流灾害模型
1.1 Flow3D物理模型

Flow3D模型是基于库仑阻力的滑块模型,该模型通过构建一个不规则三角网格或TIN法的3D数字高程模型来计算滑块速度的变化。模型通常是记录块体的轨迹沿着很小的时间增量直到停止,每一时间步长的块体的速度和位置被记录下来(图1)。由于长白山天池火山地区三角网格计算复杂,极易进入死循环,因此,本文对滑块模型进行了修改,在长白山天池火山千年大喷发地质调查的基础上,结合Energycone模型和滑块模型两种模型,对火山碎屑流喷发和流动过程进行优化,把计算路线分为若干段,根据数学模型分别赋予参数、计算移动距离、模拟火山碎屑流流动范围。

图 1 Flow3D滑块模型示意图 Fig. 1Diagrammatic sketch showing movement of sliding particles in pyroclastic flow
1.2 数学模型

Flow3D数学模型主要包含能量守恒方程和阻力计算方程两部分。其中能量守恒方程属性原理如下:

(1)
mgH ρ空气gmH
ρ碎屑
=0.5mv02
mgh+0.5mv02=F1l1+……Fnln+0.5mv02
(2)
(3)
θn=arccos hn
ln

式中,m为碎屑质量,一般看作单位质量;g为重力加速度;v0为初始速度;H为喷发柱高度;h为火口与停止面的垂直高度;ρ空气为空气密度;ρ碎屑为碎屑的密度;θn为某一分段的平均坡度;F为阻力,其中F1、2、3……为某一分段的摩擦阻力;Vn为碎屑某一分段的最终速度。

从以上公式中可以看出,在摩擦力相同的条件下,决定碎屑搬运距离的主要是喷发柱高度。

而阻力计算方程采用的是Mellor的雪崩方程:

τ=a0+a1ν+a2ν2
(4)
a0=μ
(5)

式中,v是速度;a0a1a2分别为基本摩擦力、内摩擦力和湍流的阻力系数。根据物理模型,可将火山碎屑流当作块体,所以湍流可忽略不计,这样湍流的阻力a1a2为零。

在火山碎屑流模拟中,μ是一个非常重要的系数,直接决定着摩擦阻力的大小。天池火山地区的计算路线主要是高山草甸、河谷裸露的岩石和水泥路面等,摩擦系数可根据地物不同而改变,下面重点介绍μ值的选取。

1.3 计算路线及参数的选取

Flow3D模型的一个重要基础就是基于DEM网格的计算,通过生成的tin网格,计算不同节点的能量损耗,确定移动距离。其中μ值是在tin网格两个节点之间垂直方向位移H和实际移动位移L的比值,μ值的选取主要根据地貌的情况确定,具体如表1所示。而本文是通过对DEM进行分段化计算,即简化了Flow3D模型的计算过程。

表1 不同类型岩石地表动摩擦系数μ的取值表 Table 1 Friction coefficients of basements with different rock types

根据长白山的地貌特征和所选路线的情况,在锥体附近(海拔大于1500m)、峡谷底部堆积垮塌的岩石碎块以及气候干冷处,μ值取0.7;而在盾体和台地气候湿润、岩石表面多为潮湿处,μ值取0.5。

Flow3D是一维模型,只能沿着给定的路线进行运动。因此在计算路线选取时主要遵循以下两个原则:①近源、中源一般选取峡谷通道;②远源选取平坦低洼地区。同时选择的依据应与地质调查的认识一致,近源、中源峡谷一般为火山碎屑流主要流动通道,而远源一般在低洼处堆积。根据上述原则,笔者按以下步骤对计算路线进行了选取:

(1)对长白山地区数字地形图进行大地坐标系校正;

(2)在3D analyses中,把长白山数字地形图转化为tin格式;

(3)把tin格式转化为DEM图;

(4)在DEM图上选取计算路线,如图2所示;

图 2 天池火山潜在火山碎屑流路线图 Fig. 2Flow paths of potential pyroclastic flows around Tianchi caldera

(5)根据坡度的不同对计算路线进行分段,结果如图3所示。

图 3 计算路线及分段示意图 Fig. 3Cross-sections and segmentation of flow paths
2 天池火山碎屑流灾害区划图

根据上述所选的计算路线和参数,可计算出在不同喷发规模条件下沿着各线路火山碎屑流的最大流动距离,计算结果如表2所示。

根据地形走势来连接各个喷发规模的流动距离,由此可形成不同喷发柱高度下火山碎屑流的前缘等值线。对于10km高的喷发柱,在火山碎屑流范围内基本上呈圆形分布,并以火口为中心半径约14km,但是在南坡小白山由于受到高山的阻挡,在连接时应顺着山角连接;对于20km高的喷发柱,以火口为中心,约36km为半径分布,火山碎屑流范围基本上呈圆形分布,依旧是南坡高山地貌对火山碎屑流分布影响较大,考虑到火山碎屑流具有较强的爬坡能力,在第二个高山坡脚处受到阻挡,所以这种规模的火山碎屑流分布形态类似于千年大喷发火山碎屑流分布;对于30km高的喷发柱,以火口为中心,约60km为半径分布,考虑到火山碎屑流分布范围广泛并且受地形制约明显,尤其是受到南部和东北部山脉的影响,所以火山碎屑流的前缘连线基本上沿着这些山脉的坡度梯度带分布。

表2 计算结果表 Table 2 Input parameters and calculating results

在完成上述工作的基础上,以天池火山千年大喷发火山碎屑流分布格局为参考并根据天池火山周边地区的地貌特征,笔者编制了天池火山火山碎屑流不同规模喷发的灾害区划图,如图4所示。从图中可以看出:以天池火山为中心,在喷发柱高度为10km的情况下,灾害区划最大半径为13.7km;在喷发柱高度为20km的情况下,灾害区划最大半径为35.4km;在喷发柱高度为30km的情况下,灾害区划最大半径为57.8km。

图 4 长白山天池火山火山碎屑流灾害区划图 Fig. 4Pyroclastic flow hazard map of Changbaishan Tianchi volcano
3 结论与讨论

本文设定10km、20km、30km喷发柱高度来代表天池火山不同规模的爆发,并使用国际上通用的一维滑块模型进行路径计算,同时以地貌为参考把各个路线的前缘点连接起来,从而完成了国内第一幅具有预测含义的火山碎屑流灾害区划图。与天池火山千年大喷发火山碎屑流分布范围图及前人所做的天池火山火山碎屑流历史灾害区划图相比较,本研究中技术路线的选择基本以峡谷和河谷为主,火山碎屑流主流动单元主要沿峡谷流动,远源部分堆积于低洼,这与野外调查结果基本一致。此外,前人研究成果认为天池火山千年大喷发火山喷发柱高度在25—30km之间(魏海泉等,2004),天池火山40km范围内为火山碎屑流重点防护区(刘若新等,1998)。而在本文编制的长白山天池火山火山碎屑流灾害区划图中,在喷发柱高度为20km的情况下,灾害区划最大半径为35.42km,这与前人的结果非常接近,也间接地反映出本次灾害区划模型中所选取的参数具有一定的合理性。

参考文献
[1]刘若新,魏海泉,李继泰,1998.长白山天池火山近代喷发. 北京:科学出版社 1—15.[本文引用:1次]
[2]魏海泉,金伯禄,刘永顺,2004.长白山天池火山地质学研究的若干进展与灾害分析. 岩石矿物学杂志 23(4):305—312.[本文引用:1次]
[3]杨清福,2005.长白山天池火山碎屑流分布图报告. 吉林省地震局 [本文引用:1次]
[4]余斌,2005.美国纽约州立大学布法罗分校火山碎屑流和泥石流数学模型研究近况. 山地学报 23(1):126—128.[本文引用:1次]
[5]于红梅,许建东,赵谊,2007.长白山天池火山千年大喷发空降碎屑物的数值模拟. 地震地质 29(3):522—534.[本文引用:1次]
[6]于红梅,许建东,林传勇,2011.长白山天池火山千年大喷发空降浮岩碎屑的形貌特征分析和最终沉降速度. 地震地质 33(2):440—451.[本文引用:1次]
[7]赵波,许建东,于红梅,2010.长白山地区火山碎屑粒度特征研究. 地震地质 32(2):233—242.[本文引用:1次]
[8]赵波,2011.长白山天池火山千年大喷发火山碎屑流相模式及灾害区划研究. 国际地震动态 389:44—46.[本文引用:1次]
[9]Kover T.P., 1995.Application of a digital terrain model for the modeling of volcanic flows: a tool for volcanic hazard determination.MA thesis, University at Buffalo, Buffalo, New York[本文引用:1次]
[10]Malin M.C., Sheridan M.F., 1982.Computer-assisted mapping of pyroclastic surges.Science217: 637—639.[本文引用:1次]
[11]Saucedo R., Macias J.L., Sheridan M.F., Bursik M.I., Komorowski J.C., 2005.Modeling of pyroclastic flows of Colima Volcano, Mexico: implications for hazard assessment.Journal of Volcanology and Geothermal Research139: 103—115.[本文引用:1次]
[12]Sheridan M.F., Macías J.L., 1992.PC software for 2-dimensional gravity-driven flows: Application to Colima and El Chichón Volcanoes, México.Second International Meeting on Volcanology, Colima,México5.[本文引用:1次]
[13]Sheridan M.F., Macías J.L., 1995.Estimation of risk probability for gravity-driven pyroclastic flows at Volcán Colima, México.Journal of Volcanology and Geothermal Research66: 251—256.[本文引用:1次]
[14]Sheridan M.F., Hubbard B. et al., 2004.Pyroclastic Flow Hazard at Volcán Citlaltépetl.Natural Hazards33: 209—221.[本文引用:2次]
[15]Todesco M., Neri A., Esposti Ongaro T. et al., 2002.Pyroclastic flow hazard assessment at Vesuvius (Italy) by using numerical modeling. I. Large-scale dynamics.Bulletin of Volcanology64: 155—177.[本文引用:2次]
[16]Wei H.Q., Liu G.M., Gill J., 2013.Review of eruptive activity at Tianchi volcano, Changbaishan, northeast China: implications for possible future eruptions.Bulletin of Volcanology75, doi: 10.1007/s00445-013- 0706-5.[本文引用:1次]
[17]Xu J., Pan B., Liu T., Hajdas I., Zhao B., Yu H., Liu R. and Zhao P., 2013.Climatic impact of the millennium eruption of Changbaishan volcano in China: New insights from high-precision radiocarbon wiggle-match dating.Geophys. Res. Lett.40: 54—59. doi: 10.1029/2012GL054246.[本文引用:1次]
[18]Yu H., Xu J., Luan P., Zhao B., Pan B., 2013.Probabilistic assessment of tephra fallout hazard at Changbaishan volcano, Northeast China.Natural Hazards69: 1369—1388, doi: 10.1007/s11069-013-0683-1.[本文引用:1次]
[19]Zhao B., Xu J., Lin C., 2013.Study of Distal Pyroclastic-flow Stratum from Tianchi Volcano in 1215 (±15) Eruption: Pyroclastic-flow Over Water: ACTA GEOLOGICA SINICA (English Edition). 87 (1): 801—809 [本文引用:1次]