溢流式熔岩喷发是火山喷发的一种常见形式，也是造成火山灾害最频繁的一种类型.从火口或裂隙呈液态喷发出的熔岩沿地表向前流动，流动中熔岩温度逐渐降低，黏度也随之增大，流动速度减慢，直至最后停滞凝结.炙热熔岩流在流动过程中，会对沿途及周边植被、建筑物和人员等造成破坏和伤亡^{[1, 2]}.因此，对熔岩流的流动过程研究成为了火山灾害研究的一个重点.

熔岩流模拟是利用计算机技术对熔岩流的流动过程和灾害预测的一种有效方法.目前在熔岩流模拟方面主要存在两种模型，一种是Harris, Rowland等通过对夏威夷等火山熔岩流喷发流动过程的研究，提出的热流变运动学模型^[3]；另一种是基于GIS技术，利用有限差分和有限元方法，并依据地形数据所开展的熔岩流模拟^[4~6]，例如，Negro, Fortuna等所用的Magflow 模型.国外学者分别利用两种模型对强活动性的夏威夷莫纳罗亚火山和意大利的埃特纳火山所喷发的熔岩流进行过模拟，模拟结果与实际熔岩流的流动过程、覆盖范围等都具有良好的吻合性^{[7, 8]}.国内穆大卫等^[9]利用元胞自动机法对五大连池新期喷发的熔岩流进行了模拟试验，试验结果与实际情况存在一定的相似性，说明了模拟方法的可行.但由于采用了虚拟的数字地形，与喷发实际地形存在偏差，其工作有待进一步开展.

熔岩流的流动分别受到自身物理性质和地形两方面的影响.上述两种模型都考虑了熔岩自身物理性质和地形两方面因素，但其各有偏重.热流变运动学模型更多的考虑流动过程中物理性质的变化，如热量的散失，结晶的变化等^[3]；而基于GIS 等技术的模拟方法，偏重考虑地形地势方面的影响，适用于复杂地形条件下的熔岩流模拟^[4~6].长白山气象站期碱流岩流动单元单一，所处地形简单(如图 1).所以本研究拟采用热流变运动学模型，着重对熔岩流流动过程中的物理性质方面开展模拟研究.

图 1

Fig. 1

图 1 气象站期碱流岩分布图与中线高程 (a)Landsat ETM 影像；(b)数字地形图；(c)熔岩流中线高程. Fig. 1 Flow map of the Qixiangzhan period pantellerite lava flow and the altitude of channel center line (a) Landsat ETM image； (b) Digital topographic map； (c) Altitude of channel center line.

长白山天池火山是我国最具潜在喷发危险的活动火山.其由于历史上的多次大规模喷发和潜在危险，引起政府和学者的高度关注，并开展了大量的研究工作，取得了相应的研究成果和相关数据^[10~13].熔岩流广泛分布于长白山火山早期造盾、中期造锥和晚期碎屑流堆积的各阶段中，是长白山地区火山喷发与形成的主要形式.另外，长白山地区的熔岩流喷发也经历了从早期的碱性橄榄玄武岩至粗面岩，到最后的碱流岩的演化过程^{[14, 15]}.气象站期碱流岩是从气象站寄生火山口喷发的一套熔岩流，测年结果显示距今约4000a, 为目前所知长白山火山最新一期的熔岩流^[16~20].喷发的熔岩流受地形影响，向北流动约5.4km, 其在LandsatETM 影像和地形图上清晰可辨(如图 1).

本文拟基于热流变运动学模型，利用对此区域和气象站期碱流岩的研究成果，建立相应的模型并进行模拟计算.通过对模拟结果的分析，研究熔岩流流动过程中的热量散失，以及温度、速度、黏度和结晶等变化情况.

炙热熔岩自火口喷溢出能够流动多远?流动过程中受哪些因素影响与控制?这是模型提出并回答的问题，也是模型建立的主要思路.Harris、Rowland等通过对莫那罗亚、基拉韦厄和埃特纳火山的多次熔岩流喷发与流动过程研究，发现熔岩流的流动过程与河流相似，在两侧较高的堤岸控制下沿低洼地形向前向下流动，但其流动过程中存在热量的散失与增加，这会影响熔岩流的速度、温度、黏度和结晶变化，并最终使其停止流动^[3].基于上述的研究发现，提出并建立了熔岩流的热流变运动学模型，并对流动过程进行了方程解析(如图 2).

图 2

Fig. 2

图 2 热流变运动学模型示意图[3] Fig. 2 Schematic diagram of the kinematic thermo-rheological model[3]

热流变运动学模型是以火口喷溢出的一股熔岩流为研究对象，不考虑后续喷溢出的熔岩对其影响，并假定熔岩在设定渠道中流动，且质量守恒，不发生质量亏损.在热流变运动学模型中，熔岩流速度与温度是控制其向前流动的主要参数.速度v=0或温度低于凝结温度时，则熔岩流停滞凝结，否则会继续向前流动.热流变运动学模型中主要控制方程如下:

速度方程:

(1)

其中d为熔岩厚度，ρ_lava为熔岩密度，g为重力加速度，θ 为地形坡度角，η 为黏度，YS_core/YS_base为熔岩流核部和底部的屈服强度比.根据式(1)中可知，θ、η与YS_core/YS_base是影响速度的三个变参数.坡度角θ随地形变化而变化，其体现了熔岩流自身重力在坡度面上对速度的影响；黏度η 与YS_core/YS_base随着熔岩的温度与结晶的变化而变化，其计算公式如下:

(2)

(3)

其中R、a、b、c为常系数，Φ_total 为晶体总含量，Φ_base为底部晶体含量，Terupt为喷出口温度，T_core为核部温度、T_base为底部温度.黏度是描述流体流动性的重要参数，公式(2)表明黏度随熔岩流晶体含量增加和温度降低逐渐上升，其流动性越来越差，速度也逐渐减慢.同样，核部屈服强度YS_core随晶体含量增加和温度降低逐渐上升，但底部屈服强度YS_base基本不变，则YS_core/YS_base逐渐增大.

熔岩流温度的变化与热量系统的变化相关，热量主要包括热损失与热增加两部分，其系统方程为

(4)

式中各热参数的含义与计算方法如表 1和图 2.

表 1(Table 1)

表 1 熔岩流热量系统参数与计算公式[21~28] Table 1 Parameters and formule used in heat budget terms of the lava flow

表 1 熔岩流热量系统参数与计算公式[21~28] Table 1 Parameters and formule used in heat budget terms of the lava flow

则温度的变化为

(5)

式中E_r 为溢出率，ρ_lava为熔岩密度，L_cryst为结晶潜热，δΦ_micro/δx为结晶变化量.

各火山区喷溢的熔岩流在规模、成分、流经地形和环境上各不相同，其流动过程和特征也各不相同.热流变运动学模型通过对喷发规模、岩浆性质、地形和周边环境等参数的设定来适应不同情况下的熔岩流模拟.热流变运动学模型的参数主要分为三部分:地形参数、熔岩参数和环境参数.以下将分别对各类参数的意义和长白山气象站期碱流岩的相关参数选取进行叙述.

地形参数是描述熔岩流流动路径情况的量，其以坡度角的大小来计算.由于所研究的气象站期熔岩流无喷发前的地形资料，且现研究资料也无法反演得到喷发前的地形，因此我们假定熔岩流厚度为常量，通过GIS技术在高精度数字地形图上提取气象站期碱流岩中线高程数据，并计算10 m 间距内的地形坡度角来作为模拟所用的地形参数(如图 1, 4).气象站期碱流岩位于以粗面岩为主的火山锥体北侧，其地形坡度大，最大坡度角达23.4°，其平均坡度角为8.98°.

熔岩参数是描述熔岩流喷发规模、岩浆成分与流动过程的参数，其主要分为两类，即定参数与变参数.定参数主要用来描述熔岩流的喷发初始参数和岩浆性质，其在模拟过程中值不变，如熔岩的密度，喷发量和近地面温度等.变参数则用来描述流动过程中的熔岩性质变化，是控制模拟过程的主要参数，如流动过程中变化的黏度、核部温度和结晶率等.各参数含义及赋值见表 2.

表 2(Table 2)

表 2 长白山气象站期碱流岩模拟参数及赋值表[13, 19, 20] Table 2 Parameters and values used for modeling the Qixiangzhan period pantellerite lava flow

表 2 长白山气象站期碱流岩模拟参数及赋值表[13, 19, 20] Table 2 Parameters and values used for modeling the Qixiangzhan period pantellerite lava flow

环境参数是描述熔岩流喷发与流动过程中环境对其影响的变量，主要包括研究区的降水、风力和气温情况，这些会对熔岩流的热量散失快慢产生影响.各环境参数的含义和长白山地区环境参数的赋值见表 2.

根据热流变运动学模型及其控制方程，本研究使用C 语言编制了模拟程序，模拟计算以1m 为步长，计算1m 内的速度、黏度、热量系统和温度变化等参数.具体计算流程如图 3.

图 3

Fig. 3

图 3 热流变运动学模型的计算流程[3] Fig. 3 Flow chart of the kinematic thermo-rheological model[3]

首先，程序通过文本的读取进行参数的赋值，读入的参数主要包括各方程的系数、地形坡度角和各参数值.依次计算1m 步长内的黏度、屈服强度、速度和各类热损失与热增加.热系统的计算是整个模拟计算的关键，因其影响着熔岩流的冷却速度.根据热系统的损失和增加值，计算1 m 内温度的变化量，并更新温度，同时计算结晶变化量和流动时间等参数.然后通过对速度v和核部温度T_core的判断来决定是否继续运算，如果符合条件，则继续下一步长的计算，否则结束运算，得到流动距离等各参数，并将其写入文本.

长白山气象站期碱流岩模拟共计算5170步，即熔岩流流动5.17km, 这与气象站期碱流岩实测约5.4km 的流动长度基本吻合.在模拟过程中，计算得到每步长内熔岩流的速度、黏度、热量变化及温度等参数，这些参数反映了熔岩流的流动特征和熔岩的物性变化.以下将分别讨论这些参数反映的熔岩流流动特征和参数间的关系.

熔岩流的流动速度是反映熔岩流成灾水平的主要指标，流动速度慢的熔岩流给人们更多的规避机会.图 4绘制了坡度角、黏度、速度随流动距离的变化曲线.由于气象站期碱流岩为一套高黏度和高斑晶含量的酸性熔岩，流动性相对较差，流动速度较慢，模拟计算的流动速度集中在0.05~0.3m·s^-1之间，速度值较小但符合碱流岩流动特征.另外，速度的变化与地形坡度角存在很大的关联，坡度大时速度明显增大，但随着流动距离的增大，黏度对速度的影响也越来越大，速度受地形影响的波动也逐渐减小.黏度值随热量耗散和距离增大呈对数形式逐渐增大，但其波动受地形影响较小.模拟计算的速度变化情况与野外观测研究的熔岩流速度变化情况吻合^{[29, 30]}.

图 4

Fig. 4

图 4 地形坡度角和模拟的速度与黏度变化 Fig. 4 Down flow variation in slope and modeled velocity and viscosity

热量系统的计算是熔岩流模拟的关键部分，因热量散失多少决定了熔岩流冷却速度，影响着黏度、速度、温度以及最终的流动距离.图 5绘制了各种热量随流动距离的变化曲线，其中热损失以实线绘制，包括了辐射热、传导热、风热损、大气传导热和降水热；热增加以虚线表示，包括微晶结晶热和黏度热.与图 4中的速度曲线对比发现，除黏度热与速度的变化趋势一致外，其余热变化与速度变化相反，说明熔岩流速度大时，快速通过此段距离，其热变化值小.另外，在总体趋势上，除黏度热外，随着流动距离的增大和速度的减慢，各种热变化逐渐增大，并在最后速度近为零时，热变化呈近垂直上升趋势.模拟结果反映的热变化情况，与国外火山观测的熔岩流热变化趋势吻合^[31].

图 5

Fig. 5

图 5 模拟的热系统变化 Fig. 5 Down flow variation in modeled heat budget terms

在热量系统中，辐射热Q_rad值约为10⁶w·m^-1，为最主要热损失，其约占整个热损失的50%.而微晶结晶热Q_cryst 是最主要的热增加，值约为10⁵ ~10⁶W·m^-1，约占热增加的95%以上.另外，风热损Qforce约为10⁵W·m^-1，约占热损失系统的30%，而降水热为量最少的热损失，约为10⁴ W·m^-1，仅占热损失系统的1%.

熔岩流随热系统的变化，温度与晶体也随之变化，图 6与图 7分别给出温度和晶体与流动距离和时间的变化曲线.温度的变化是随总热量ΔH的变化而变化，在其与流动距离和时间的曲线上可以看出，ΔH随距离的变化与速度的变化趋势相反，而在时间变化曲线上与速度变化趋势相一致，且总体上呈现自然对数上升.另外，核部温度T_core随距离的增大逐渐降低，在近停滞的5000 m 附近温度降低速率明显加快，而在与时间的变化关系上呈近线性递减的变化，说明温度的变化主要与时间有关.晶体的结晶变化与时间的关系图上，晶体随时间呈现近线性的增长，说明晶体的生长随时间变化.这些参数的变化与理论和实验研究结果相吻合，验证了模拟的可信性^[24].

图 6

Fig. 6

图 6 模拟的总热量与冷却速率变化 Fig. 6 Down flow variation in modeled total heat and coolingrate

图 7

Fig. 7

图 7 模拟的核部温度与晶体变化 Fig. 7 Down flow variation in modeledcore temperature and volume fraction of microlites

在热流变运动学模型中，地形坡度角和热量是控制熔岩流流动的主要参数，也是自然界中喷发熔岩流流动过程中的主要控制因素.然而，实际的熔岩流的流动过程更为复杂，需考虑的因素更多.例如熔岩流的宽度，在模拟中熔岩流呈线性前进，其宽度受速度控制；而实际中熔岩流的宽度不仅受到速度控制，还会受到两侧地形的影响.另外，模型中热量系统的计算公式基本都依据实验或经验的研究成果，当模拟对象和条件改变时，公式及系数也需相应改变，但如何改变更为完善还值得进一步讨论.因此，模型的完善将成为今后讨论研究的重要部分.

参数的设定是影响模拟结果合理与否的重要因素.热流变运动学模型通过对熔岩和环境参数的设定来适用不同地区的熔岩流模拟.这在结果的分析中也知各参数对控制熔岩流的流动起有重要作用.例如，增大黏度将导致熔岩流的流动速度减慢和热量散失加快，并影响最终的流动距离.前人虽对长白山火山区和气象站期碱流岩开展过相关研究工作，但涉及相关参数赋值的研究工作相对不足.因此，开展深入的野外地质工作和参数的提取也将成为下一步工作的重点.

本文基于热流变运动学模型对长白山气象站期碱流岩进行了模拟，模拟的流动长度与测量长度相近，计算得到的模拟数据与观测和研究资料相吻合，这说明了模拟的合理性.通过模拟结果的分析可知，地形坡度角和热量变化量是影响熔岩流流动过程的主要因素，这为应对熔岩流灾害提供参考或依据.例如，面对熔岩流可以筑堤拦挡，或喷水冷却使其停滞.另外，此模拟方法还可为活动火山区的熔岩流灾害区划提供借鉴，为地方规划提供参考数据.