, (1)北京大学学报(自然科学版) 第62卷 第2期 2026年3月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 62, No. 2 (Mar. 2026)
doi: 10.13209/j.0479-8023.2025.028
国家重点研发计划(2022YFF0800301)资助
收稿日期: 2025–02–19;
修回日期: 2025–05–06
摘要 从熔岩流数值模拟的理论基础出发, 介绍随机模型、一维管道流模型、元胞自动机模型、基于深度平均的浅水方程模型和三维模型, 探讨不同模型的建模假设、计算效率以及适用场景。结果表明, 随机模型的计算速度快, 但缺乏时间演化数据; 一维管道流模型可快速评估熔岩流的热损失过程; 元胞自动机模型计算效率高且能提供温度和空间形态数据; 深度平均的浅水方程模型适用于平坦地形的熔岩流模拟; 三维模型提供了最详尽的动力学和热学信息, 但计算成本高。最后, 对熔岩流模拟方法的发展趋势做出展望, 包括进一步改进对熔岩流物理性质(黏度和屈服强度)的约束、更精确地描述熔岩流内部的速度和温度分布以及引入流体力学模拟研究的新方法。
关键词 熔岩流; 数值模拟; 随机模型; 一维管道流模型; 元胞自动机模型; 深度平均的浅水方程模型; 计算流体力学(CFD)
溢流式火山喷发表现为岩浆从火山口持续、平稳地溢出, 所形成的熔岩流是一个包含液态熔体、固态晶体与结壳以及气态挥发组分等多相相互作用的复杂系统[1]。研究熔岩流的流动过程不仅对评估地球上的火山灾害至关重要(例如通过预测流动路径来规避基础设施的破坏[2]), 而且可为火星、金星等类地行星表面熔岩地貌的成因解析提供关键约束。20 世纪 70 年代以来, 熔岩流数值模拟方法已从一维稳态模型逐渐发展成二维、三维的高精度模型, 其核心在于通过量化质量、动量与能量的传输过程, 揭示熔岩流的动力学与热力学演化机制[3]。
熔岩流建模的主要挑战源于其非牛顿流体特性(如屈服强度)[4]、相变过程(冷却导致的结晶与表面结壳形成)以及熔体与环境的动态热交换[5–6]。为在计算效率与物理真实性之间取得平衡, 一些模型通常基于质量、动量和能量守恒方程构建, 并需要引入关键假设。例如, 动量方程需结合宾汉姆流体本构关系来描述屈服强度效应, 能量方程则需考虑辐射、对流和传导等多种散热机制。这些方程的非线性耦合特性使得数值求解过程高度依赖简化的物理假设与高效的计算方法。
常用的熔岩流数值模拟方法包括 5 类。其中, 随机模型(如 DOWNFLOW)通过随机扰动地形高程, 生成多条可能的流动路径, 计算速度快, 适用于快速灾害风险评估, 但缺乏动力学细节刻画[5]; 一维管道流模型(如 FLOWGO)聚焦热损失与流变学演化过程, 将几何形态简化为稳态通道流, 能够高效地预测熔岩流延伸长度, 但忽略了横向扩展与地形影响[7]; 元胞自动机模型(如 MAGFLOW)将空间离散为二维单元, 基于简化的纳维–斯托克斯稳态解更新厚度与温度分布, 兼顾计算效率与空间形态信息, 但难以准确地模拟惯性效应与垂向热分层[8]; 深度平均浅水方程模型通过垂直积分, 将三维问题降维为二维, 适用于平坦地形中的黏性流动动力学模拟, 但目前对速度与温度剖面的假设尚需完善[2–3]; 三维计算流体力学模型(如 LavaSIM)完整地求解三维守恒方程, 能够提供最详尽的流场与温度场信息, 但计算成本极高, 依赖高性能计算资源[6]。
各模型在计算效率与物理真实性之间存在显著差异, 应用的场景也有所不同。随机模型和一维模型适用于应急响应与参数敏感性分析, 元胞自动机模型在中等计算成本下提供二维形态数据, 浅水方程模型适用于中等复杂度地形, 三维模型则专门用于高精度灾害评估与机理研究。
基于对现有理论和模型的分析, 我们认为熔岩流模拟方法未来的发展将围绕两个方面: 首先是物理参数约束的精细化, 完善屈服强度与黏度模型, 尤其是其随温度、剪切速率与晶体含量的动态变化, 需要结合实验与现场观测数据, 构建更普适的本构关系, 如 Herschel-Bulkley 模型[3,9]; 其次是加强对熔岩流内部温度和速度分布耦合的约束, 分析熔岩流内部速度与温度梯度的相互影响机制, 从而提升浅水方程与元胞自动机模型的物理真实性[10]。
熔岩流是一个复杂的多相系统, 包含液态(熔岩的熔融部分)、固态(如晶体和结壳)和气态(所有挥发性组分)[11], 这使得熔岩流的建模存在众多难点。熔岩流侵位过程中发生的相变也带来额外的复杂性。熔岩从火山口喷发后开始冷却, 这个过程中发生晶体的成核和生长以及气体的出溶和损失[7,12], 并出现冷却分层现象[13]。因此, 熔岩流属于开放系统, 会与周围环境发生物质和热量交换, 其物理和流变特征在时间和空间上持续演化[14–15]。同时, 熔岩流的物质供应通常不稳定, 在整个喷发过程中可能出现显著的波动[16–18]。
我们将溢流相喷发产生的熔岩流的流动过程归纳为物质、动量和能量传输 3 个方面, 这与地球表面河流的流动有相似之处, 但也表现出一些特殊的性质, 例如非牛顿流体[19–20]、表面结壳和熔岩管道的形成[13,21]、对下伏地表的侵蚀作用[22]和活跃熔岩流系统(active lava flow system, ALFS)的发育[15]等。
若要准确地建立熔岩流动的数学模型, 通常需要基于质量、动量和能量守恒定律来分析。对于连续流体, 其质量守恒方程可表述为
, (1)其中, ρ是流体密度, u是流体速度矢量。在本研究关注的温度与压力范围内, 熔岩流的密度变化不大, 因此可将熔岩流视为不可压缩流体, 即
。 (2)
动量守恒方程的一般形式为
, (3)其中, f是单位质量流体所受的合外力(一般是重力), p是流体内压力, τ是偏应力张量。为得到偏应力张量与速度的关系, 我们需要讨论熔岩流的力学性质, 得到其本构关系。
对于最简单的黏性流体——牛顿流体, 式(3)可以转化为纳维–斯托克斯方程:
(4)纳维–斯托克斯方程是描述黏性流体动量守恒最常用的方程，大量流体模型均基于纳维–斯托克斯方程的合理简化而建立。
熔岩流的物质传输过程比牛顿流体更加复杂。Shaw 等[23]通过测量夏威夷玄武岩的流变学数据, 首次提出采用屈服强度(yield strength)描述其非牛顿流体特性, 即当岩浆所受剪切力超过某一阈值时, 才会发生剪切流动, 该阈值称为屈服强度。Pinker-ton 等[24]在测量埃特纳火山岩浆的流变学性质时, 也得到相似的实验结果。具有屈服强度的黏性物质模型的最简单描述是宾汉姆流体(kingham fluid)模型。宾汉姆流体是一种非牛顿流体, 其内部受到的剪切力 σ达到某一阈值(即屈服强度)后, 流体才能发生剪切流动, 此后剪应力与剪切速率之间呈线性关系:
, (5)
为流体的剪切速率, η为流体黏度。参考文献[25], 我们可以将该本构关系推广至高维情况。宾汉姆流体的本构关系为
(6)
其中, 
, 是偏应力张量t的第二不变量, 我们可以将其视为偏应力张量的一种“模长”; 
, 是剪切速率。如果只考虑单一方向的剪切, 该本构关系可简化为式(5)的形式。与牛顿流体相比, 可将宾汉姆流体的表观黏度视为在牛顿黏度基础上增加一个与剪切速率和屈服强度相关的附加项
。
约束屈服强度对流体流变学行为的影响研究面临两方面的挑战。一方面, 与较为成熟的岩浆黏度模型[26–29]相比, 屈服强度模型的建立仍然不够完善[30–31]。另一方面, 尽管黏度模型更加完善, 但现有模型多基于固定剪切速率下表观黏度的实验结果(例如 Di Fiore 等[32]采用中心圆筒实验), 均未能充分揭示表观黏度与剪切速率之间的关系。这使得准确地构建流体的本构关系变得困难。
总之, 屈服强度对熔岩流的流变行为具有显著影响, 但其约束条件尚不明确。按照对流体的处理方式, 目前的模拟方法可分为两类。
第一类方法基于适当的假设, 将模型简化为只需求解一维稳态问题。我们考虑熔岩流沿具有一定坡度的斜面运动, 形成具有固定宽度和厚度的稳态层流, 如图 1 所示。假设熔岩流的宽度远大于厚度, 可以忽略侧面边界对速度场的影响, 且流动方向始终沿斜面向下。在此条件下, 速度向量 u 可简化为仅沿 x 方向的分量(可以记为标量
), 并且 u 仅在 z方向存在变化。因此, 在稳态层流假设下, 式(3)中的时间导数项满足
。我们进一步展开对流项
:
熔岩流沿着斜面下降方向(x轴正方向)流动, z轴垂直于斜面, y轴沿着熔岩流宽度方向。假设宽度w远大于熔岩流厚度
图1 稳态层流示意图
Fig. 1 Schematic of steady laminar flow
(7)由于速度只有沿流动方向的分量
不为零, 且仅依赖于 z, 因此代入式(7)可得
(8)
本构关系为
, (9)动量方程为
。 (10)
将本构关系带入动量方程, 得到
。 (11)式(11)是一维稳态问题的纳维–斯托克斯方程, 屈服强度通过该方程的边界条件体现, 通常在低剪应力处(如熔岩流表面)形成不发生剪切运动的区域。
第二类方法直接基于三维的纳维–斯托克斯方程进行分析, 但在判断熔岩流停止条件时引入屈服强度的影响。熔岩流的能量传输过程是复杂的, 包括岩浆与大气的热辐射、熔岩流内部和熔岩流表层大气的对流散热以及熔岩流下底面向地层的传导散热[33–34]。若熔岩流表面形成结壳, 其传热模式也会发生很大的改变[13]。我们可以采用与分析熔岩流的动量类似的思路, 得出如下能量守恒方程:
, (12)其中, U 是单位质量熔岩流的内能, k 是热传导系数, T 是温度, 式(12)描述熔岩流内部的热传导、热对流和摩擦生热。
式(2)、(3)和(12)构成一个描述熔岩流流体动力学和热力学过程的偏微分方程组(partial differential equations, PDE)。但是, 熔岩流是一个高度开放且复杂的系统, 其运动与热演化过程的最重要影响因素是熔岩流与外界的相互作用。因此, 仅仅依靠上述 3 个公式描述的守恒关系, 难以构建具有实用价值的熔岩流模拟方法。例如, 影响熔岩流流动的主要因素包括地形坡度和熔岩体表部的高差[5], 而控制熔岩流热演化的主要因素是熔岩流向上部空气和下部基底的热散失[35]。此外, 熔岩流模型的建立还需考虑其他因素, 如初始地形的变化[36]、喷发速率的约束[37]以及熔岩流的流变学性质(例如酸性与基性岩浆的差异)[27]等。
根据所需解决的问题, 不同的熔岩流模拟方法对熔岩流运动提出合理的假设并建立模型, 再选用适当的数值方法进行求解, 最终从模拟结果中获得答案。相比于实际的物理过程, 模型进行了基于假设的简化。简化程度越高, 所依赖的先验约束(主要指熔岩流的物理性质)越少, 消耗的计算资源也相对越少, 同时计算结果中包含的信息也越少。由于计算能力和时间的限制, 在需要快速预测熔岩流运动的场景中, 需采用充分简化的模型, 例如随机模型 DOWNFLOW[5]和一维模型 FLOWGO[7]。如果需要获得熔岩流活动的详细特征, 则需选用简化程度低, 且能充分描述物理过程的计算流体动力学(computational fluid dynamics, CFD)模型, 例如基于深度平均浅水方程的二维模型[2]和三维模型 Lava-SIM[6]。元胞自动机模型(cellular automata, CA)[8,38]是一种广义的二维建模方法, 相比于基于偏微分方程的二维模型, 能够以较低的计算成本获得较好的模拟效果, 是一种兼具效率与实用性的常用模拟手段。上述模型方法在物理实现、数值精度和计算效率方面存在差异, 表 1 总结了各模型的优缺点。
采用随机过程替代复杂的方程求解, 是一种能够快速获取简单结果的思路。其中, 最简单的一类随机模型将熔岩流的侵位过程简化为重力流[5], 其流动依靠重力驱使, 以数字高程模型(digital eleva-tion model, DEM)作为地形数据, 判断熔岩流流动方向。DOWNFLOW 模型是一种采用这种思路的熔岩流模拟方法[5,11], 类似的模型还有 ELFM[39]。
DOWNFLOW 模型[5,11]基于随机方法对熔岩流路径进行快速预测, 并生成相应的熔岩流流域图。为实现最大程度的简化, 该模型将地形重力视为控制熔岩流运动的关键因素, 未对其时间演化过程进行建模, 并忽略流体内部因素的影响[5]。熔岩流的流动路径无法沿坡度增大方向传播, 最有可能沿坡度最大下降方向推进。DOWNFLOW 模型以 DEM为输入, 生成不规则三角网络(triangular irregular network),用来模拟实际地形。在此基础上, 进一步评估该模拟地形上的最陡下降路径(steepest des-cent path, SDP)。对于一个确定的模拟地形, 只能得到一条最陡下降路径。考虑到熔岩流的宽度、高度和跨越地形障碍的能力, 我们在模拟地形的每个网格点引入幅度为 Δh的随机扰动(图 2(a)), 即微小地调整地形形态, 再重新计算最陡下降路径。每次扰动都生成一条新路径, 该过程循环 N次。将扰动幅度 Δh 作为外部输入参数, 用于表征熔岩流对地形障碍的跨越能力。这一系列生成的路径会覆盖一定的区域(图 2(b)), 在路径分布较为密集的区域, 熔岩流流经的概率较高, 反之则较低。熔岩流路径的停止遵循两个条件: 一是当路径进入盆地后, 会以随机方式传播和扩散, 直到接触盆地边缘而停止; 二是基于历史火山喷发以及熔岩流长度的统计资料, 对熔岩流的流动长度设定上限, 当路径达到最大长度时停止延伸[5]。
表1 几种模拟方法的对比
Table 1 Comparison of several simulation methods
方法优点缺点文献 随机模型计算速度快, 及时生成预期流域图无法给出熔岩流厚度数据和随时间演化的情况[5] 一维管道流计算速度快无法给出熔岩流垂向剖面数据[7] 元胞自动机可得到二维区域的厚度和温度分布无法给出熔岩流垂向剖面数据[8] 深度平均的浅水方程理论更可靠, 可得到二维区域的厚度和温度分布需要预设熔岩流垂向剖面数据[2] 三维模型数据丰富全面运算速度缓慢[6]
DOWNFLOW 模型的优势在于计算效率高且数据量需求低, 单次模拟只需数秒即可完成, 且只需DEM 和两个经验参数(Δh和 N), 无需熔岩温度、黏度等复杂参数。模型输出的概率结果可以直观地展示熔岩流的多路径可能性, 为应急疏散和基础设施规划提供参考。尤其在需要大量模拟的情景下, 其快速生成结果的优势更加突出[11]。然而, DOWN-FLOW 模型也具有明显的局限性。首先, 模型完全缺失动力学过程(如流速、厚度和温度演化等), 无法反映熔岩流的时间演化规律; 第二, 参数 Δh的取值依赖经验, 并且未在平坦地形(无坡度驱动)中测试过; 第三, 忽略了冷却过程对流动路径的反馈, 可能错误地估计实际流动范围。
Favalli 等[5]报道了 DOWNFLOW 模型在 2002和 2004 年埃特纳火山喷发事件中的应用, 2004 年 9月 16 日观测到的真实熔岩流场与模型预测结果高度吻合。该模型基于通过激光扫描获得的 Valle del Bove 山谷上部高分辨率 DEM 数据建立。在 66%的区域内, 模型预测的熔岩厚度误差小于 40cm[40]。
(a)在原始地形上施加±Δh的扰动, 得到多条新的高程剖面; (b)在三角网格上生成的一系列最陡下降路径集合; (c)模型中的三角网格和其中一条最陡下降路径示意图
图2 DOWNFLOW 模型的模拟示意图(修改自文献[11])
Fig. 2 Schematic of the simulation of DOWNFLOW model (modified from Ref. [11])
提升 DOWNFLOW 功能的改进方向包括动态地校准 Δh参数(关联其与熔岩流温度或黏度变化)、与热力学模型(例如 FLOWGO)耦合[41]以及利用机器学习技术优化参数经验依赖。尽管 DOWNFLOW更适用于初步风险评估, 但其在快速灾害响应中的作用不可替代, 常作为一维 FLOWGO 或三维 Lava-SIM 等复杂模型的辅助工具, 共同支撑火山灾害的全链条管理。
熔岩流的降温过程是决定熔岩流活动寿命的最主要过程。因此, 我们对熔岩流的能量转移过程进行建模, 主要考虑熔岩流的热力学过程, 尽可能简化熔岩流的几何模型和的动力学过程。Harris 等[7]提出一种基于热力学–流变学耦合的一维熔岩流模拟模型 FLOWGO, 通过简化流体动力学方程, 揭示熔岩流冷却、结晶和停止前进的物理机制, 分析受限通道中熔岩流的热演化过程及其对流动距离的影响。
FLOWGO 模型假设熔岩流在具有一定坡度的扁平笔直的通道中流动, 通道底部和两侧均不发生侵蚀, 顶部是自由表面(图 3)。将通道流中某一位置(滞后于熔岩流前端)的熔岩流单元视为一维空间中的一个控制体单元。由于冷却作用, 熔岩流随延伸距离增加而减速, 当控制体速度降为零时, 模拟终止。模型假设控制体在流动和演化过程中质量保持不变, 其初始速度取决于喷发率和火口处通道的尺寸。
模型假设表层为冷却形成的结壳, 在冷却过程中, 结壳的面积覆盖率和厚度也不断增加, 若与更高温物质接触, 结壳可能发生破裂。熔岩流核心为高温均质流体, 在流动过程中持续冷却, 基底则由温度较低的底部地壳构成(图 3)。模型通过控制体分析 3 种热量传递过程: 经熔岩流表面的辐射和对流散热、通过基底的传导以及内部粘性摩擦生热, 并同时考虑结壳形成对热耗散的影响, 以此约束熔岩流的冷却速率和结晶速率。计算公式如下:
。 (13)在式(13)中, 计算的热量项包括辐射散热(
)、空气对流(
)、降雨损失(
)、底部热扩散损失(
)和内部摩擦生热(
), Er是喷发率, L是结晶潜热, 
是结晶度的温度梯度。FLOWGO 模型假设在熔岩流从喷发到停止流动的过程中, 结晶程度对流速的影响是线性的, 从而可以直接计算出的值。
流动速度可以通过与晶体和温度相关的流变学参数(黏度和屈服强度)计算[31]得出:
图3 FLOWGO 模型考虑的能量传输过程(修改自文献[7])
Fig. 3 Energy transfer process considered by FLOWGO model (modified from Ref. [7])
, (14)其中, θ为地形坡度, d为通道深度, η为岩浆黏度, σcore和 σbase分别为岩浆核心和底部的屈服强度。黏度和屈服强度与温度有关, 可由经验公式[7,27,42]计算得到。
FLOWGO 模型在熔岩流模拟以及火山灾害评估方面表现出显著的能力。基于 Harris 等[7]对 3 处真实案例(包括 1984 年夏威夷冒纳罗亚火山、1997年 5 月夏威夷基拉韦厄火山的 Pu’u’O’o 火山口以及1998 年 10 月意大利埃特纳火山的 SE 火山口的熔岩流)的模拟, 该模型成功地重现熔岩流的关键参数(如流动长度、速度、宽度、冷却速率、结晶度、黏度和屈服强度), 模拟结果与观测数据高度吻合。例如, 在冒纳罗亚案例中, 模型模拟的熔岩流长度为 25～27km, 与观测结果完全一致; 在 Pu’u’O’o 案例中, 模型模拟的熔岩流长度为 2.5～5.7km, 与观测结果 4km 相近。这些成果有力地验证了该模型的准确性与可靠性。
FLOWGO 模型通过分析熔岩流的热流变关系, 能够确定受限于通道的熔岩流延伸长度。作为一种简单高效的模型, 其突出优势在于可快速评估熔岩流的热损失过程, 进而预测熔岩流的寿命和最大流动长度。已有研究者将 DOWNFLOW 和 FLOWGO结合使用[41]。近期, 还有学者利用该模型研究金星通道熔岩流的平均侵位时间和冷却速度[43]。
2018年, Chevrel 等[44]报道了 FLOWGO 的 Py-thon 的实现版 PyFLOWGO。这个工具具备开源与跨平台的特点, 支持用户选择不同的热通量和流变模型来模拟熔岩流的热力学以及流变学行为。Py-FLOWGO 在设计上注重灵活性, 允许用户在不修改代码架构的情况下切换或扩展模型。目前, Py-FLOWGO 已打包为开源 Python 库, 可便捷地导入各类 Python 项目中重复调用。PyFLOWGO 在模拟熔岩流的热力学和流变学性质时表现出良好的准确性和收敛性, 适用于火山灾害评估和熔岩流路径预测等实际场景。PyFLOWGO 的源代码和相关文档发布于 GitHub(https://github.com/pyflowgo/pyflowgo), 用户可通过 Anaconda 或 pip 直接安装使用。
元胞自动机模型是一种描述二维空间内演化过程的模型, 该模型将实际物理空间抽象为二维动态单元网络, 每个单元包含若干状态函数, 用于表征物理空间中各个位置的物理量。状态函数根据特定的演化函数同步更新, 对应物理空间中的演化过程。给定模型中各单元的初始状态, 并根据物理和化学规则设定演化函数, 随着时间推移, 各单元状态将依据该函数持续更新, 从而实现对复杂空间或结构的动态模拟。Crisci 等[45]首次将元胞自动机的方法运用于熔岩流的模拟过程, Dragoni 等[19]和 Mi-yamoto 等[46]先后提出并完善基于纳维–斯托克斯方程的熔岩流运动数值模拟方法。结合上述熔岩流运动数值模拟方法和元胞自动机的方法，Vicari 等[8]提出 MAGLOW 模型, 并将其成功地应用于埃特纳火山喷发的模拟。该模型能够解释质量守恒、与大气的辐射热交换、黏度随温度的变化以及岩浆的非牛顿流体行为, 并具备较快的计算速度, 在地球上一些熔岩流流动预测中的应用比较成熟[8,19,38,45–49]。
我们参考 Vicari 等[8]提出的 MAGLOW 模型, 阐述基于元胞自动机的熔岩流模型的基本原理。为了描述熔岩流的流动状态和温度, 将熔岩流的厚度H和温度 T作为每个单元的两个状态函数, 厚度的演化函数通过纳维–斯托克斯方程的稳态解确定。Vicari 等[8]将单元之间的熔岩传递视为倾斜平面上不可压缩的稳态扁平层流运动。假设流体为不可压缩的稳定扁平层流[19,42], 且只考虑重力作为外力作用, 则式(11)可简化为
, (15)其中, x轴平行于流动方向, z轴垂直于坡面向上; u是沿 x方向的速度, H是流体厚度, a是坡度。式(15)是简化版的纳维–斯托克斯方程, 描述稳定扁平层流在垂直于地面方向的动量守恒。
如前文所述, 熔岩流的重要特征是具有屈服强度, 只有当流体受到的剪切应力超过屈服强度 σy时, 它才开始流动。在稳定扁平层流中, 流体上部的剪应力较小, 因此有一部分不发生剪切变形(如图 4(a)所示), 这一部分称为塞子[19]。塞子的厚度称为临界厚度(critical thickness)(Hc), 即只有当流体厚度大于 Hc时才会流动。临界厚度计算公式如下:
(16)当流体内部高度在
之上时, 不发生剪切变形(¶u/¶z = 0), 并且流体底部的速度为 0。将二者作为边界条件, 对式(15)在 z方向进行两次积分, 得到流体在 z方向的速度剖面。该速度剖面在z方向再次积分, 得到流体的体积流量 q的表达式(与式(14)等价[31]):
(a)流体分为两部分: 塞子和流动部分[46]; (b)临界厚度的计算: Δz 和 Δh 分别是两个单元间基底高程和熔岩流厚度的差值; (c)单元间的热通量计算。辐射热损失是流向环境的表面, 传导热在单元的接触面之间传递, 对流热由从中心单元移动到邻居单元的浅粉红色部分携带
图4 元胞自动机模型截面示意图
Fig. 4 Schematic cross-section of the cellular automaton model
。 (17)熔岩流的热通量考虑了其与大气的辐射散热、与地面的传导散热和自身的热对流(图 4(c))。模型忽略了熔岩流厚度方向上的温度差异, 将每个动态单元内的熔岩视为一个整体。该热模型基于充分混合假设, 即表面的快速冷却效应由整个岩浆层均匀地分担, 从而加快其整体冷却速率[12]。然而, 对于热分层现象较明显的玄武岩(如夏威夷)[13], 这一假设尚存争议。热分层会导致岩浆表面形成冷凝壳, 而内核仍然保持高温并持续流动。例如, 前面介绍的 FLOWGO 模型即考虑了熔岩流表面结壳的形成机制。这种结壳实际上起到保温作用, 可显著地延长熔岩的流动距离[7]。
Vicari 等[8]将 MAGFLOW 模型应用于 2001 年埃特纳火山喷发事件, 结果显示, MAGFLOW 模型在陡坡区域能够较准确地模拟熔岩流的传播时间和路径, 但是在低坡度区域存在高估覆盖面积的问题。此外, 模型未考虑短暂喷发口的形成, 这在喷发后期对熔岩流场的形态产生显著影响。尽管如此, MAGFLOW 模型仍然被证明是一个有效且稳健的工具, 可用于预测熔岩流可能影响的区域, 为火山监测和灾害预警提供支持。
MAGFLOW 模型也存在局限性。首先, 该模型主要关注熔岩流最终沉积的范围和厚度, 将熔岩流的流动简化为就位过程, 忽略流体动力学和运动学特性, 无法准确地描述熔岩流的复杂流动行为[50]。其次, 模型未考虑内部通道的形成——熔岩流成熟过程中的关键现象。流体通道常在喷发几天后出现, 只有部分区域(通道)允许岩浆从中穿过, 到达熔岩流前端, 这一特性对熔岩流的演化有重要影响[51]。此外, 在纳维–斯托克斯方程的简化过程中忽略了惯性项, 这是基于低流速假设而进行的处理, 因此该简化模型无法有效地模拟高流速条件下熔岩的流动行为。同时, 该模型采用斜坡上的扁平稳态层流作为演化函数, 虽然在形式上应用了纳维–斯托克斯方程, 但并未准确地描述熔岩流在实际二维地形上的运动过程。由于熔岩流动并非多个方向一维流动的简单叠加, 这种简化处理可能导致模拟结果与实际流动行为之间的偏差。
因此, 未来有必要对元胞自动机模型开展更广泛的测试, 通过与实际观测数据以及其他模型的结果进行比较, 验证模型的准确性和可靠性。同时, 可考虑基于深度平均的浅水方程模型, 完善演化函数的理论基础, 从而更有效地描述熔岩流在复杂地形上的运动特性, 提高模型的适用性和预测精度。
目前, 已经有基于元胞自动机算法的熔岩流模拟工具发布, 如 MOLASSES[52]。MOLASSES 主要用于评估熔岩流的潜在影响范围和灾害风险, 该工具基于 DEM 和熔岩流的物理参数(如体积和厚度等), 模拟熔岩在地形上的流动路径, 并通过迭代计算生成熔岩流的覆盖范围。MOLASSES 的模拟结果可进一步用于概率分析(例如蒙特卡洛模拟), 从而评估特定区域面临的熔岩流灾害风险。MOLA-SSES 的源代码可以从 https://github.com/geoscience-community-codes/MOLASSES 获取。
深度平均方法是一种使用深度平均流量变量的二维建模方法, 通过浅水近似假设, 有效地避免三维模型计算量过大的问题。该假设要求流动的深度(厚度)远小于流动的水平尺度。较小的深宽比意味着垂直速度远小于水平速度, 可以忽略。深度平均方法已广泛应用于洪水、海啸和泥石流等灾害预测中[2–3]。经典浅水模型基于不可压缩、等温和非黏性流体假设, 但可通过引入黏度和附加输运方程进行扩展。我们基于浅水方程(shallow water equa-tions, SWE)的二维模型, 利用深度平均的速度、温度和水深的守恒方程, 导出一阶偏微分方程组[2]。
经典浅水方程的假设如下: 1)垂直速度分量可以忽略不计; 2)压力仅包含静水压; 3)水平速度分量恒定的垂直剖面(源于无黏性流假设); 4)等温 流体。
当熔岩流被视为与环境存在热交换且具有自由表面的黏性流体时, 经典浅水方程的假设 1 和 2 不再适用。对于黏性流体(如宾汉姆流体), 垂直剪切效应显著, 且通常温度分布不均匀。因此, 在保留前两条假设的基础上, 将速度变量沿深度进行平均, 可得到如下控制方程组:
, (18)
(19)
(20)
其中, h是流体深度; 
是深度平均下的流体速度分量; βij是 Boussinesq 动量系数; γ是摩擦系数, 取决于流体流变学以及基底和流体的性质。对于牛顿流体, 
[2], 对于宾汉姆流体, γ的形式更复杂[3]。
Costa 等[2]假设速度剖面呈抛物线型, 并由此推导出动量方程。相比理想流体的浅水方程, 该方程被修正。由于速度平方的均值不等于平均速度的平方, 对流项中引入一个额外的系数b, 称为 Boussi-nesq 动量系数或形状因子。该系数取决于速度剖面的具体形式, 反映水平流体速度剖面的剪切分布。该系数通常被简化为 1, 但对于半抛物线速度剖面, 其值为 1.2[53]。
我们设定温度为剖面分段的线性函数, 并计算因下表面发生热传导作用引起的地表附近温度的变化。
为深度平均温度, 由此可得到如下能量方程:
(21)
其中, 右侧各项依次代表辐射、对流散热、传导散热和黏性耗散热, 它们依赖于 E, W, H和 K这 4 个半经验校正因子。bTx 和bTy 是是形状因子，其值与速度和温度的垂直分布有关。
基于深度平均浅水方程的熔岩流模型是一类稳健的物理模型, 它在充分刻画物理过程与保持高效之间实现良好的平衡。Costa 等[2]对 1991—1993 年埃特纳火山爆发期间若干熔岩流事件的模拟研究结果表明, 该模型能够半定量地再现真实熔岩流的流动行为, 并与实地观测结果吻合较好。Biagioli等[3]将该模型应用于 2014—2015 年佛戈火山(Pico do Fogo, Cape Verde)喷发事件。该喷发事件发生在地形相对平坦的区域, 且资料记录详尽, 有利于模型验证。他们采用分辨率为 12m 的 DEM, 设定两个可能的喷发口位置, 并引入平均喷发率和记录的喷发温度等参数。模型采用温度依赖的宾汉塑性流变学模型, 并考虑熔岩与环境之间的热交换。通过网格分辨率测试与参数敏感性分析, 发现 40m×40m的网格能够在计算效率与结果准确性之间取得较好的平衡。在合理地调整黏度、屈服应力和温度依赖性参数后, 模拟所得熔岩流扩展范围与实际观测结果表现出良好的一致性。该研究验证了模型在实际火山事件中模拟的准确性, 表明其在熔岩流灾害评估和风险管理中具有较好的适用性。
对于现有基于深度平均的浅水方程的熔岩流模型, 需要对速度和温度剖面进行约束, 但这些约束在目前模型中过于简化。首先, 以 Costa 等[2]的研究为例, 他们没有考虑熔岩流的屈服强度, 其速度剖面仍然采用简单的抛物线形式, 而非具有无剪切塞子层的分段剖面。其次, 由于温度对黏度和屈服强度具有显著的控制作用, 熔岩流内部温度分布不均匀会导致黏度变化, 进而影响速度剖面结构。这表明温度垂直剖面与速度垂直剖面之间存在耦合关系[10], 而现有浅水模型尚未充分考虑这一耦合效应。Ortega-Moya 等[10]分析了温度与速度剖面通过黏度相互耦合的机制, 其结果对完善熔岩流模型具有重要参考价值。
此外, 在忽略垂直速度和静水压力分布的假设下, 浅水模型仍然需进一步改进, 以便完整地捕捉高梯度地形上发展的快速动力学过程。例如, 熔岩流在陡峭地形上可能形成类似雪崩、山体滑坡和泥石流的快速流动现象。目前已有研究对该问题进行探讨, 并提出改进方案[54–56]。
已经有许多方法可用于实现基于深度平均浅水方程的熔岩流模拟。其中, VolcFlow 是一个采用MATLAB 编写的二维欧拉模型。通过深度平均近似求解质量和动量守恒方程, 能够模拟等温条件下的牛顿流体和宾汉姆流体[39]。此外, 也可借助成熟的多物理场有限元商用软件 COMSOL, 求解基于深度平均的浅水方程[57]。
一维模型无法模拟地形起伏导致的不规则流动, 熔岩流被限制在管道中; 二维模型在模拟高速流体(如 Kilauea 火山)时也存在局限, 因其难以考虑惯性项的影响, 且无法刻画相变过程等结构特征。三维熔岩流模拟代码 LavaSIM[6]采用三维分析和自由表面评估, 通过合适的关系式计算熔岩与地面、空气和水体以及熔体与地壳之间的热传递, 完整地描述熔岩流在三维空间中的质量、动量和能量传递过程, 因此适用于模拟各种类型的熔岩流和流动行为, 例如溢流玄武岩和水下熔岩流。
LavaSIM 模型由 Hidaka 等[6]建立, 方法源于核工业中用于严重灾害分析的 SAMPSON 模型[58–59]。LavaSIM 模型考虑熔岩流的惯性、压力、黏度、液位驱动的动量以及浮力等因素, 并基于地面高程和重力方向(地面坡度)描述熔岩在起伏地形上的流动行为。模型将熔岩流视为化学均质、可以凝固并形成结壳的流体, 其液态部分按牛顿流体处理。La-vaSIM 的运行流程包括以下步骤: 确定岩浆喷发口位置(及对应焓值); 计算屈服强度和最小流层厚度; 分析各单元内部对流界面的流动; 计算对流界面的进出流量; 求解热通量; 重新计算岩浆质量和焓, 更新熔岩流分布状态; 计算固结分数和结壳分布。
LavaSIM 模型将空间分割为若干正方形单元(纵截面如图 5 所示), 包含熔岩流的单元按其状态分为对流单元和自由表面单元两类。对流单元完全充满熔岩流物质(可为固态或液态); 自由表面单元则未被熔岩流物质充满, 可位于地表、熔岩顶面和结壳下方等位置。针对内部(对流单元)和外部(自由表面单元)采用不同的方法求解, 对流单元采用简化标记单元法(simplified marker and cell, SMAC)[60], 自由表面单元采用液面高度函数法(height function method)[61]。两种方法相结合, 满足 SMAC 方法约束自由液面性质的处理需求。
熔岩流在以下两种情况下会停止流动: 1)单元内熔岩流完全凝固, 模型根据能量守恒方程计算单元焓值, 判断单元是否已全部凝固; 2)熔岩流在某个平面位置的厚度低于临界厚度 Hc, 该临界厚度可通过屈服强度理论公式(式(16))计算, 也可基于具体火山的观测数据设定为经验常数。
LavaSIM 可对熔岩流的能量传递过程进行细致的建模, 包括熔岩内部的热量传递和熔岩流向外界的热量散失, 单元之间的热传递方案如表 2 所示。
LavaSIM 将能量传递的求解过程分为以下 3 步(图 6)。1)在求解能量守恒方程之前, 使用一般的传热相关性或热导率方程计算液体、固体、大气流体和地面界面处的热通量。当漂浮在液态熔岩上的结壳下方存在自由表面单元时, 计算结壳与自由表面单元下方液态熔岩之间的热通量。此步骤的特点在于, 除非自由表面单元位于顶部, 否则其热通量不纳入计算。2)结合液态熔岩动量计算的结果求解能量传输问题, 并计算固体之间的热传导。此步骤中, 传热计算获得的热通量将作为边界条件。3)在重新设置单元属性后, 计算对流单元(完整液态熔岩单元或完整结壳单元)与自由表面单元之间的能量传输(与地面直接接触的自由表面单元除外)。该步骤假定对流单元中的熔岩与自由表面单元之间瞬间达到热平衡。
图5 LavaSIM 中各类单元的示意图(修改自文献[6])
Fig. 5 Schematic of the various types of units in LavaSIM (modified from Ref. [6])
表2 LavaSIM 中单元之间的热传递方案
Table 2 Heat transfer scheme between cells in LavaSIM
受热单元传热单元液态熔岩熔岩结壳地表结构 液态熔岩熔岩传输和热传导自然对流自然对流 熔岩结壳自然对流热传导热传导 地表结构自然对流热传导热传导 地面自然对流热传导热传导 大气自然对流和热辐射自然对流和热辐射自然对流和热辐射
为了验证 LavaSIM 模型的有效性, Hidaka 等[6]尝试通过模拟, 重现 1986 年伊豆大岛熔岩流 LBIII(席状熔岩流)和 LCI(线性熔岩流)的观测结果。模拟结果显示, 淹没面积与 LBIII 实际观测数据定性地一致, 流动剖面和厚度分布则与 LCI 的实际熔岩特征定性地一致。这表明该模型能够适用于不同坡度(缓坡或陡坡)的熔岩流动模拟。
图6 LavaSIM 中处理能量传递的分步示意图(修改自文献[59])
Fig. 6 Schematic of the step-by-step energy transfer process in LavaSIM (modified from Ref. [59])
相较于其他熔岩流模拟方法, 三维模型的物理学基础更可靠, 并能提供更详细的熔岩流动力学和热力学信息。然而, 三维模型也存在计算成本较高的问题。此外, 由于模型依赖于多个描述熔岩流的物理参数(密度、比热容、热导率、黏度和屈服强度等), 对这些参数的准确约束直接影响模型模拟的可靠性。
各学科对流体三维模拟的需求日益增长, 现有三维模拟软件多为成熟的商用工具, 可用于模拟三维熔岩流等复杂流体现象。其中, OpenFOAM[62]是基于有限体积法的数值建模工具箱, 内置求解复杂流体、化学反应、湍流、传热、固体力学和电磁学等问题的多种求解器。此外, 用户可根据需要添加新的方程、求解器和应用程序。该软件基于 Open-MPI 实现完全并行化计算, 并具有可集成离散元法(DEM)的外部网格划分接口以及相应的预处理和后处理工具。可访问 OpenFOAM 官方网站(http://open foam.org)了解详细信息。
目前成熟的熔岩流模拟方法基本上可归结为上述方法的一种近似或改进。总的来说, 上述方法涵盖从最简单快捷到最复杂详尽的模拟方式, 不同方法在计算效率和应用场景上各有侧重。例如, 以DOWNFLOW 为代表的随机模型和以 FLOWGO 为代表的一维模型计算效率极高, 能够快速输出模拟结果。采用此类模型开展研究工作时, 模型计算所需时间相较于整体研究周期可以忽略不计。DO-WNFLOW 模型常用于火山喷发期间快速预判熔岩流可能淹没的危险区域[11]。而 FLOWGO 等一维模型由于能够较为细致地刻画温度演化过程, 适用于关注热流变行为的模拟研究, 以便揭示温度演化对熔岩流运动的约束机制[63]。
元胞自动机模型是受研究者青睐的熔岩流模拟工具, 这类模型支持在 GPU 平台并行运行[64], 利用普通个人显卡即可在数小时内完成计算, 并可输出熔岩流的温度分布和空间形态等数据。因此, 元胞自动机模型的效果可与基于深度平均浅水方程的二维模型乃至三维 CFD 模型媲美。就计算速度而言, Dietterich 等[57]的基准测试结果表明, 元胞自动机模型 MOLASSES 比二维模型 VolcFLOW 和 COM-SOL 快 2～3 个数量级, 比三维模型 FLOW-3D 和OpenFOAM 快 4～5 个数量级(表 3), 而模拟效果未见显著差异。
相比于三维 CFD 模型, 基于深度平均浅水方程的二维模型具有更高的计算效率, 且在描述物理过程方面比元胞自动机模型更准确。因此, 它适用于模拟无复杂地形障碍的熔岩流, 或用于地形相对平坦且障碍物较少的场景。对于地形复杂且需准确考虑动力学与热力学演化过程的熔岩流模拟(例如评估熔岩流对植被和建筑物的热影响, 或设计熔岩流干预措施), 可采用三维 CFD 模型进行更详尽的计算, 以便获得更精细更准确的结果。此类模型复杂度较高, 通常需借助高性能计算机进行大规模并行计算。
本文介绍几种溢流相火山活动的数值模拟方法, 讨论其理论基础、建模假设、计算思路以及各自的优缺点。不同的数值模拟方法适用于不同的应用场景, 也存在相互联系。受计算资源的限制, 数值模拟不能片面地追求物理图景的精准再现, 而需基于对物理性质和熔岩流运动状态的理解, 对影响熔岩流运动中的多种因素进行合理的取舍。
通过本文的梳理与探讨, 我们认为熔岩流数值模拟方法主要朝两个方面发展。首先, 需加强对熔岩流物理性质(包括热学性质和本构关系等)的约束与改进。一些研究者探讨将 Herschel-Bulkley 本构关系应用于熔岩流模拟中[3,9]。Herschel-Bulkley 本构关系是对 Bingham 本构关系的扩展, Bingham 本构关系是 Herschel-Bulkley 本构关系的一阶特例。其次, 需进一步刻画熔岩流内部速度和温度的分布特征, 尤其是沿垂直方向的变化规律。Ortega-Moya等[10]分析熔岩流内部温度和速度的垂直分布, 并探讨温度分布对速度分布的影响。他们通过引入 3 种更接近实际情况的温度分布模式, 揭示温度对黏度分布的影响, 并计算熔岩流在垂直方向上的速度分布。结果表明, 这种情况下的速度分布与等黏度情况下的速度分布存在显著的差异, 并且直接影响熔岩流的运动行为。这一研究结果说明, 在模拟中考虑温度和速度分布的相互耦合, 对改进熔岩流模拟方法(如基于深度平均的浅水方程方法)有一定的价值。未来可以从上述两个方面持续推进, 发展更准确更实用的熔岩流数值模拟方法。
表3 计算成本对比[57]
Table 3 Computational cost comparison[57]
测试编号CPU小时数OpenFOAMFLOW-3DVolcFlowCOMSOLMOLASSES BM-A (等温, 等黏倾斜流动) 90.45 20.830.3608 1.52.778×10–4 BM-C (冷却, 等黏性轴对称流动)628.02312.971.527812.00.2744 BM-D (冷却, 凝固, 倾斜流动) 20.09114.670.244412.01.944×10–3
参考文献
[1] Sigurdsson H, Houghton B, McNutt S, et al. The encyclopedia of volcanoes. Amsterdam: Elsevier, 2015
[2] Costa A, Macedonio G. Numerical simulation of lava flows based on depth-averaged equation. Geophysical Research Letters, 2005, 32(5): L05304
[3] Biagioli E, Vitturi M D M, Di Benedetto F, et al. Benchmarking a new 2.5 D shallow water model for lava flows. Journal of Volcanology and Geothermal Research, 2023, 444: 107935
[4] Bagdassarov N, Pinkerton H. Transient phenomena in vesicular lava flows based on laboratory experiments with analogue materials. Journal of Volcanology and Geothermal Research, 2004, 132(2/3): 115–36
[5] Favalli M, Pareschi M T, Neri A, et al. Forecasting lava flow paths by a stochastic approach. Geophysical Research Letters, 2005, 32(3): L03305
[6] Hidaka M, Goto A, Umino S, et al. VTFS project: de-velopment of the lava flow simulation code LavaSIM with a model for three-dimensional convection, sprea-ding, and solidification. Geochemistry, Geophysics, Geosystems, 2005, 6(7): Q07008
[7] Harris A J, Rowland S. FLOWGO: a kinematic ther-morheological model for lava flowing in a channel. Bulletin of Volcanology, 2001, 63: 20–44
[8] Vicari A, Alexis H, Del Negro C, et al. Modeling of the 2001 lava flow at Etna volcano by a cellular automata approach. Environmental Modelling & Software, 2007, 22(10): 1465–1471
[9] Tang H S, Kalyon D M. Estimation of the parameters of Herschel-Bulkley fluid under wall slip using a combination of capillary and squeeze flow viscome-ters. Rheologica Acta, 2004, 43: 80–88
[10] Ortega-Moya J, Martínez-Aranda S, Fernández-Pato J, et al. A vertically non-uniform temperature approach for the friction term computation in depth-averaged viscoplastic lava flows. Journal of Computational Physics, 2024, 519: 113378
[11] Tarquini S, Favalli M. Simulating the area covered by lava flows using the DOWNFLOW code // Harris A J L, De Groeve T, Garel F, et al. Detecting, Modelling and Responding to Effusive Eruptions. London: Geolo-gical Society, 2015, 293–312
[12] Cashman K V, Thornber C, Kauahikaua J P. Cooling and crystallization of lava in open channels, and the transition of Pāhoehoe lava to ’A’ā. Bulletin of Volca-nology, 1999, 61: 306–323
[13] Ball M, Pinkerton H, Harris A. Surface cooling, ad-vection and the development of different surface tex-tures on active lavas on Kilauea, Hawai’i. Journal of Volcanology and Geothermal Research, 2008, 173(1/ 2): 148–156
[14] Robert B, Harris A, Gurioli L, et al. Textural and rheo-logical evolution of basalt flowing down a lava chan-nel. Bulletin of Volcanology, 2014, 76: 1–21
[15] Tarquini S. A review of mass and energy flow through a lava flow system: insights provided from a non-equi-librium perspective. Bulletin of Volcanology, 2017, 79(8): 64
[16] Wadge G. The variation of magma discharge during basaltic eruptions. Journal of Volcanology and Geo-thermal Research, 1981, 11: 139–68
[17] Thordarson T, Larsen G. Volcanism in Iceland in histo-rical time: volcano types, eruption styles and eruptive history. Journal of Geodynamics, 2007, 43(1): 118–152
[18] Gudmundsson A. Volcanotectonics: understanding the structure, deformation and dynamics of volcanoes. Cambridge: Cambridge University Press, 2020
[19] Dragoni M, Bonafede M, Boschi E. Downslope flow models of a Bingham liquid: implications for lava flows. Journal of Volcanology and Geothermal Re-search, 1986, 30(3/4): 305–325
[20] Hulme G. The interpretation of lava flow morphology. Geophysical Journal International, 1974, 39(2): 361–383
[21] Harris A, Rowland S. Effusion rate controls on lava flow length and the role of heat loss: a review. Studies in Volcanology: The Legacy of George Walker, 2009, 2(2): 33–51
[22] Hurwitz D M, Head J W, Wilson L, et al. Origin of lunar sinuous rilles: modeling effects of gravity, sur-face slope, and lava composition on erosion rates during the formation of Rima Prinz. Journal of Geo-physical Research: Planets, 2012, 117(E12): E00H12
[23] Shaw H, Wright T, Peck D, et al. The viscosity of basaltic magma: an analysis of field measurements in Makaopuhi lava lake, Hawaii. American Journal of Science, 1968, 266(4): 225–264
[24] Pinkerton H, Sparks R. Field measurements of the rheology of lava. Nature, 1978, 276(5686): 383–385
[25] Jahromi H A, Raisi A, Ghasemi B, et al. Numerical study of the Bingham fluid flow in a cylindrical en-closure with exact Bingham model. Journal of the Bra-zilian Society of Mechanical Sciences and Enginee-ring, 2020, 42(4): 162
[26] Hui H, Zhang Y. Toward a general viscosity equation for natural anhydrous and hydrous silicate melts. Geo-chimica et Cosmochimica Acta, 2007, 71(2): 403–416
[27] Giordano D, Russell J K, Dingwell D B. Viscosity of magmatic liquids: a model. Earth and Planetary Scien-ce Letters, 2008, 271: 123–134
[28] Mauro J C, Yue Y, Ellison A J, et al. Viscosity of glass-forming liquids. Proceedings of the National Academy of Sciences, 2009, 106(47): 19780–19784
[29] Langhammer D, Steinle-Neumann G, Di Genova D. Modelling viscosity of volcanic melts with artificial neural networks. Authorea Preprints, 2022, 23(12), e2022GC010673
[30] Chevrel M, Platz T, Hauber E, et al. Lava flow rheo-logy: a comparison of morphological and petrological methods. Earth and Planetary Science Letters, 2013, 384: 109–120
[31] Moore H J. Preliminary estimates of the rheological properties of 1984 Mauna Loa lava. US Geol Surv Prof Pap, 1987, 1350: 1569–1588
[32] Di Fiore F, Vona A, Scarani A, et al. Experimental constraints on the rheology of lavas from 2021 Cumbre Vieja eruption (La Palma, Spain). Geophysical Re-search Letters, 2023, 50(4): e2022GL100970
[33] Griffiths R W, Fink J H. The morphology of lava flows in planetary environments: predictions from analog experiments. Journal of Geophysical Research: Solid Earth, 1992, 97(B13): 19739–19748
[34] Griffiths R W. The dynamics of lava flows. Annual Review of Fluid Mechanics, 2000, 32(1): 477–518
[35] 支羽萧, 石耀霖. 月球和地球上玄武岩浆冷却速率差异的数值模拟. 地球物理学报, 2025, 68(2): 595–608
[36] Flynn I T, Chevrel M O, Crown D A, et al. The effects of digital elevation model resolution on the PyFLOW-GO thermorheological lava flow model. Environmen-tal Modelling & Software, 2023, 167: 105768
[37] Harris A J, Dehn J, Calvari S. Lava effusion rate definition and measurement: a review. Bulletin of Volcanology, 2007, 70: 1–22
[38] Bilotta E, Pantano P. Artificial micro-worlds part II: cellular automata growth dynamics. International Jour-nal of Bifurcation and Chaos, 2011, 21(3): 619–645
[39] Damiani M L, Gropelli G, Norini G, et al. A lava flow simulation model for the development of volcanic hazard maps for Mount Etna (Italy). Computers & Geo-sciences, 2006, 32(4): 512–526
[40] Mazzarini F, Pareschi M T, Favalli M, et al. Morpho-logy of basaltic lava channels during the Mt. Etna September 2004 eruption from airborne laser altimeter data. Geophysical Research Letters, 2005, 32(4): L04305
[41] Harris A, Chevrel O, Coppola D, et al. Validation of an integrated satellite-data-driven response to an effusive crisis: the April–May 2018 eruption of Piton de la Fournaise. Annals of Geophysics, 2019, 62(2): VO230
[42] Dragoni M. A dynamical model of lava flows cooling by radiation. Bulletin of Volcanology, 1989, 51: 88–95
[43] Flynn I T, Chevrel M O, Ramsey M S. Adaptation of a thermorheological lava flow model for Venus condi-tions. Journal of Geophysical Research: Planets, 2023, 128(7): e2022JE007710
[44] Chevrel M O, Labroquère J, Harris A J, et al. PyFLO-WGO: an open-source platform for simulation of chan-nelized lava thermo-rheological properties. Computers & Geosciences, 2018, 111: 167–180
[45] Crisci G M, Gregorio D, Pindaro O, et al. Lava flow simulation by a discrete cellular model: first imple-mentation. International Journal of Modelling and Simulation, 1986, 6(4): 137–140
[46] Miyamoto H, Sasaki S. Simulating lava flows by an improved cellular automata method. Computers & Geosciences, 1997, 23(3): 283–292
[47] Barca D, Crisci G, Di Gregorio S, et al. Cellular auto-mata methods for modelling lava flows: simulation of the 1986–1987 eruption, Mount Etna, Sicily. Active Lavas. Routledge, 2022: 291–309
[48] Tarquini S, Favalli M. A physics-based model for the dynamics of lava-flow emplacement // Harris A J L, De Groeve T, Garel F, et al. Detecting, Modelling and Responding to Effusive Eruptions. London: Geological Society, 2015, 293-312
[49] Ishihara K, Iguchi M, Kamo K. Numerical simulation of lava flows on some volcanoes in Japan. Lava flows and domes: emplacement mechanisms and hazard implications. Berlin Springer, 1990: 174–207
[50] Hyman D M, Dietterich H R, Patrick M R. Toward next‐generation lava flow forecasting: development of a fast, physics‐based lava propagation model. Journal of Geophysical Research: Solid Earth, 2022, 127(10): e2022JB024998
[51] Wadge G, Young P, McKendrick I. Mapping lava flow hazards using computer simulation. Journal of Geo-physical Research: Solid Earth, 1994, 99(B1): 489–504
[52] Connor L J, Connor C B, Meliksetian K, et al. Pro-babilistic approach to modeling lava flow inundation: a lava flow hazard assessment for a nuclear facility in Armenia. Journal of Applied Volcanology, 2012, 1 (1): 3
[53] Biagioli E, Vitturi M D M, Di Benedetto F. Modified shallow water model for viscous fluids and positivity preserving numerical approximation. Applied Mathe-matical Modelling, 2021, 94: 482–505
[54] Hergarten S, Robl J. Modelling rapid mass movements using the shallow water equations in Cartesian coordi-nates. Natural Hazards and Earth System Sciences, 2015, 15(3): 671–685
[55] Xia Xilin, Liang Qiuhua. A new efficient implicit scheme for discretising the stiff friction terms in the shallow water equations. Advances in Water Resour-ces, 2018, 117: 87–97
[56] Bachini E, Putti M. Geometrically intrinsic modeling of shallow water flows. ESAIM: Mathematical Model-ling and Numerical Analysis, 2020, 54(6): 2125–2157
[57] Dietterich H R, Lev E, Chen J, et al. Benchmarking computational fluid dynamics models of lava flow simulation for hazard assessment, forecasting, and risk management. Journal of Applied Volcanology, 2017, 6: 1–14
[58] Hidaka M, Ujita H. Verification for flow analysis capa-bility in the model of three-dimensional natural con-vection with simultaneous spreading, melting and soli-dification for the debris coolability analysis module in the severe accident analysis code ‘SAMPSON’ (I). Journal of Nuclear Science and Technology, 2001, 38(9): 745–756
[59] Hidaka M, Sato N, Ujita H. Verification for flow analysis capability in the model of three-dimensional natural convection with simultaneous spreading, mel-ting and solidification for the debris coolability ana-lysis module in the severe accident analysis code ‘SAMPSON’ (II). Journal of Nuclear Science and Technology, 2002, 39(5): 520–530
[60] Amsden A A. The SMAC method: a numerical techni-que for calculating incompressible fluid flows. Los Alamos Sci Lab Rep, 1970
[61] Shimegi N. Suuchi ryuutai rikigaku. Tokyo: Asakura Publ, 1994
[62] Jasak H, Jemcov A, Tukovic Z. OpenFOAM: A C++ library for complex physics simulations // Internatio-nal Workshop on Coupled Methods in Numerical Dy-namics. Dubrovnik, 2007: doi: 10.1186/s13617-017-0061-x
[63] 潘波, 许建东, 林旭东, 等. 基于热流变运动学模型的长白山熔岩流数值模拟. 地球物理学报, 2011, 54(9): 2317–2324
[64] 高超, 孟宪海, 李吉刚, 等. 基于 GPU 的元胞自动机熔岩流动模拟. 计算机工程与设计, 2015, 36(3): 793–796
Numerical Simulation Methods for Effusive Lava Flows
Abstracts Based on the theoretical foundations of lava-flow numerical simulation, we first review five classes of models — stochastic, 1D channel flow, cellular automata, depth-averaged shallow water models, and 3D models —highlighting their underlying assumptions, computational efficiency, and applicable scenarios. Stochastic approaches deliver rapid run-times but supply no temporal evolution; 1D channel models quickly quantify heat loss; cellular automata models balance high efficiency with simultaneous temperature and morphological data; depth averaged models are best suited to low-relief terrains; while 3D models provide the most complete dynamical and thermal information at the price of high computational cost. We conclude by outlining future directions: tighter constraints on lava rheology (viscosity and yield strength), more accurate representations of internal velocity and temperature fields, and the adoption of cutting-edge techniques from broader computational fluid dynamics.
Key words lava flow; numerical simulation; stochastic model; one-dimensional pipe flow model; cellular automata; depth-averaged shallow water equation; computational fluid dynamics (CFD)