SPH Simulation of Molten-Fluid Flows with a Plastic Surface Skin: A Lava-Flow-Oriented Model Study

1. Introduction

Lava flows represent complex thermofluid phenomena in which a solidified surface layer forms as a result of surface cooling during flow. Understanding the influence of such a solidified surface layer on fluid flow is a fundamental problem, not only in geophysical contexts but also in a wide range of engineering applications. In mechanical engineering, similar phenomena are exemplified by the formation, deformation, and fracture of solidification shells and thin surface skins in continuous casting. Related surface-solidification phenomena are also observed in welding processes, where thin solidified layers can appear on the surface of molten pools and are associated with weld surface formation [1]. These geophysical and engineering phenomena share essential characteristics.

In lava flows, it is known that two types of solidified surface layers are formed depending on the surface cooling conditions [2]. When the cooling rate is high, the lava surface undergoes supercooling and becomes an amorphous glass. As a result, a thin solidified layer (skin) is formed, which deforms plastically. On the other hand, when the cooling rate is low, crystallization progresses and a rigid solidified layer (crust) is formed.

Based on observations of flowing lava, Hon et al. [3] described the advance of lava accompanied by the skin and the crust. The advance is interpreted as a repeated sequence of the processes outlined below. Initially, a plastically deformable skin with a thickness of approximately 1–2 mm forms on the lava surface. This skin retains the high-temperature lava inside and suppresses its motion. As the flow velocity decreases and cooling progresses further, a rigid crust begins to grow from the skin. When the crust thickness reaches several centimeters, it becomes strong enough to halt the incoming lava. Subsequently, when the crust is locally broken due to an increase in internal pressure, the interior lava flows out again, and the surface returns to a state without a solidified layer.

Furthermore, the rheological properties of lava flows are known to change from Newtonian to non-Newtonian behavior due to the presence of crystals and bubbles within the lava [4]. It has been reported that silicate melts and lavas with a low crystal fraction can be treated as Newtonian fluids over a wide range of temperatures and strain rates [5,6]. On the other hand, lavas with a high crystal fraction exhibit Bingham behavior with a finite yield stress [2]. This change in rheological properties has a significant influence on lava flow dynamics. In addition, lava is also a glass-forming liquid that forms an amorphous glass upon rapid cooling. Its viscosity shows a strong temperature dependence near the glass transition temperature [7].

For the mitigation of damage caused by lava flows, accurate prediction of lava flow behavior is an important social demand. Many previous studies have focused on lava properties and characteristic flow phenomena. These factors must be appropriately incorporated into predictive models. In the context, numerous studies have been conducted on lava flows using a variety of approaches, including field studies [3,5,6,7,8,9,10,11], experimental modeling [4,12,13,14], and numerical simulations [15,16,17,18,19,20,21,22,23,24,25].

Field studies include observations of lava flow behavior and cooling processes [3], as well as analyses of the rheological properties of remelted samples of solidified lava [7,8]. In addition, direct measurements of flowing lava have been conducted to evaluate rheological properties, surface temperature, and flow velocity [9]. However, due to the low frequency of eruptions accompanied by lava flows and the complexity of the phenomena, there are limitations to achieving a comprehensive understanding based only on field observations.

Therefore, analog experiments and numerical simulations play an important role in the study of lava flows. In analog experiments, polyethylene glycol (PEG) is often used to model lava. The materials can reproduce the formation of a rigid crust associated with surface cooling [12,13,14]. Based on these analog experiments, Lyman et al. [13,14] extended the similarity solution for viscosity-dominated gravity currents proposed by Huppert [26]. They incorporated non-Newtonian behavior and introduced the formation of a surface crust as a stopping condition for the flow. On the other hand, these models do not account for the influence of a plastic surface skin on flow behavior.

In numerical simulations, models based on the shallow water equation are widely used in operational forecasting of lava inflow areas [15]. Because these equation treats velocity and temperature distributions as depth-averaged quantities, it is difficult to directly represent phenomena that involve vertical structure, such as viscosity changes associated with surface cooling and the formation of a solidified surface layer.

Phase changes by cooling and deformation of solid–liquid interfaces complicate numerical simulations. For this reason, several previous studies have proposed simplified models for heat transport and phase change. For example, Starodubtsev et al. [18,19] represented solidification behavior indirectly by introducing a time-dependent viscosity model, without solving heat transport. In addition, mesh-based computational methods such as finite volume and finite element methods have the disadvantage of requiring remeshing to accurately treat the deformation of the solid–liquid interface [16,17]. In contrast, smoothed particle hydrodynamics (SPH) which is one of mesh-free particle methods, is well suited for simulating flows with phase change and free surfaces. Despite this advantage, only a few studies have reported lava flow simulations that consider heat transport and solidification using particle methods [20,21,22,24,25].

Hérault et al. [24] conducted lava flow simulations for the actual topography of Mt. Etna, considering heat transport, temperature-dependent viscosity, and non-Newtonian rheology. Through qualitative comparisons with observed volcanic landforms and numerical results, they demonstrated the necessity of modeling lava flows as non-Newtonian fluids with temperature-dependent viscosity. However, this study did not provide a detailed discussion of the specific shapes or formation processes of characteristic lava flow structures, such as lava levees and lava accumulations. Tomita et al. [25] conducted SPH simulations of lava flows using molten metal as an analog material. Under simplified conditions, they successfully reproduced the morphology of characteristic lava flow structures and discussed their geometrical features. Nevertheless, previous studies [20,21,22,24,25] have not sufficiently considered the formation of a surface skin or its influence on flow behavior.

Therefore, the objective of this study is to clarify the influence of the plastic surface skin on lava flow behavior. In this study, two-dimensional SPH simulations are conducted to analyze the flow of a virtual fluid representing lava. The fluid flows down on a horizontal plate while forming a plastically deformable surface skin. The plastic surface skin is represented as a low-temperature, high-viscosity region, through the introduction of a temperature-dependent viscosity.

2. Computational Method

2.1. Governing Equations

In this study, the computational object is the lava flow over a horizontal flat plate with cooling from the surroundings. The thermofluid behavior is described by the continuity equation, the Navier–Stokes equation, and the enthalpy transport equation. The continuity equation for a flow with constant density is given by

$$\nabla ·\overrightarrow{u}=0,$$

(1)

where $\overrightarrow{u}$ is the velocity. For a fluid with constant density, the Navier–Stokes equation in Lagrangian form is written as

$${\displaystyle \frac{D\overrightarrow{u}}{Dt}}=-{\displaystyle \frac{1}{\rho }}\nabla p+{\displaystyle \frac{1}{\rho }}\nabla ·\left\{\mu \left[(\nabla \otimes \overrightarrow{u})+{(\nabla \otimes \overrightarrow{u})}^{\mathrm{T}}\right]\right\}+\overrightarrow{g},$$

(2)

where t is the time, $\rho$ is the density, p is the pressure, $\mu$ is the viscosity, and $\overrightarrow{g}$ denotes the gravitational acceleration. Here, the superscript ^T denotes the transpose of a tensor, and the operator ⊗ represents the tensor product. The Lagrangian form of the enthalpy transport equation is given by

$$\rho {\displaystyle \frac{Dh}{Dt}}=\nabla ·(\kappa \nabla T)-Q,$$

(3)

where h is the specific enthalpy, $\kappa$ is the thermal conductivity, T is the temperature, and Q denotes the heat loss term, which includes radiative heat loss and heat loss to the ambient air. Assuming that the specific enthalpy depends only on temperature and using specific heat at constant pressure $c_p=dh/dT$, Equation (3) can be rewritten as the following temperature evolution equation:

$${\displaystyle \frac{DT}{Dt}}={\displaystyle \frac{1}{\rho c_p}}\nabla ·(\kappa \nabla T)-{\displaystyle \frac{Q}{\rho c_p}}.$$

(4)

2.2. Discretization of Governing Equations in SPH Method

2.2.1. Treatment of Physical Quantities and Differential Operators in SPH

In the SPH method, a physical quantity $\varphi$ at the position $\overrightarrow{r}={\overrightarrow{r}}_a$ is evaluated by taking into account the contributions from the neighboring particles b. Defining the relative position vector from position a to position of particle b as ${\overrightarrow{r}}_{ba}={\overrightarrow{r}}_b-{\overrightarrow{r}}_a$, the $\varphi \left({\overrightarrow{r}}_a\right)$ can be expressed as follows [27]:

$$\varphi \left(\overrightarrow{r}={\overrightarrow{r}}_a\right)=\sum _b{m}_b{\displaystyle \frac{{\varphi }_b}{{\rho }_b}}W\left({\overrightarrow{r}}_{ba},H\right),$$

(5)

where m is the particle mass, and W is the kernel function that represents the spatial distribution of the physical quantity carried by each particle. The kernel value is determined by the distance $|{\overrightarrow{r}}_{ba}|$ between points a and b and the particle diameter H. In the following, the notation $W_{ba}=W({\overrightarrow{r}}_{ba},H)$ is used. In this study, the M4 spline function for a two-dimensional space was used as the kernel function [27].

Furthermore, the most fundamental model for the gradient of a physical quantity is obtained by differentiating the kernel function in Equation (5), and is expressed as follows [27]:

$$\nabla \varphi \left(\overrightarrow{r}={\overrightarrow{r}}_a\right)=\sum _{b\ne a}{m}_b{\displaystyle \frac{{\varphi }_b}{{\rho }_b}}\nabla W_{ba},$$

(6)

where $\nabla W_{ba}$ denotes the gradient of the kernel function, which is given by

$$\nabla W_{ba}=-{\displaystyle \frac{{\overrightarrow{r}}_{ba}}{\left|{\overrightarrow{r}}_{ba}\right|}}{\displaystyle \frac{d{W}_{ba}}{d\left|{\overrightarrow{r}}_{ba}\right|}}.$$

(7)

In addition, the Laplacian of a physical quantity is modeled using the gradient of the kernel function, $\nabla W_{ba}$, as follows [28]:

$${\nabla }^2\varphi \left(\overrightarrow{r}={\overrightarrow{r}}_a\right)=2\sum _{b\ne a}{\displaystyle \frac{m_b}{{\rho }_b}}{\displaystyle \frac{{\varphi }_b-{\varphi }_a}{|{\overrightarrow{r}}_{ba}|}}\left({\displaystyle \frac{{\overrightarrow{r}}_{ba}}{\left|{\overrightarrow{r}}_{ba}\right|}}·\nabla W_{ba}\right)=-2\sum _{b\ne a}{\displaystyle \frac{m_b}{{\rho }_b}}{\displaystyle \frac{{\varphi }_b-{\varphi }_a}{|{\overrightarrow{r}}_{ba}|}}{\displaystyle \frac{d{W}_{ba}}{d|{\overrightarrow{r}}_{ba}|}}.$$

(8)

2.2.2. Discretization of Navier–Stokes Equation

In the SPH method used in this study, the continuity equation does not need to be explicitly solved because mass conservation is satisfied when the particle masses and the total number of particles in the computational domain are constant. Therefore, only the Navier–Stokes equation and the enthalpy transport equation were solved.

By applying the SPH discretization to Equation (2), the following form is obtained:

$${\displaystyle \frac{D{\overrightarrow{u}}_a}{Dt}}=-\sum _{b\ne a}{m}_b\left({\displaystyle \frac{p_a}{{\rho }_a^2}}+{\displaystyle \frac{p_b}{{\rho }_b^2}}\right)\nabla W_{ba}+{\displaystyle \frac{2}{{\rho }_a}}\sum _{b\ne a}{\displaystyle \frac{2{\mu }_a{\mu }_b}{{\mu }_a+{\mu }_b}}{\displaystyle \frac{m_b}{{\rho }_b}}{\displaystyle \frac{{\overrightarrow{u}}_b-{\overrightarrow{u}}_a}{|{\overrightarrow{r}}_{ba}|}}\left({\displaystyle \frac{{\overrightarrow{r}}_{ba}}{\left|{\overrightarrow{r}}_{ba}\right|}}·\nabla W_{ba}\right)+{\overrightarrow{g}}_a.$$

(9)

The first term on the right-hand side of Equation (9) represents the pressure gradient, which is discretized in SPH to satisfy the action–reaction law for each interacting particle pair. This term was replaced by the incompressibility algorithm proposed by Shigeta et al. [29], in order to represent incompressible fluid behavior.

The second term on the right-hand side is the viscosity term. In this simulation, the viscosity is varied in both space and time. In such cases, the Navier–Stokes equation must be discretized without assuming a constant viscosity over the computational domain. However, this leads to more complicated discretization. Therefore, this study adopted the scheme proposed by Morris et al. [30], which can take into account for spatio-temporal variations in viscosity.

In this scheme, the viscosity is initially assumed to be constant, and the term is discretized using the SPH Laplacian model presented in Equation (8),

$${\displaystyle \frac{{\mu }_a}{{\rho }_a}}{\nabla }^2{\overrightarrow{u}}_a={\displaystyle \frac{2\mu }{{\rho }_a}}\sum _{b\ne a}{\displaystyle \frac{m_b}{{\rho }_b}}{\displaystyle \frac{{\overrightarrow{u}}_b-{\overrightarrow{u}}_a}{|{\overrightarrow{r}}_{ba}|}}\left({\displaystyle \frac{{\overrightarrow{r}}_{ba}}{\left|{\overrightarrow{r}}_{ba}\right|}}·\nabla W_{ba}\right).$$

(10)

Subsequently, the viscosity $\mu$ in Equation (10) is redefined as the mean viscosity of the two interacting particles.

$${\displaystyle \frac{1}{{\rho }_a}}\nabla ·\left\{{\mu }_a\left[(\nabla \otimes {\overrightarrow{u}}_a)+{(\nabla \otimes {\overrightarrow{u}}_a)}^{\mathrm{T}}\right]\right\}={\displaystyle \frac{2}{{\rho }_a}}\sum _{b\ne a}{\displaystyle \frac{2{\mu }_a{\mu }_b}{{\mu }_a+{\mu }_b}}{\displaystyle \frac{m_b}{{\rho }_b}}{\displaystyle \frac{{\overrightarrow{u}}_b-{\overrightarrow{u}}_a}{|{\overrightarrow{r}}_{ba}|}}\left({\displaystyle \frac{{\overrightarrow{r}}_{ba}}{\left|{\overrightarrow{r}}_{ba}\right|}}·\nabla W_{ba}\right).$$

(11)

The mean value was evaluated using the harmonic mean, in order to treat the large spatial gradients of viscosity.

2.2.3. Discretization of Enthalpy Transport Equation

By applying the SPH discretization to Equation (4), the following form is obtained:

$${\displaystyle \frac{D{T}_a}{Dt}}={\displaystyle \frac{2\delta }{{\lambda }_a{n}_a{\rho }_a{c}_{p_a}}}\sum _b{\displaystyle \frac{{\kappa }_a+{\kappa }_b}{2}}\left(T_b-T_a\right)W_{ba}-{\displaystyle \frac{Q_a}{{\rho }_a{c}_{p_a}}},$$

(12)

where $\delta$ is the number of spatial dimensions, $\lambda$ is the weighted average of the squared values of distances to neighboring particles within the influence radius, and n denotes the particle number density. The first term on the right-hand side represents heat conduction, which is discretized using the Laplacian model of the Moving Particle Semi-implicit (MPS) method [31]. In the SPH method, the Laplacian is represented as a summation of the gradients of kernel functions, which tends to underestimate its value near the surface. On the other hand, in the MPS formulation, the Laplacian is evaluated by taking a weighted average of the kernel functions as a weighting function, which is expected to reduce the underestimation near the surface. The second term on the right-hand side is the heat loss term. According to the previous study by Keszthelyi et al. [10], heat flux due to atmospheric convection and thermal radiation play a dominant role in the cooling of the lava surface. Therefore, likewise in this study, the heat flux to the ambient air $Q_a^{\mathrm{FLUX}}$ and the radiation on the surface $Q_a^{\mathrm{RAD}}$ are considered as the heat loss rate $Q_a$

$$Q_a=Q_a^{\mathrm{FLUX}}+Q_a^{\mathrm{RAD}}.$$

(13)

These effects were applied only to the surface particles. To identify the surface particles, the method proposed by Ito et al. [32] was used. For a surface particle, neighboring particles are unevenly distributed within the influence radius. As a result, the number-density-based center shifts from the actual particle position toward the interior of the fluid. The magnitude of this displacement is computed as

$$\Delta {\overrightarrow{r}}_a^{\mathrm{wc}}={\displaystyle \frac{{\displaystyle \sum _b{\overrightarrow{r}}_b{W}_{ba}}}{{\displaystyle \sum _b{W}_{ba}}}}-{\overrightarrow{r}}_a.$$

(14)

In this simulation, particles with $|\Delta {\overrightarrow{r}}_a^{\mathrm{wc}}|>0.3H$, were classified as surface particles.

The heat flux to the ambient air is formulated as follows:

$$Q_a^{\mathrm{FLUX}}={\displaystyle \frac{S_a}{V_a}}{\displaystyle \frac{{\kappa }_a{\kappa }_{\mathrm{air}}\left(T_a-T_{\mathrm{air}}\right)}{{\kappa }_{a}d{z}^{\prime }+{\kappa }_{\mathrm{air}}dz}}$$

(15)

where S is the cross-sectional length of a particle, V is the area of the computational particle in the 2D simulation, $d{z}^{\prime }$ denotes the distance from the surface to air phase, and $dz$ is the radius of a computational particle. The subscript "$\mathrm{air}$" represents a physical quantity associated with the ambient air. Figure 1 illustrates the concept of this formulation. This equation is derived by assuming that the heat flux from the particle center (point C in Figure 1) to the particle surface (point B) is equal to the heat flux from the particle surface (point B) to the air phase (point A) , and by formulating these conditions as simultaneous equations for the particle surface temperature.

The heat loss due to radiation is expressed by the Stefan–Boltzmann law as follows:

$$Q_a^{\mathrm{RAD}}={\displaystyle \frac{S_a}{V_a}}\epsilon \sigma \left(T_a^4-T_{\mathrm{air}}^4\right)$$

(16)

where $\epsilon$ is the emissivity and $\sigma$ is the Stefan–Boltzmann constant.

In this simulation, solidification was modeled by increasing the viscosity depending on its temperature. This treatment captures the plastic behavior of the cooling-induced amorphous glassy surface skin. The formulation and material parameters of the temperature-dependent viscosity are described in detail in Section 2.3.1.

2.3. Computational Conditions

Figure 2 shows the computational domain. In this figure, brown and red areas represent the locations of solid and liquid particles, respectively. A rectangular container with outer dimensions of $8.0\times {10}^{-1}\mathrm{m}\left(x\right)\times 3.2\times {10}^{-1}\mathrm{m}\left(y\right)$ and a thickness of $4.0\times {10}^{-3}\mathrm{m}$, was prepared. The origin of the coordinate system was defined at a position $4.0\times {10}^{-3}\mathrm{m}$ in the x-direction and $4.0\times {10}^{-3}\mathrm{m}$ in the y-direction from the lower-left corner of the container. In addition, a horizontal plate was placed in the region $0.8\mathrm{m}\le x\le 4.0\mathrm{m}$ and $-4.0\times {10}^{-3}\mathrm{m}\le y\le 0.0\mathrm{m}$. Then, $129,600$ liquid particles were arranged in a regular lattice within the region $0.0\mathrm{m}\le x\le 5.40\times {10}^{-1}\mathrm{m}$ and $0.0\mathrm{m}\le y\le 2.40\times {10}^{-1}\mathrm{m}$. Figure 2 shows a schematic of the particle configuration at this step. As shown here, the liquid particles are initially accumulated in the lower-left region of the container, which form a dam. Subsequently, a dam-break simulation in which the dam collapsed under gravity acting in the $-y$-direction was conducted as a preliminary computation. When the x-directional velocities of all fluid particles converged to $|u_x|<1.0\times {10}^{-3}\mathrm{m}/\mathrm{s}$, this state was regarded as the initial condition of main simulation. Figure 2 illustrates a schematic of the resulting particle configuration. Finally, in the main computation, the right wall of the container was removed, and the fluid particles flowed onto a horizontal plate, as shown in Figure 2. Inside the container, the fluid was treated as an isothermal fluid at a temperature of $T=1,500\mathrm{K}$. Heat transport was considered only for particles located on the horizontal plate ($x>0.8\mathrm{m}$).

In actual lava flows, the thickness of the surface skin is typically $1.0\times {10}^{-3}$–$2.0\times {10}^{-3}\mathrm{m}$ [3]. If the particle diameter was set sufficiently small to resolve the detailed formation process of such a thin solidified layer, the number of computational particles would become excessively large. Therefore, in this simulation, the effect of skin formation on the overall flow was focused and the particle diameter was set to $1.0\times {10}^{-3}\mathrm{m}$.

Due to this setting, the particle resolution is relatively coarse. Because heat fluxes due to atmospheric convection and radiative cooling are averaged over the particle area (volume in three-dimensional simulations) at the free surface, the resulting temperature decrease of surface particles cannot be resolved accurately. Accordingly, in order to enhance the temperature response of surface particles and to reproduce skin formation within a feasible computational time, the effects of atmospheric convection and radiation at the surface particles were artificially amplified by a factor of ten.

Furthermore, sufficient temporal resolution was ensured under a feasible computational cost. A time step of $\Delta t=1.0\times {10}^{-4}\mathrm{s}$ was employed throughout the simulation.

In addition, the ambient air temperature was kept constant at $T_{\mathrm{air}}=300\mathrm{K}$, and its thermal conductivity was set to ${\kappa }_{\mathrm{air}}=2.63\times {10}^{-2}\mathrm{W}/\left(\mathrm{m}\mathrm{K}\right)$ [33]. The distance from the lava surface to the air phase was set equal to the particle radius ($d{z}^{\prime }=0.5\times {10}^{-3}\mathrm{m}$).

Table 1 summarizes the computational conditions. The solid phase was assumed to be soil, and its physical properties were taken from the literature [34].

2.3.1. Modeling of Lava Properties

Because lava effusion is commonly observed for low-viscosity basaltic magma, the material properties of liquid particles in this simulation were based on data for basaltic lava. Comprehensive surveys of physical properties conducted within a single volcano are limited. Therefore, data reported in previous studies from various volcanoes were collected.

For the lava density, a constant value of $\rho =1,370\mathrm{kg}/{\mathrm{m}}^3$ was adopted. It is based on measurements of flowing basaltic lava sampled at Piton de la Fournaise [9].

For the viscosity of lava, experimentally measured values obtained by remelting solidified basaltic lava samples from Mauna Ulu, Kīlauea volcano, were used [7]. In this study [7], the measured viscosity values were fitted using the Vogel–Fulcher–Tammann (VFT) equation and expressed as a function of temperature as follows:

$${log}_{10}{\mu }_{\mathrm{basalt}}\left(T\right)=-5.08+{\displaystyle \frac{6140.5}{T-558.8}},$$

(17)

where ${\mu }_{\mathrm{basalt}}\left(T\right)$ and T are in $\mathrm{Pa}\mathrm{s}$ and $\mathrm{K}$ units, respectively. Figure 3 shows the temperature dependence of the viscosity. The blue curve in this figure represents the viscosity calculated based on Equation 17.

Since the rheological behavior of lava flows is consistent with that of a Bingham fluid. As described above, a model used in the simulation should possess the expression of yield stress. Accordingly, an apparent viscosity model [35] was introduced, in which the apparent viscosity ${\mu }_{\mathrm{app}}$ is defined as a function of the second invariant of the strain rate tensor ${ϵ}_{\mathrm{I}\mathrm{I}}$ given by,

$${\displaystyle {ϵ}_{\mathrm{I}\mathrm{I}}=\frac{\partial u_x}{\partial x}\frac{\partial u_y}{\partial y}-{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{\partial u_y}{\partial x}+\frac{\partial u_x}{\partial y}}\right)}^2,}$$

(18)

where $\overrightarrow{u}=(u_x,u_y)$. The apparent viscosity ${\mu }_{\mathrm{app}}$ is expressed as a function of ${ϵ}_{\mathrm{I}\mathrm{I}}$ as follows:

$${\mu }_{\mathrm{app}}={\mu }_0+{\tau }_0\left[{\displaystyle \frac{1-\exp \left(-\xi \sqrt{|{ϵ}_{\mathrm{I}\mathrm{I}}|}\right)}{\sqrt{|{ϵ}_{\mathrm{I}\mathrm{I}}|}}}\right],$$

(19)

where ${\mu }_0$ is the reference viscosity, ${\tau }_0$ is the yield stress, and $\xi$ is a fitting parameter. This model can express not only the non-Newtonian behavior, but also temperature dependence of viscosity by treating the reference viscosity as a function of temperature.

However, in SPH simulations of viscous fluids, the following stability condition associated with viscosity must be satisfied [28]:

$$C{\displaystyle \frac{4H^2\rho }{\mu }}\ge \Delta t,$$

(20)

where C denotes a safety factor, and a value of $C=0.25$ was adopted in this simulation. Under the particle diameter and time step conditions listed in Table 1, the range of viscosity is estimated as $\mu \le 13.7\mathrm{Pa}\mathrm{s}$ from Equation 20.

To satisfy this constraint, in this simulation, the apparent viscosity given by Equation 19 was scaled by a factor of $1/100$ as follows:

$${\mu }_{\mathrm{app}}^{\prime }={\displaystyle \frac{{\mu }_0\left(T\right)}{100}}+{\displaystyle \frac{{\tau }_0}{100}}\left[{\displaystyle \frac{1-\exp \left(-\xi \sqrt{|{ϵ}_{\mathrm{I}\mathrm{I}}|}\right)}{\sqrt{|{ϵ}_{\mathrm{I}\mathrm{I}}|}}}\right],$$

(21)

where ${\mu }_{\mathrm{app}}^{\prime }={\mu }_{\mathrm{app}}/100$. However, even after this scaling, the apparent viscosity can still exceed the stability threshold by temperature decrease. Therefore, an upper limit of ${\mu }_{\mathrm{app}}^{\prime }\le 13.7\mathrm{Pa}\mathrm{s}$ was imposed. In addition, large spatial gradients of viscosity between neighboring particles can lead to numerical instability. Therefore, a lower limit was also introduced, and the viscosity was restricted to the range $0.685\mathrm{Pa}\mathrm{s}\le {\mu }_{\mathrm{app}}^{\prime }\le 13.7\mathrm{Pa}\mathrm{s}$.

The reference viscosity ${\mu }_0\left(T\right)$ was represented by the temperature-dependent viscosity of basaltic lava ${\mu }_{\mathrm{basalt}}\left(T\right)$. In Figure 3, the purple curve shows the scaled viscosity, ${\mu }_{\mathrm{basalt}}\left(T\right)/100$. The scaled reference viscosity ${\mu }_0^{\prime }\left(T\right)={\mu }_0\left(T\right)/100$ was also restricted by the same upper and lower limits, $0.685\mathrm{Pa}\mathrm{s}\le {\mu }_0^{\prime }\left(T\right)\le 13.7\mathrm{Pa}\mathrm{s}$. In Figure 3, the red curve denotes the temperature dependent reference viscosity ${\mu }_0^{\prime }\left(T\right)$ applied in the simulations with prescribed upper and lower limits. In this simulation, regions where the temperature dependent viscosity ${\mu }_0^{\prime }\left(T\right)=13.7\mathrm{Pa}\mathrm{s}$ were treated as the surface skin. The corresponding temperature for this condition is $T\le 1,306\mathrm{K}$.

In addition, the yield stress was assumed to be constant at ${\tau }_0=1,000\mathrm{Pa}$ based on values reported in a previous study [24].

For the thermal conductivity and specific heat, temperature-dependent properties of basaltic lava were adopted based on previous studies [10,11].

$$\kappa =0.848+1.1\times {10}^{-3}T$$

(22)

$$c_p=\left\{\begin{array}{cc}1100& \left(1010<T\right)\hfill \\ 1211-{\displaystyle \frac{1.12\times {10}^5}{T}}& \left(1010\ge T\right)\hfill \end{array}\right.$$

(23)

where $\kappa$, T and $c_p$ are in $\mathrm{W}/\left(\mathrm{m}\mathrm{K}\right)$, $\mathrm{K}$ and $\mathrm{J}/\left(\mathrm{kg}\mathrm{K}\right)$ units, respectively.

Based on the formulation of the heat flux to the ambient air described in Equation 15, the following equation corresponds to the heat transfer coefficient $\alpha$:

$$\alpha ={\displaystyle \frac{{\kappa }_{\mathrm{particle}}{\kappa }_{\mathrm{air}}}{{\kappa }_{\mathrm{particle}}d{z}^{\prime }+{\kappa }_{\mathrm{air}}dz}}.$$

(24)

By substituting the computational conditions and the physical properties of lava and air described above into this equation, the heat transfer coefficient is evaluated as $\alpha \approx 50\mathrm{W}/\left({\mathrm{m}}^2\mathrm{K}\right)$. This value is within the range of $70\pm 25\mathrm{W}/\left({\mathrm{m}}^2\mathrm{K}\right)$ experimentally obtained for forced convection in a previous study [10].

The emissivity of lava is generally considered to be $\epsilon =0.9$ [2], and this value was also used in the present study.

Table 2 summarizes the physical properties of basaltic lava used in this study.

In this simulation, four different rheological conditions were considered to evaluate the effects of non-Newtonian behavior and formation of the surface skin on lava flow dynamics. Specifically, four rheological conditions were examined, combining Newtonian or Bingham models with constant or temperature-dependent viscosity (N, TN, B, and TB).

• N: Newtonian fluid with a constant viscosity of ${\mu }_a=0.685\mathrm{Pa}\mathrm{s}$
• TN: Newtonian fluid with temperature-dependent viscosity (${\tau }_0=0\mathrm{Pa}$)
• B: Bingham fluid with a constant reference viscosity of ${\mu }_0^{\prime }=0.685\mathrm{Pa}\mathrm{s}$
• TB: Bingham fluid with temperature-dependent viscosity

3. Results and Discussion

3.1. Validation of Numerical Simulations via Comparison with Similarity Solutions

In this study, numerical simulations were considered for a Newtonian fluid with constant viscosity that is instantaneously released onto a horizontal plate. The simulation results are compared with known theoretical scaling laws formulated for viscous-dominated flows [26], as follows:

$$L\left(t\right)=\beta {\left({\displaystyle \frac{1}{3}}{\displaystyle \frac{g\Delta \rho q^3}{\mu }}\right)}^{1/5}{t}^{1/5}$$

(25)

where L is the flow length, g is the gravitational acceleration, q represents the cross-sectional area of the released fluid. The constant $\beta$ is a dimensionless coefficient determined by the similarity solution [26]. In the case of a two-dimensional flow, $\beta =1.41$. $\Delta \rho$ denotes the density difference between the fluid and the ambient medium. In the present simulations, the ambient medium is not modeled. Therefore, $\Delta \rho$ is taken as $\Delta \rho ={\rho }_{\mathrm{fluid}}-{\rho }_{\mathrm{air}}\approx {\rho }_{\mathrm{fluid}}$. This assumption is reasonable for lava flows, since ${\rho }_{\mathrm{lava}}\gg {\rho }_{\mathrm{air}}$.

Figure 4 compares the simulated flow length of condition N with the theoretical solution proposed by Huppert [26]. To compare the simulation results with the theoretical scaling, it is necessary to identify the time range in which the flow is dominated by viscous effects. In actual physical processes, the early stage flow behavior is influenced by inertial slumping, which precedes the viscous-dominated flow [13]. Therefore, this initial transient stage should be excluded from the quantitative comparison.

According to the scaling relation given by Equation 25, the flow length follows $L\propto t^{1/5}$ in the viscous-dominated flow. When the flow follows the scaling, the quantity $L/t^{1/5}$ becomes constant. Therefore, the compensated quantity $\tilde{L}=L/t^{1/5}$ is introduced and evaluated to diagnose the transition to the viscous-dominated flow. Figure 5 shows the time evolution of the relative change of the compensated quantity, described as $\left|{\tilde{L}}_{n+1}-{\tilde{L}}_n\right|/{\tilde{L}}_n$. As time progresses, the magnitude of the relative change decreases. After approximately $t=5.0\mathrm{s}$, the relative change remains at a low level of order ${10}^{-4}$. Based on Figure 5, $t\ge 5.0\mathrm{s}$ is adopted as the viscous-dominated regime.

For the data within the viscous-dominated regime, the relative root-mean-square error (rRMSE) between the simulation result of condition N and the theoretical solution proposed by Huppert [26] is evaluated. The resulting rRMSE is $5.36\times {10}^{-2}$. This result indicates that the simulation results agree with the theoretical scaling within a relative error of approximately $5\%$.

In the following analysis, the time exponent is fixed to $1/5$, and only the constant $\tilde{L}=\beta {\left[g\Delta \rho q^3/\left(3\mu \right)\right]}^{1/5}$ is treated as a fitting parameter. The value of $\tilde{L}$ is determined by minimizing the rRMSE between the simulation result and the theoretical scaling. As a result, the optimal value is found to be $\tilde{L}=2.53$, while the theoretical value is $\tilde{L}=2.39$. With this optimized prefactor, the rRMSE is reduced to $1.25\times {10}^{-2}$. This result indicates that the temporal scaling predicted by the theory is well reproduced by the simulation.

These results indicate that the numerical simulations show reasonable agreement with the known theoretical scaling laws. This agreement supports the validity of the simulations for the Newtonian condition.

For Bingham fluids, Lyman [13] proposed a theoretical model that describes the final stopping position of a flow front. This model is derived from the static balance between the yield stress and the gravitational driving force. The theoretical stopping length $L_y$ is given by

$$L_y={\left({\displaystyle \frac{9g\Delta \rho q^2}{8{\tau }_0}}\right)}^{1/3}.$$

(26)

In this study, the numerical simulation results of condition B are compared with this theoretical prediction, with a focus on the final flow length at stopping.

Figure 6 shows the temporal evolution of the flow front position obtained from the numerical simulation for condition B. The final stopping position predicted by the theoretical model of Lyman [13] is also shown for comparison.

As shown here, the stopping length measured from the origin in the simulation is approximately $L\approx 2.6\mathrm{m}$. In contrast, the theoretical model predicts a stopping length of $L=2.94\mathrm{m}$. The stopping length obtained from the numerical simulation corresponds to approximately $88\%$ of the theoretical prediction. Lyman [13] reported that, for Bingham fluids, the experimentally observed stopping lengths range from $79\%$ to $93\%$ of the theoretical values. In Figure 6, the thin dashed lines indicate this experimentally observed variability range. The result of simulation is within the variability range reported in previous experimental validations.

The result of simulation is within the experimentally observed variability range indicated by the thin dashed lines in Figure 6.

Based on the comparison described above, these results indicate that the numerical simulations show reasonable agreement with the theoretical model for Bingham fluids. This agreement supports the validity of the simulations in capturing the stopping behavior both qualitatively and quantitatively.

3.2. Effects of Bingham Behavior on Flow Deceleration and Distance

Figure 7 shows the time evolution of the flow front position for each condition. The blue, green, orange, and red lines correspond to the conditions N, TN, B, and TB, respectively.

As shown in Figure 7, a clear difference is observed in the approach to the stopped state between Newtonian fluids (N and TN) and Bingham fluids (B and TB). For conditions N and TN, the flow velocity decreases gradually. In contrast, for conditions B and TB, the flow undergoes a rapid decrease in velocity before reaching the stopped state.

Figure 8 shows the velocity distributions in the x-direction for each condition when the flow front reaches $x=1.5\mathrm{m}$. It should be noted that the aspect ratio is changed for visibility in the figures. A difference is observed in the spatial extent of the region with velocities higher than $1.8\mathrm{m}/\mathrm{s}$ between Newtonian fluids (conditions N and TN) and Bingham fluids (conditions B and TB). On the other hand, no significant difference is found in this high-velocity region between the two Newtonian conditions (N and TN). A similar result is obtained between the two Bingham conditions (B and TB).

Figure 9 shows the velocity distributions in the x-direction for each condition when the flow front reaches $x=2.0\mathrm{m}$. The Bingham fluids (Figure 9 and Figure 9) exhibit lower velocities than the Newtonian fluids (Figure 9 and Figure 9). The region with velocities higher than $1.2\mathrm{m}/\mathrm{s}$ disappears in the Bingham conditions.

Figure 10 shows the velocity distributions in the x-direction for each condition when the flow front reaches $x=3.0\mathrm{m}$ for Newtonian conditions (N and TN) and $x=2.4\mathrm{m}$ for Bingham conditions (B and TB). In addition, Figure 11 shows the velocity distributions at $t=10.83\mathrm{s}$ when the flow under condition N reaches the domain boundary. A comparison between Figure 10 and Figure 11 for condition TB shows that the flow front exhibits almost no further displacement. This observation allows the Bingham flows under condition TB to be regarded as stopped at the final stage. In particular, as shown in Figure 10, the flow under condition TB has already lost most of its velocity at $x=2.4\mathrm{m}$.

Figure 12 shows the apparent viscosity distributions for conditions B and TB when the flow reaches at $x=2.0\mathrm{m}$ and $x=2.4\mathrm{m}$. Here, the apparent viscosity distribution under condition TB at $x=2.4\mathrm{m}$ (Figure 12) is compared with that at $x=2.0\mathrm{m}$ (Figure 12), where the flow is still advancing. This comparison shows that the apparent viscosity increases over the entire flow region during this time. A similar comparison is conducted for condition B. The apparent viscosity distribution at $x=2.4\mathrm{m}$ (Figure 12) is compared with that at $x=2.0\mathrm{m}$ (Figure 12). The increase in the apparent viscosity is found to start from the upstream region.

Figure 12 and Figure 12 are compared with the velocity distributions at the same times (Figure 9 and Figure 10). The upstream region exhibits lower velocities than the downstream region. This observation indicates that the apparent viscosity increases in regions where the flow velocity decreases.

Therefore, a feedback process within the Bingham fluid model is suggested. A decrease in flow velocity reduces the velocity gradient. This reduction leads to an increase in the apparent viscosity. The increased apparent viscosity further decreases the flow velocity. The repetition of this process can cause a rapid decrease in velocity and eventually lead to a stopping of the flow.

Through this feedback process, the simulation captures behavior that is consistent with the characteristics of a Bingham fluid. In actual Bingham fluids, a decrease in the velocity gradient causes the shear stress to fall below the yield stress. As a result, the fluid exhibits solid-like behavior.

In this simulation, the Bingham fluid behavior was confirmed to have an influence on the flow deceleration and on the flowing distance. These results are consistent with the understanding [24] that non-Newtonian fluid models are required for accurate prediction of lava flow behavior.

3.3. Effects of Skin Formation on Flow Deceleration and Distance

As shown in Figure 7, the time evolution of the flow front position indicates that conditions N and TN exhibit similar gradual deceleration trends. However, at the same elapsed time for both conditions, condition TN exhibits a shorter flowing distance than condition N. This difference in flowing distance is found to increase as time progresses.

Similarly, conditions B and TB exhibit similar abrupt deceleration trends. However, at the same elapsed time for both conditions, condition TB exhibits a shorter flowing distance than condition B. As a result, the final stopping distance under condition TB is smaller than that under condition B.

Figure 8, Figure 9, Figure 10 and Figure 11 show the velocity distributions in the x-direction for four conditions. First, the velocity distributions of the two Newtonian fluids under conditions N and TN are compared. When the flow front reaches $x=1.5\mathrm{m}$ (Figure 8 and Figure 8), the velocity distributions under conditions N and TN are similar. When the flow front reaches $x=2.0\mathrm{m}$ (Figure 9 and Figure 9), the spatial extent of the region with velocities between $1.4\mathrm{m}/\mathrm{s}$ and $1.6\mathrm{m}/\mathrm{s}$ is smaller under condition TN than under condition N. When the flow front reaches $x=3.0\mathrm{m}$ (Figure 10 and Figure 10), the difference in velocity between conditions N and TN becomes more pronounced.

Similarly, the velocity distributions of the two Bingham fluids under conditions B and TB are compared. When the flow front reaches $x=1.5\mathrm{m}$ (Figure 8 and Figure 8), the velocity distributions under conditions B and TB are similar, as in the Newtonian conditions. When the flow front reaches $x=2.0\mathrm{m}$ (Figure 9 and Figure 9), the velocity under condition TB is lower overall than that under condition B. When the flow front reaches $x=2.4\mathrm{m}$ (Figure 10 and Figure 10), the flow under condition TB has lost most of its velocity, and the front position shows almost no further change, even at $t=10.83\mathrm{s}$ (Figure 11). On the other hand, at $x=2.4\mathrm{m}$ (Figure 10), the flow under condition B continues to advance. At $t=10.83\mathrm{s}$, the front position under condition B (Figure 11) is located further downstream than that under condition TB (Figure 11).

Figure 13 shows the apparent viscosity distributions at the time when flows reach at $x=2.0\mathrm{m}$ for each condition. Under the temperature-dependent conditions TN (Figure 13) and TB (Figure 13), a high-viscosity region with ${\mu }_{\mathrm{app}}^{\prime }=13.7\mathrm{Pa}\mathrm{s}$ is observed at the flow surface.

Figure 14 shows the temperature distributions at the time when flows reach at $x=2.0\mathrm{m}$ for conditions TN and TB. A comparison between Figure 13 and Figure 14 indicates that this surface region satisfies ${\mu }_{\mathrm{app}}^{\prime }=13.7\mathrm{Pa}\mathrm{s}$ and $T\le 1306\mathrm{K}$. A similar relationship between the apparent viscosity (Figure 13) and temperature (Figure 14) is also observed in the condition TB. This result shows that the surface layer is represented as a low-temperature, high-viscosity skin.

Figure 13 for condition N and Figure 13 for condition TN show that the scaled viscosity ${\mu }_{\mathrm{app}}^{\prime }$ inside the flow is mostly lower than $0.7\mathrm{Pa}\mathrm{s}$. This observation indicates that cooling-induced skin formation occurs only at the flow surface.

Similarly, Figure 13 for condition B and Figure 13 for condition TB show comparable viscosity distributions inside the flow. This result indicates that skin formation due to cooling is confined to the flow surface, even for the Bingham fluid conditions.

On the other hand, Figure 12 and Figure 12 show the apparent viscosity distributions when the flow fronts reach $x=2.4\mathrm{m}$. At this stage, Figure 12 and Figure 12 indicate that the apparent viscosity inside the flow is overall higher under condition TB than under condition B. This result suggests that the velocity reduction process associated with Bingham behavior occurs at a shorter flowing distance under condition TB than under condition B.

From the above comparisons of the velocity distributions and the time evolution of the flow front positions, the effects of the surface skin are examined for both the Newtonian conditions (N and TN) and the Bingham conditions (B and TB). These results indicate that the cooling-induced surface skin suppresses the flow and reduces the flow velocity. This interpretation is consistent with qualitative observations of lava flows. In such flows, a plastically deformable surface skin formed at the flow surface retains the internal lava near the flow front [3].

Under condition TN, the suppression of the flow by the surface skin causes the flow to exhibit longer arrival times at corresponding positions than under condition N. Similarly, under condition TB, the surface skin causes the Bingham related velocity reduction to occur at a shorter flowing distance than under condition B. As a result, this behavior leads to a difference in the final stopping position of the flow.

These results suggest that surface skin formation plays a role in the flow characteristics of lava flows, and that accurate simulation of lava flow behavior requires viscosity models that appropriately account for both non-Newtonian behavior and skin formation.

4. Conclusions

The objective of this study was to clarify the influence of the plastic surface skin on lava flow behavior. To this end, two-dimensional SPH simulations were conducted to analyze the flow of a lava-inspired virtual fluid on a horizontal plate. The surface skin was represented as a low-temperature, high-viscosity region through the introduction of a temperature-dependent viscosity. The flow behavior was evaluated by comparing conditions with and without non-Newtonian (Bingham) behavior and with and without surface skin formation.

The main conclusions of this study are summarized as follows:

1.

The numerical results obtained for the isothermal Newtonian condition showed reasonable agreement with established theoretical solution [26] in terms of the time evolution of the flow. In addition, the results for the isothermal Bingham condition were consistent with theoretical prediction [13] regarding the final stopping position of the flow. These agreements support the validity of the present simulation framework for describing viscous-dominated and yield-stress-controlled flow behavior.

2.

By introducing the apparent viscosity model, it was confirmed that the apparent viscosity increases as the flow velocity decreases. This increase in apparent viscosity further reduces the flow velocity and eventually leads to the stopping of the lava flow. This behavior qualitatively reproduces the characteristics of a Bingham fluid. As a result, the flow velocity and the flowing distance are affected. These results are consistent with the prevailing understanding that non-Newtonian fluid models are required for accurate prediction of lava flowing distances.

3.

By considering the temperature dependence of viscosity, a plastic surface skin was successfully reproduced in the numerical simulations. Comparison with simulations without surface skin formation revealed that the surface skin reduces the flow velocity. This behavior is consistent with qualitative observations [3]. In particular, when Bingham behavior is taken into account, the suppressive effect of the surface skin causes the successive deceleration characteristic in Bingham fluids, as described above, to emerge at a shorter flowing distance. Therefore, these results suggest that accurate prediction of lava flowing distances requires models that incorporate both non-Newtonian behavior and surface skin formation.

These present findings highlight the importance of surface skin formation in lava flow dynamics. Similar physical considerations are also relevant to engineering and other thermofluid systems involving surface solidification. Future work will focus on extending the present model to incorporate realistic thermophysical properties of actual lava and to evaluate the quantitative influence of surface skin formation at natural flow scales. In particular, a semi-implicit treatment [22,36] of the viscous term will be introduced to relax the severe time step restriction caused by high viscosity. Such developments are expected to contribute to more accurate prediction of lava flow behavior.