GEOPHYSICAL RESEARCH LETTERS
Supporting Information for “Some lava flows may not
have been as thick as they appear”
Jonas Katona
1
,
2
, Xiaojing Fu
1
,
3
, Tushar Mittal
1
,
4
, Michael Manga
1
, and
Stephen Self
1
.
1
University of California Berkeley, Berkeley, CA, United States
2
Yale University, New Haven, CT, United States
3
California Institute of Technology, Pasadena, CA, United States
4
Massachusetts Institute of Technology, Cambridge, MA, United States
Contents of this file
Introduction
Texts S1 to S7
Figures S1 to S11
Table S1
Additional Supporting Information (Files uploaded separately)
Captions for Movies S1 to S3
November 28, 2021, 8:23am
X - 2
:
Introduction
In text S1, we describe how we derived or calibrated the phase-field parameters
τ
,
ω
2
φ
,
and
H
from either known quantities or the interfacial width
d
, which numerically acts
somewhat like a viscosity/smoothing term; Figure S1 gives relevant details as to how
we calibrated
d
. In text S2, we describe the mathematical formulation for the initial
and boundary conditions in the phase-field model, as well as further explain some math-
ematical notation. In text S3, we describe the numerical scheme we used to integrate
the phase-field equations given in the main text for the initial and boundary conditions
given in text S2. In text S4, we define and explain the relative
L
2
error, which is used
in Figure S4 and for Milne’s device in our adaptive time-stepping scheme. In text S5,
we describe the choice of grid size ∆
x
for each simulation. In text S6, we explain the
process of columnar joints which appear as a result of cooling, and how different cooling
mechanisms generally result in different columnar jointing patterns, which is relevant for
describing flows in the main text (e.g., the Cohassett flow). Finally, text S7 provides
relevant details for how we compared our model with the empirical data given in Hon et
al. (1994) and explains in greater detail the successes and shortfalls of our model as shown
via this comparison than in the main text.
Figures S1-S3 contain extra figures that we used to compare our phase-field model to
results elsewhere in the literature — namely, Wright and Okamura (1977); Wright et al.
(1976); Wright and Marsh (2016) for Figure S1 and Hon et al. (1994) for Figures S2-S3.
Figures S4-S9 showcase several compelling yet slightly tangential properties we observed
in our simulations of lava emplacement dynamics, as well as how the qualitative aspects
November 28, 2021, 8:23am
:
X - 3
of lava solidification dynamics change with the emplacement time interval
t
emp
and lobe
height
h
. The captions for Figures S5-S9 are contained in the body of the supplement,
while the captions for for all other figures (including for Figures S10 and S11, which are
high-resolution versions of the last two examples provided in Figure 2 of the main text),
have captions attached directly below each figure. Table S1 contains the values of the
parameters we used in our model. Finally, we provide captions for three movies which
correspond to the three cases featured in Figure 1 in the main text.
Text S1. Formulation of parameters
In practice, any phase-field model has an acceptable range of values for its phase-field
modeling parameters which primarily corresponds to how much we want to “smooth out”
the interface between the different phases we are modeling, i.e., how much we want to
reduce the discontinuity. Our model loses precision if we smooth the interface out too
heavily, but our model will also become stiff and therefore difficult to compute if we do not
“smooth out” the interface(s) sufficiently. But if we are less concerned about the size of
the interface and more about the macroscopic properties that a given multiphase system
has, then we generally do not lose much if we smooth out the interface, even considerably
so, as we see in Figure S1.
However, there are still limits to this smoothing, since the interface could become large
enough that it affects the thermal evolution sufficiently elsewhere in each phase. Further-
more, if there are multiple phase-field parameters to consider (and in our model, there
are), then one should relate these together via physical considerations and measurable pa-
rameters, as we do below. See Kim and Kim (2005); Provatas and Elder (2010) for more
November 28, 2021, 8:23am
X - 4
:
details on the implementation of phase-field models and how to best choose reasonable
phase-field parameters.
To reiterate for the reader’s convenience, we have the following PDE which describes
the evolution of the phase variable,
φ
(where
φ
= 1 for melt and
φ
= 0 for solid phase),
and the temperature,
T
:
τ
∂φ
∂t
+
∇·
(
−
ω
2
φ
∇
φ
)
=
−
d
Γ
dφ
−
L
H
(
T
−
T
m
)
T
m
d
Ψ
dφ
,
(1)
∂T
∂t
+
∇·
(
−
α
∇
T
) =
L
c
p
d
Ψ
dφ
∂φ
∂t
,
(2)
where the relevant parameters are described in Table S1. Equations (1) and (2) above are
completed with the following auxiliary functions: Γ(
φ
) =
φ
2
(1
−
φ
)
2
; Ψ(
φ
) = (3
−
2
φ
)
φ
2
(Provatas & Elder, 2010). The phase-field modeling parameters in our model are
d
,
ω
φ
,
and
τ
.
Rewriting the parameters in Kim and Kim (2005) in terms of the parameters in our
model, we have that
ω
=
H
,
M
φ
=
M
,
ε
=
ε
φ
,
D
T
=
α
, ∆
H
m
=
L
, Γ (
φ
) =
φ
2
(1
−
φ
)
2
,
f
c
= 0, and
f
φ
(
φ,T
) =
(
T
−
T
m
)
L
T
m
d
Ψ
dφ
, where Ψ (
φ
) = (3
−
2
φ
)
φ
2
. Then, from here, we go
through the same derivations in (Kim & Kim, 2005) to derive the interfacial width
d
= 2
ξ
and the interface energy
σ
.
Consider a partially-solidified lava system at equilibrium where we have a 1D interface
between solid
φ
= 1 at
x
=
d
and liquid
φ
= 0 at
x
= 0. Since this system is at equilibrium
and we assume the equal temperature condition for pure substances,
∂φ
∂t
=
∂T
∂t
= 0 and
T
=
T
m
, such that equation (1) from the supplement can be integrated for the equilibrium
November 28, 2021, 8:23am
:
X - 5
phase-field profile
φ
0
(
x
).
ω
2
φ
∂
2
φ
0
∂x
2
−
d
Γ (
φ
0
)
dφ
0
(
T
m
−
T
m
)
T
m
d
Ψ (
φ
0
)
dφ
0
= 0
⇒
ω
2
φ
∂φ
0
∂x
∂
2
φ
0
∂x
2
−
∂φ
0
∂x
d
Γ (
φ
0
)
dφ
0
= 0
⇒
d
dx
[
1
2
ω
2
φ
(
∂φ
0
∂x
)
2
−
Γ (
φ
0
)
]
= 0
⇒
1
2
ω
2
φ
(
∂φ
0
∂x
)
2
−
Γ (
φ
0
) =
1
2
ω
2
φ
(
∂φ
0
∂x
∣
∣
∣
∣
x
=
x
0
)
2
−
Γ (
φ
0
(
x
0
))
,
(3)
where we assume that
ω
φ
is a constant and
x
0
∈
R
is some reference position. Finally,
without loss of generality, we can let
x
0
= 0 and put a Dirichlet boundary condition here
(we would expect one anyways if the lava is fully liquid there), such that Γ (
φ
0
(
x
0
)) =
Γ (
φ
0
(0)) = Γ (0) = 0 and
∂φ
0
∂x
∣
∣
x
=
x
0
=
∂φ
0
∂x
∣
∣
x
=0
= 0. Hence,
1
2
ω
2
φ
(
∂φ
0
∂x
∣
∣
x
=
x
0
)
2
−
Γ (
φ
0
(
x
0
)) =
1
2
ω
2
φ
(0)
2
−
Γ (0) = 0 in (3), in which case we can integrate (3) as follows:
1
2
ω
2
φ
(
∂φ
0
∂x
)
2
= Γ (
φ
0
)
⇒
∂φ
0
∂x
=
√
2
ω
2
φ
Γ (
φ
0
)
⇒
d
=
∫
φ
b
φ
a
dφ
0
√
2
ω
2
φ
Γ (
φ
0
)
=
ω
φ
√
2
∫
φ
b
φ
a
dφ
0
|
φ
0
||
1
−
φ
0
|
.
(4)
As in Kim and Kim (2005), we use
φ
a
= 0
.
1 and
φ
b
= 0
.
9 to integrate (4), whence
d
=
ω
φ
2
√
2 ln 3
,
(5)
which is essentially the same result derived in Kim and Kim (2005).
Next, to obtain the interface energy, we again repeat the steps in Kim and Kim (2005)
by considering an equilibrium system with a cylindrical solid in liquid matrix while main-
taining a diffuse interface between them. Since
ε
2
φ
(which is in units of newtons) effectively
scales with the force field associated with the interface, we get the following:
σ
=
ε
2
φ
∫
∞
−∞
(
dφ
0
dr
)
2
dr
=
√
2
ε
φ
∫
1
0
√
H
Γ (
φ
0
)
dφ
0
=
√
2
ε
φ
√
H
∫
1
0
|
φ
0
||
1
−
φ
0
|
dφ
0
=
ε
φ
3
√
H
2
.
(6)
November 28, 2021, 8:23am
X - 6
:
Making necessary assumptions in the thin interface limit, equation (22) from Kim and
Kim (2005) gives us that
J
=
∫
1
0
h
p
(
φ
) [1
−
h
d
(
φ
)]
√
Γ (
φ
)
dφ
=
∫
1
0
φ
3
(6
φ
2
−
15
φ
+ 10) [1
−
φ
3
(6
φ
2
−
15
φ
+ 10)]
|
φ
||
1
−
φ
|
dφ
=
209
420
.
(7)
Thus, using (7), equation (21) in Kim and Kim implies that
β
=
1
3
√
2
T
m
√
H
ε
φ
LM
−
L
αc
p
ε
φ
√
2
H
J
=
1
3
√
2
T
m
√
H
ε
φ
LM
−
209
420
L
αc
p
ε
φ
√
2
H
.
(8)
Our only unknown parameter is
ω
φ
∼
d
, which we have to adjust as we run simulations
to match known solidification data, but once given
ω
φ
and
d
, we can derive
H
,
τ
, and
ε
φ
.
Therefore, our parameter search is only
one-dimensional
, since once we choose a value of
d
or
ω
φ
, all other parameters can be immediately determined.
Using equations
ω
φ
=
ε
φ
/
√
H
and
τ
= 1
/
(
HM
) from Provatas and Elder (2010) along
with (5), (6), and (8) above, we can rewrite all unmeasured parameters in terms of mea-
surable quantities and
d
as follows:
ω
φ
=
√
2
4 ln 3
d, H
= 12 ln 3
σ
d
, τ
=
1
8 ln
2
3
d
2
L
σT
m
(
β
+
209
1680 ln 3
dL
αc
p
)
.
(9)
Using the sample parameters from Table S1, the last two equations in (9) become
H
≈
6
.
592
d
J m
−
2
, τ
≈
2
.
899
×
10
6
d
3
s m
−
3
+ 3
.
454
×
10
−
6
d
2
s m
−
2
.
(10)
Even if the solid-liquid interface width were microscopic, i.e.,
d
∼
10
−
9
m, the second
equation in (10) would still imply that 3
.
454
×
10
−
6
d
2
s m
−
2
2
.
899
×
10
6
d
3
s m
−
3
,
since in that case,
3
.
454
×
10
−
6
d
2
s m
−
2
2
.
899
×
10
6
d
3
s m
−
3
∼
10
−
3
. Thus, we can further make the simplification
τ
≈
2
.
899
×
10
6
d
3
s m
−
3
, and in general, for parameters similar to basalt lava, the third
November 28, 2021, 8:23am
:
X - 7
equation in (9) can be simplified to
τ
=
209
13440 ln
3
3
d
3
L
2
αc
p
σT
m
.
(11)
Finally, by the above considerations and equations, we can also derive the following
informative scaling properties:
ω
φ
∼
d, H
∼
σ
d
, τ
∼
d
3
L
2
αc
p
σT
m
, M
∼
αc
p
T
m
d
2
L
2
, ε
φ
∼
√
σ
√
d.
(12)
The scaling relationships in (12) provide a physical interpretation for these phase-field
parameters and simple sanity-checks for the validity of the assumptions we made via a
given choice of parameters. And as mentioned in Kim and Kim (2005), the analysis above
should hold in theory as long as
d
α/V
and
d
R
, where
V
and
R
are the velocity
of the solidification front and the local radius of curvature for the solid-liquid interface,
respectively.
Text S2. Initial and boundary conditions and other mathematical details
The domain initially consists of a substrate that is 4
×
h
thick with a uniform ground
temperature of
T
g
= 20
◦
C
. The total domain grows dynamically as lava lobes are emplaced
at temperature
T
0
= 1200
◦
C
:
φ
(
t
= 0) =
{
1
z
∈
[0
,
4
h
)
0
z
∈
[4
h,
5
h
]
,
T
(
t
= 0) =
{
T
g
z
∈
[0
,
4
h
)
T
0
z
∈
[4
h,
5
h
]
The bottom boundary condition is set to a constant temperature of
T
g
and always solid,
assuming that the deep ground maintains a fixed temperature. We choose a domain of size
4
×
h
because the domain must be large enough such that the Dirichlet (fixed) temperature
boundary condition at the bottom does not stop the propagation of the temperature profile
to a non-negligible extent. Hence, in principle, we could use a larger domain, but this
November 28, 2021, 8:23am
X - 8
:
would significantly slow down our numerical simulations with negligible improvement to
accuracy. The top boundary condition is set to lose heat due to black-body radiation,
convection by a fixed wind speed and conduction. We also assume that a crust readily
forms at the top of lava surface:
Bottom boundary :
φ
= 1;
T
=
T
g
;
(13)
Top boundary :
φ
= 1;
k
∂T
∂z
=
−
h
c
(
T
−
T
s
)
−
σ
s
ε
(
T
4
−
T
4
s
)
.
(14)
Here,
T
s
= 30
◦
C,
h
c
,
σ
s
, and
ε
are the surface air temperature, convective heat transfer
coefficient for air flow, Stefan-Boltzmann constant, and the emissivity of the lava surface,
respectively. In practice,
h
c
depends on the wind speed and angle at which it travels with
respect to the lava. However, considering that fluctuations in the external environment
are on a much smaller timescale compared to the solidification timescale, we assume a
constant
h
c
as shown in Table S1, which corresponds to a wind speed of roughly 2 m
/
s
(Patrick et al., 2004).
∇
is defined as
(
∂
∂x
1
,...,
∂
∂x
n
)
in an
n
-dimensional Cartesian coordinate system, and if
we were to use our phase-field model to simulate a solid-liquid system over an arbitrary,
n
-
dimensional domain, then we would have to work with this in full generality. However, for
the simulations in this paper, we are working in one dimension, i.e.,
n
= 1. Thus,
∇
=
∂
∂x
in our case, and since
ω
φ
and
α
are assumed to be constant,
∇·
(
−
ω
2
φ
∇
φ
)
=
−
ω
2
φ
∂
2
φ
∂x
2
and
∇·
(
−
α
∇
T
) =
−
α
∂
2
T
∂x
2
for this paper, respectively. Equations (1) and (2) in the main
paper are written in a more general form for sake of completeness and to demonstrate the
generality of this phase-field model.
November 28, 2021, 8:23am
:
X - 9
Text S3. Numerical scheme
We perform the numerical simulations with a 4th-order
centered difference discretization in space to properly resolve the phase boundary. The
spatial grid size we use ∆
x
roughly scales with
h
, which balances computational efficiency
with numerical precision relative to the lobe size (see text S5). Moreover, we use the
AB4-AM4 predictor–corrector method (Atkinson, 1988; Zlatev, 1985) to integrate in time,
which allows us to increase time step size while ensuring accuracy for the highly-resolved
grid. Because our simulations need to capture temporal dynamics that span from the
order of seconds (initial cooling) to years, we also implement adaptive time-stepping as
monitored with Milne’s device (Atkinson, 1988; Zlatev, 1985; Fujii, 1991). Milne’s device
requires one to prescribe an error tolerance (for our simulations,
tol
= 10
−
11
) as well as
some appropriate norm to measure the error in our grid solution at each timestep; for
our simulations, we use the relative
L
2
error as described in text S4 below. We also use
Ralston’s 4th-order Runge-Kutta method (Ralston, 1962) to predict the first four time
steps after each change in time step size.
Text S4. Definition of relative L
2
error
Say we have
n
data points in space. Suppose that
f
(
x
) is our exact function and
̂
f
(
x
)
is an approximation for
f
. Then, the exact
L
2
error on [0
,L
] would be
e
=
[
∫
L
0
(
f
(
x
)
−
̂
f
(
x
)
)
2
dx
]
1
/
2
.
However, given that
̂
f
lives on a grid with
n
points and spatial intervals of size ∆
x
, we
have to approximate
e
as follows, using a Riemann sum:
e
≈
[
∆
x
n
∑
i
=1
(
f
(
x
i
)
−
̂
f
(
x
i
)
)
2
]
1
/
2
,
November 28, 2021, 8:23am
X - 10
:
where
x
n
=
L
and
x
1
= ∆
x
. Finally, to compute the relative
L
2
error,
e
rel
, we divide
e
by
the
L
2
norm of
f
, i.e.,
e
rel
≈
[
∆
x
∑
n
i
=1
(
f
(
x
i
)
−
̂
f
(
x
i
)
)
2
]
1
/
2
[
∆
x
∑
n
i
=1
f
(
x
i
)
2
]
1
/
2
=
∑
n
i
=1
(
f
(
x
i
)
−
̂
f
(
x
i
)
)
2
∑
n
i
=1
f
(
x
i
)
2
1
/
2
.
The relative
L
2
error is also what we use to control the timestep for Milne’s device, in
which one lowers or increases the timestep to keep the prescribed error within a certain
range. The motivation is that if the timestep is unnecessarily small, then the numerical so-
lution will evolve too slowly and become computationally inefficient, but if the timestep is
too large, then discretization errors will begin to dominate the solution and the numerical
solution could also blow up. For further details, see Atkinson (1988).
Text S5. Spatial grid size
∆
x
The spatial grid size we use is
∆
x
= 10
−
3
min
{
nint
(
√
10
h
)
,
10
}
,
(15)
where nint is the function which rounds its argument to the nearest integer. Intuitively,
we can think of equation 15 as interpolating ∆
x
from 10
−
3
(for
h
= 0
.
1
,
0
.
2) to 10
−
2
(for
h
= 10
,
15
,
20) using the square root function, except rounding each value to the third
decimal place for simplicity’s sake. That way, ∆
x
roughly scales with
h
, which balances
computational efficiency with numerical precision.
Text S6: Different cooling features in lava flows
As a lava flow solidifies, it devel-
ops contraction fractures called columnar joints which propagate inward from the cooling
front (Spry, 1962; Long & Wood, 1986; Phillips et al., 2013; Forbes et al., 2014). The
shape of these columnar joints can vary from the classical regularly spaced hexagonal
November 28, 2021, 8:23am
:
X - 11
shape to more irregular shapes. If the joints are dominantly vertical (indicative of nearly
horizontal isotherms), and have a coarse, quasi-regular pattern, they are called a colon-
nade. In contrast, if the joints are chaotic/irregular, finely spaced, and/or have multiple
orientations, they are called entablature. Typically, these are indicative of faster cooling
rates with complex cooling front shapes due to rainfall or stream flow (for more discussion
see Spry, 1962; Long & Wood, 1986; Phillips et al., 2013; Forbes et al., 2014). A key
characteristic of potentially fused flow lobes is that they have multiple zones of entabla-
ture/colonnade within a single thick flow lobe suggesting that multiple parts of this flow
lobe experienced faster/slower cooling. This observation can be naturally explained if the
present day single flow lobe was in fact emplaced as multiple lobes.
Text S7: Comparing the model with data from Hon et al. (1994).
Figures S2-S3
show how well our model fits the data given in Figures 8-10 from Hon et al. (1994). Figure
8 in Hon et al. (1994) shows the temperatures at fixed depths, flow thickness, and crack
depth vs. time measured from an actively inflating sheet flow emplaced on 17 April 1990
on K ̄ılauea for up to just over 100 hours, but we only compared with the temperatures
here because our model does not include inflation. Figure 9 in Hon et al. (1994) tracks the
thickness of the solidified crust as a function of the square root of time at isotherms 800
◦
C
and 1070
◦
C for up to around 60 hours, using data from both the K ̄ılauea flow data from
Hon et al. (1994) and the Makaopuhi lava lake data from Wright and Okamura (1977).
Finally, Figure 10 in Hon et al. (1994) shows the same data as in Figure 8, except with a
log scaled horizontal time axis and only tracking temperature. For the simulations shown
in both Figures S2 and S3, we use an initial temperature of
T
0
= 1142
◦
C as given in Hon
November 28, 2021, 8:23am
X - 12
:
et al. (1994) (as opposed to
T
0
= 1200
◦
C mentioned in text S2 and used in the rest of our
simulations), take the values of all physical parameters aside from
α
to be those in Table
S1, and set the ground temperature
T
g
= 1200
◦
C, since the temperature propagation at
the bottom boundary does not need to be tracked for comparisons.
Figures S2 and S3 are the same, except that we use a different value of
α
in the model for
generating our results; Figure S2 uses the same thermal diffusivity
α
= 3
.
75
×
10
−
7
m
2
/s
from Table S1 and in the main results from our paper, while Figure S3 uses
α
= 2
.
9
×
10
−
7
m
2
/s. We do this because, as mentioned in section 2.2 of the main text, near the surface,
there are various physical processes which could lead to a lower diffusivity/conductivity in
the melt, such as surface effects and the presence of bubbles. Furthermore, as mentioned
in Touloukian et al. (1971), the value of
α
for basalt could range by around 50% regardless
of such effects, depending on temperature or chemical composition of the particular basalt
melt in question. For the value of
α
used in Figure S3, our model agrees noticeably better
with the empirical results given Figures 8-10 from Hon et al. (1994) than the one used in
Figure S2, and the value of
α
in the former, while lower than what many sources use, is
still within the reasonable range described in Touloukian et al. (1971).
In either case for
α
, our phase-field model breaks down within
≈
10 cm of the lava’s
surface; the concave computed temperature curves in Figures S2-S3 gradually approach
the convex shape of the one at the surface, while the measured temperature curves right
below the surface remain concave and more discontinuously become convex at the surface.
Furthermore, since the thermocouples used to take measurements in Hon et al. (1994) were
also within only a few centimeters of each other in the inflating lobe, it is also plausible
November 28, 2021, 8:23am
:
X - 13
that measurement errors could be responsible for some discrepancies between our model
and the data. Finally, we note that we are only considering lava flow cooling in the limit
of low undercooling, since we do not consider cooling rates fast enough such that flash
heating and recalescence are important. We expect that this simplification should not
significantly affect our conclusions, since the impact of this process, e.g., as described
in Whittington and Sehlke (2021), is only relevant for very fast cooling rates (
∼
10s of
◦
C/s) and when the melt becomes undercooled by tens to hundreds of degrees before
crystallization begins. For lava flows, this process would only be relevant for at most the
skin of lava flows. The phase field model provides a natural framework to incorporate the
disequilibrium crystallization processes in future work.
Figure S4 Caption.
Using nonlinear least squares, we fit the solidification time for a
single
lobe to cool,
t
h
, by itself to the function
t
h
=
h
2
α
[
A
h
B
+
C
exp (
−
Dh
) +
E
]
(16)
which heuristically models the nonlinear trend as deviations from the conventional cooling
estimate
t
h
∼
h
2
derived from solving the Stefan problem. The best fit parameters we
find are
A
≈
0
.
0110,
B
≈
0
.
2294,
C
≈
0
.
3346,
D
≈
24
.
8922, and
E
≈
0
.
0320. With these
parameters, the relative
L
2
error (as defined in text S4) between
t
h
α/h
2
as fitted above and
the actual data is
≈
9
.
8237
×
10
−
3
<
1%, which shows very respectable agreement. Hence,
(16) could be a starting point for modeling
t
h
with more general physical parameters and
initial conditions.
By using the term “strong nonlinearity” in Figure S4, we are referring to how there is
a qualitative difference in the curve for small enough lobe sizes. This difference is best
November 28, 2021, 8:23am
X - 14
:
explained by the quick decay of the exponential term in our curve fit: For
h
not too large,
the exponential term quickly disappears and the trend becomes primarily dominated by
the power law term. Hence, motivated by how the relative
L
2
error between our best-
fit curve and numerical solution is just under 1%, we heuristically have drawn the line
between the “strongly nonlinear” and “weakly nonlinear” regions by indicating where the
relative error between the fit with and without the exponential term falls below 1%. That
point is at roughly
h
= 0
.
26344, after which the exponential term contributes an error
which is below 1% and decreases further as
h
increases.
We label these two regions in Figure S4 to give a rough estimate of where the usual
t
h
∼
h
2
scaling relationships are relatively valid, and show how for small enough lobe sizes,
deviations from this trend begin to dominate significantly. The physical interpretation
of these regions is as follows: As we work with smaller and smaller lobes, the nonlinear
effects of convective cooling and radiative heat loss at the lava surface begin to dominate
the time it takes for a lobe of that size to cool. The usual Stefan problem formulation
often ignores these nonlinear effects in the boundary condition at the lava-air interface,
but based off of our results here, we suggest that these will contribute a non-negligible
effect to the solution when the lobe size is too small.
Notes for Figures S5-S9.
For Figures S5-S9, we consider the trends between different
lobe thicknesses once we weight the emplacement time by
t
h
, the time it takes for a
single lobe of thickness
h
to solidify completely. For each of these dimensionless plots,
the stars represent fused cases, the diamonds represent in parallel cases, and the boxes
represent in sequence cases, while the colors of these markers differentiate the lobe sizes
h
.
November 28, 2021, 8:23am
:
X - 15
Moreover, as the legend (which is the same for Figures S5-S9) indicates, the shape of the
markers distinguishes the type of flow which occurred (fused, in parallel, or in sequence,
as described in the main text under section 3), while the color of the markers distinguishes
h
, the thickness of both the initial lobe and the one being emplaced.
Figure S5 Caption.
(Compare with Figure S6.) The dimensionless log-log plot shows
t
solidification
, the time it takes for the entire two-lobe system to solidify (including the time
before the second lobe is emplaced), as a function of the emplacement time interval,
t
emp
,
with both axes scaled by
t
h
, the time it takes for a single lobe of thickness
h
to solidify.
Note in particular that the graph at any lobe size has a minimum near or slightly below
t
emp
=
t
h
. This minimum reflects some optimal balance between the emplacement time
and the thermal/phase interaction between the two lobes which minimizes the solidifica-
tion time across the domain; this optimal balance lies within the in parallel regime.
Figure S6 Caption.
(Compare with Figure S5.) The plot highlights an alternate in-
terpretation of the two-lobe solidification time in which we neglect the time between em-
placements; like S5, this plot also shows the time it takes for the entire two-lobe system to
solidify (this time
without
the time before the second lobe is emplaced), as a function of
the emplacement time interval,
t
emp
, with both axes scaled by
t
h
. On either plot, we note
that as
t
emp
→
0,
t
solidification
→
4
t
h
. This reflects how, since
t
h
∼
h
2
,
t
2
h
∼
(2
h
)
2
= 4
h
2
.
Meanwhile, in comparison to Figure S4, this plot better demonstrates how, as
t
emp
→∞
,
t
solidification
→
t
h
+
t
emp
.
Figure S7 Caption.
The plot only considers the time for the lower (i.e., earlier) lobe
to solidify (in this case also represented by
t
solidification
) vs.
t
emp
/t
h
, thereby highlighting
November 28, 2021, 8:23am
X - 16
:
the thermal influence of the upper lobe upon how the lower lobe solidifies relative to
t
h
.
As expected,
t
solidification
→
t
h
when
t
emp
→∞
. Physically, we can interpret this result as
follows: If the lower lobe has fully solidified before the upper lobe is emplaced, then the
upper lobe will have no influence on the solidification of the lower lobe, provided that no
remelting can occur.
Figure S8 Caption.
This plot indicates the height
z
, scaled to the lobe size, at which
solidification completed across the entire two-lobe system vs.
t
emp
/t
h
. This variable is
significant because horizontal fractures often form wherever solidification completes in a
lava lobe, i.e., where two solidifying fronts meet. Note in particular that the smaller lobe
sizes appear to have greater solidification heights in the merged and in sequence regimes,
while the opposite behavior is observed for the in parallel regime. The quantitative dif-
ferences in behavior across different lobe sizes appears to be greatest for the in parallel
regime, which we also see in Figure S9.
Figure S9 Caption.
This plot is the same as Figure S8, except that it measures the
height at which solidification occurred in the first lobe only, i.e., where the first lobe
solidified. Note that for a given height, the graph appears to increase in the fused regime,
decrease sharply in the in parallel regime, and then finally level out in the in sequence
regime. The trend in the in parallel regime appears to decrease sharper the smaller the
lobe size is, which indicates how the height where solidification occurred in the lower lobe
seems to move closer to the upper lobe the smaller the emplaced lobes were. This is
because the thermal influence of the upper lobe on the lower lobe increases as the lobe
size decreases, assuming that the lobes do not just merge entirely. As with Figure S8, the
November 28, 2021, 8:23am
:
X - 17
greatest disparity in dynamics across different lobe sizes appears to be greatest for the in
parallel regime.
Movie S1.
Under the folder
movies
in the Zenodo data repository (
https://zenodo
.org/badge/latestdoi/357729300
),
emplacementresults_10_10_843K_406hours.mp4
is a movie showing the solidification dynamics for the fused flow case shown in Figure 1.
Movie S2.
Under the folder
movies
,
emplacementresults_10_10_843K_1625hours.mp4
is a movie showing the solidification dynamics for the in parallel case shown in Figure 1.
Movie S3.
Under the folder
movies
,
emplacementresults_10_10_843K_26000hours.mp4
is a movie showing the solidification dynamics for the in sequence case shown in Figure 1.
References
Atkinson, K. E. (1988).
An introduction to numerical analysis
. New York: Wiley.
Forbes, A. E., Blake, S., & Tuffen, H. (2014). Entablature: fracture types and mecha-
nisms.
Bulletin of Volcanology
,
76
(5), 1–13. doi: 10.1007/s00445-014-0820-z
Fujii, M. (1991). An extension of Milne’s device for the Adams Predictor-Corrector
Methods.
Japan Journal of Industrial and Applied Mathematics
,
8
(1), 1–18. doi:
10.1007/BF03167183
Hon, K., Kauahikaua, J., Denlinger, R., & Mackay, K. (1994). Emplacement and inflation
of p ̄ahoehoe sheet flows: observations and measurements of active lava flows on
K ̄ılauea Volcano, Hawai‘i.
Geological Society of America Bulletin
,
106
(3), 351–370.
doi: 10.1130/0016-7606(1994)106
〈
0351:EAIOPS
〉
2.3.CO;2
Kim, S. G., & Kim, W. T. (2005). Phase-Field Modeling of Solidification. In
Handbook of
materials modeling
(pp. 2105–2116). Springer, Dordrecht. doi: 10.1007/978-1-4020
November 28, 2021, 8:23am
X - 18
:
-3286-8
109
Long, P. E., & Wood, B. J. (1986). Structures, textures, and cooling histories of columbia
river basalt flows.
Geological Society of America Bulletin
,
97
(9), 1144–1155. doi:
10.1130/0016-7606(1986)97
〈
1144:STACHO
〉
2.0.CO;2
Patrick, M. R., Dehn, J., & Dean, K. (2004). Numerical modeling of lava flow cooling
applied to the 1997 Okmok eruption: Approach and analysis.
Journal of Geophysical
Research
,
109
(B3), B03202. doi: 10.1029/2003JB002537
Phillips, J. C., Humphreys, M. C., Daniels, K., Brown, R., & Witham, F. (2013).
The formation of columnar joints produced by cooling in basalt at Staffa, Scotland.
Bulletin of Volcanology
,
75
(6), 1–17. doi: 10.1007/s00445-013-0715-4
Provatas, N., & Elder, K. (2010).
Phase-field methods in materials science and engineer-
ing
(1st ed.). Wiley. doi: 10.1002/9783527631520
Ralston, A. (1962). Runge-Kutta methods with minimum error bounds.
Mathematics of
Computation
,
16
(80), 431–437. doi: 10.1090/S0025-5718-1962-0150954-0
Sheth, H. (2017). Morphology and architecture of flood basalt lava flows and sequences.
In
A photographic atlas of flood basalt volcanism
(pp. 33–79). Springer International
Publishing. doi: 10.1007/978-3-319-67705-7
3
Spry, A. (1962). The origin of columnar jointing, particularly in basalt flows.
Journal of
the Geological Society of Australia
,
8
(2), 191–216. doi: 10.1080/14400956208527873
Touloukian, Y., Powell, R., Ho, C., & Klemens, P. (1971).
Thermophysical properties
of matter - The TPRC data series. Volume 2. Thermal conductivity - Nonmetallic
solids
. Defense Technical Information Center.
November 28, 2021, 8:23am
:
X - 19
Whittington, A. G., & Sehlke, A. (2021). Spontaneous reheating of crystallizing lava.
Geology
. doi: 10.1130/G49148.1
Wright, T. L., & Marsh, B. (2016). Quantification of the intrusion process at K ̄ılauea
Volcano, Hawai‘i.
Journal of Volcanology and Geothermal Research
,
328
, 34–44. doi:
10.1016/j.jvolgeores.2016.09.019
Wright, T. L., & Okamura, R. T. (1977).
Cooling and crystallization of tholeiitic basalt,
1965 Makaopuhi Lava Lake, Hawai‘i
(Tech. Rep.). doi: 10.3133/pp1004
Wright, T. L., Peck, D. L., & Shaw, H. R. (1976). K ̄ılauea lava lakes: Natural labo-
ratories for study of cooling, crystallization, and differentiation of basaltic magma.
In
The geophysics of the pacific ocean basin and its margin
(p. 375-390). American
Geophysical Union (AGU). doi: 10.1029/GM019p0375
Zlatev, Z. (1985). Variable stepsize variable formula methods based on predictor-
corrector schemes.
Applied Numerical Mathematics
,
1
(5), 395–416. doi: 10.1016/
0168-9274(85)90003-0
November 28, 2021, 8:23am
X - 20
:
Table S1.
Parameters used for the model.
Definition
Unit
Values used
L
latent heat of fusion
J/m
3
4
×
10
5
c
p
specific heat at constant pressure
J/(m
3
·
K)
2
.
57
×
10
6
k
thermal conductivity
J
/
(m
·
s
·
K)
9
.
64
×
10
−
1
T
m
melting temperature
o
C
1070
α
thermal diffusivity
m
2
/s
3
.
75
×
10
−
7
τ
characteristic time of solidification
s
2
.
90
×
10
6
ω
φ
interfacial width coefficient
m
3
.
22
×
10
−
1
ε
φ
interfacial force coefficient
√
N
8
.
26
×
10
−
1
d
interfacial width
m
1
.
00
σ
interfacial energy
J/m
2
5
×
10
−
1
H
energy barrier
J/m
3
6
.
59
β
kinetic coefficient
m
·
Pa/K
2
5
.
6
×
10
−
8
h
c
convective heat transfer coefficient of air W/(m
2
·
K) 2
.
62
×
10
1
σ
s
Stefan-Boltzmann constant
W/(m
2
·
K
4
) 5
.
67
×
10
−
8
ε
emissivity of the lava surface
-
0.6
November 28, 2021, 8:23am