Supporting Information: A phase-field model for
wet snow metamorphism
Adrian Moure
∗
and Xiaojing Fu
∗
Department of Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, CA, United States
E-mail: amoure@caltech.edu; rubyfu@caltech.edu
1
1 Equivalence to models of two-phase transition
2
In this section we show that our wet snow metamorphism model is equivalent to the two-phase
3
models for solidification, sublimation, and evaporation when only two phases are present in
4
the system. Only for this section, we assume that
θ
nucl
= 0 (see Eqs. (1)–(3) in the main
5
text). To prove the equivalency, we first need to show that a two-phase system remains as a
6
two-phase system and the non-existing phase does not appear in Ω. It is trivial to show that
7
the phase-fields equations in our model, Eqs. (1)–(3) in the main text, have the property
8
φ
l
(
x
,t
) = 0
∀
t
if
φ
l
(
x
,
0) = 0 and
θ
nucl
= 0, which proves this first condition. The second
9
condition constitutes the actual equivalence, which we analyze for each phase transition in
10
the following paragraphs.
11
1
1.1 Solidification
12
Solidification and melting occurs when
φ
a
= 0. We can rewrite the
φ
i
evolution equation
13
(Eq. (1) in the main text) imposing
φ
a
= 0 as
14
∂φ
i
∂t
=
−
3
M
0
1
ε
φ
i
(1
−
φ
i
)(1
−
2
φ
i
)
−
ε
∇
2
φ
i
−
α
sol
φ
2
i
φ
2
w
T
−
T
melt
L
sol
/c
p,w
.
(1)
The Lagrange multiplier in our formulation guarantees that
φ
i
+
φ
w
= 0, which implies that
15
∂φ
w
∂t
=
−
∂φ
i
∂t
. Imposing
φ
a
= 0, the temperature evolution equation (Eq. (12) in the main
16
text) is expressed as:
17
ρ
(
φ
)
c
p
(
φ
)
∂T
∂t
=
∇·
K
(
φ
)
∇
T
+
ρ
(
φ
)
L
sol
∂φ
i
∂t
.
(2)
We can disregard the vapor density equation since there is no air phase. Equations (1) and (2)
18
constitute a two-phase model for solidification (and melting) known as the generalized Stefan
19
problem.
1,2
Under certain parameter regimes, this model reproduces the kinetics defined by
20
the Gibbs-Thomson condition for the ice-liquid water interface, which reads
21
T
−
T
melt
L
sol
/c
p,w
=
−
d
sol
χ
−
β
sol
v
n
,
(3)
where
T
is the interface temperature,
β
sol
is the kinetic attachment coefficient,
d
sol
is the
22
capillary length,
χ
is the curvature of the interface (positive for spherical ice grains), and
23
v
n
is the normal velocity of the interface (positive for ice growth). Karma and Rappel
1,3
24
showed that the two-phase model captures the Gibbs-Thomson condition kinetics when (i)
25
ε
→
0 and (ii) the model parameters
M
0
and
α
sol
follow the relations
26
M
0
=
M
sol
=
ε
3
τ
sol
, α
sol
=
λ
sol
τ
sol
, d
sol
=
a
1
ε
λ
sol
, β
sol
=
a
1
τ
sol
ελ
sol
−
a
2
ε
D
∗
iw
,
(4)
2
where
a
1
≈
5,
a
2
≈
0
.
1581
1
, and
D
∗
iw
≈
(
D
i
+
D
w
)
/
2 represents thermal diffusivity, such
27
that
D
j
=
K
j
/
(
ρ
j
c
p,j
) for
j
=
{
i,w
}
. However, we note that the above relations assume
28
equal thermal diffusivity around the freezing/melting interface, which is not true for ice
29
and water in our case. Asymptotic analysis assuming asymmetric diffusivities have been
30
conducted for compositional problems,
4,5
but the same type of analysis for the thermal
31
problems remain to be developed. Nevertheless, Kaempfer and Plapp
6
has shown that the
32
phase-field solidification model recovers the sharp-interface solution if
ε < β
sol
D
∗
iw
≈
10
−
5
m,
33
which is true in the simulations presented in this paper.
34
1.2 Sublimation
35
Sublimation and deposition occurs when
φ
w
= 0. After imposing
φ
w
= 0, Eqs. (1), (12) and
36
(14) in the main text are expressed as
37
∂φ
i
∂t
=
−
3
M
0
1
ε
φ
i
(1
−
φ
i
)(1
−
2
φ
i
)
−
ε
∇
2
φ
i
+
α
sub
φ
2
i
φ
2
a
ρ
v
−
ρ
I
vs
(
T
)
ρ
i
(5)
ρ
(
φ
)
c
p
(
φ
)
∂T
∂t
=
∇·
K
(
φ
)
∇
T
−
ρ
(
φ
)
L
sub
∂φ
a
∂t
,
(6)
∂
(
φ
a
ρ
v
)
∂t
=
∇·
[
φ
a
D
v
(
T
)
∇
ρ
v
] +
ρ
SE
∂φ
a
∂t
,
(7)
with
∂φ
a
/∂t
=
−
∂φ
i
/∂t
. These equations constitute the two-phase sublimation model
6
and
38
capture the kinetics defined by the Gibbs-Thomson condition
39
ρ
v
−
ρ
I
vs
(
T
)
ρ
I
vs
(
T
)
=
d
sub
χ
+
β
sub
v
n
,
(8)
1
The
a
1
and
a
2
values are computed by following the same procedure as in Karma and Rappel
3
, applied
to Eqs. (1) and (2).
3
as long as the parameters
M
0
,
α
sub
, and
ρ
SE
follow the relations:
6
40
M
0
=
M
sub
=
ε
3
τ
sub
, α
sub
=
λ
sub
τ
sub
, ρ
SE
=
ρ
i
,
d
sub
ρ
I
vs
ρ
i
=
a
1
ε
λ
sub
, β
sub
ρ
I
vs
ρ
i
=
a
1
τ
sub
ελ
sub
−
a
2
ε
D
∗
ia
−
a
2
ε
D
v
0
,
(9)
where
β
sub
and
d
sub
are the Gibbs-Thomson coefficients for sublimation and
D
∗
ia
represents
41
the average thermal diffusivity of ice and air. Again, we note that the above relations
42
assume equal thermal diffusivity around the interface, which is not true for ice and air in
43
our case. Asymptotic analysis assuming for this exact problem remains to be developed.
44
Meanwhile, Kaempfer and Plapp
6
has shown that the above phase-field model recovers the
45
sharp-interface solution if
ε < β
sub
ρ
I
vs
ρ
i
D
∗
ia
and
ε < β
sub
ρ
I
vs
ρ
i
D
v
0
. These two conditions
46
yield the requirement that
ε <
10
−
7
m, which we do not satisfy in the simulations presented.
47
Thus, our numerical results do not strictly recover the Gibbs-Thomson condition for air-ice
48
interfaces.
49
1.3 Evaporation
50
The derivation for the evaporation case (
φ
i
= 0) is analogous to the sublimation case. The
51
final two-phase model for evaporation (without flow) is analogous to the 2-phase sublimation
52
model (Eqs. (5)–(7)), where the evaporation latent heat is expressed as (
L
eva
=
L
sub
−
L
sol
).
53
Here, we simply list the final relations between the model parameters
M
0
,
α
eva
,
ρ
SE
, and the
54
Gibbs-Thomson coefficients
β
eva
and
d
eva
, which read
55
M
0
=
M
eva
=
ε
3
τ
eva
, α
eva
=
λ
eva
τ
eva
, ρ
SE
=
ρ
w
,
d
eva
ρ
W
vs
ρ
w
=
a
1
ε
λ
eva
, β
eva
ρ
W
vs
ρ
w
=
a
1
τ
eva
ελ
eva
−
a
2
ε
D
∗
wa
−
a
2
ε
D
v
0
,
(10)
where
D
∗
wa
represents the average thermal diffusivity of water and air. We note again that
56
the above relations assume equal thermal diffusivity around the interface, which is not true
57
4
for water and air in our case. Asymptotic analysis assuming for this exact problem remains to
58
be developed. Meanwhile, Kaempfer and Plapp
6
has shown that the above phase-field model
59
recovers the sharp-interface solution if
ε < β
eva
ρ
W
vs
ρ
w
D
∗
wa
and
ε < β
eva
ρ
I
vs
ρ
w
D
v
0
. These two
60
conditions yield the requirement that
ε <
10
−
8
m, which we do not satisfy in the simulations
61
presented. Thus, our numerical results do not strictly recover the Gibbs-Thomson condition
62
for air-water interfaces.
63
2 Additional details for numerical implementation
64
2.1 Determining temporal scaling coefficients
65
Wet snow metamorphism involves processes at different time scales. The characteristic times
66
for solidification (
t
sol
), evaporation (
t
eva
), sublimation (
t
sub
), thermal diffusion (
t
T,i
, t
T,w
,
and
t
T,a
),
67
and vapor diffusion (
t
v
), assuming a characteristic length scale of
L
= 10
−
4
m, are
68
t
sol
≡
L/v
n,
sol
∼
1 s
, t
sub
≡
L/v
n,
sub
∼
10
4
s
, t
eva
≡
L/v
n,
eva
∼
10
3
s
, t
T,i
≡
L
2
/D
i
∼
10
−
2
s
,
t
T,w
≡
L
2
/D
w
∼
10
−
1
s
, t
T,a
≡
L
2
/D
a
∼
10
−
3
s
, t
v
≡
L
2
/D
v
0
∼
10
−
4
s
,
(11)
where thermal diffusivities are calculated as
D
l
=
K
l
/
(
ρ
l
c
p,l
), for phase
l
, and we consider
69
typical interface velocities
v
n,j
, for
j
=
{
sol
,
sub
,
eva
}
(see Table S1). These values show
70
that sublimation and evaporation are several orders of magnitude slower than solidification
71
and any diffusion process. A monolithic numerical solver requires the use of time steps
72
that capture the lower characteristic time (i.e.,
t
v
), which would significantly increase the
73
computational times.
74
Here, we leverage the procedure explained in Kaempfer and Plapp
6
to speed up the
75
simulations. This approach leverages the fact that
T
and
ρ
v
exhibit a quasi-steady state
76
compared to the phase transition kinetics. The procedure consists of multiplying the right-
77
hand side of the
T
and
ρ
v
equations (Eqs. (12) and (14) in the main text, respectively) by
78
5
a time-scale factor, so that the characteristic times for vapor and thermal diffusion increase
79
some orders of magnitude while
T
and
ρ
v
keep displaying a quasi-steady state. This approach
80
enables the use of larger time steps without introducing noticeable errors in the simulations.
81
We multiply the right-hand side of the
ρ
v
equation by the time-scale factor
ξ
v
<
1, which
82
provokes an increase of the vapor diffusion characteristic time to
t
∗
v
=
t
v
/ξ
v
. As long as
t
∗
v
83
does not exceed
t
sol
,
t
sub
, and
t
eva
, vapor concentration displays a quasi-steady state and
84
the overall problem kinetics is not affected. This condition leads to
ξ
v
>
10
−
4
. For the
T
85
equation we consider the largest thermal diffusion characteristic time,
t
T,w
. In Eq. (11) we
86
observe that
t
T,w
∼
t
sol
, while
t
T,w
< t
sub
and
t
T,w
< t
eva
. Thus,
T
displays a quasi-steady
87
state when solidification does not occur. We follow the same procedure and multiply the
88
right-hand side of the
T
equation by the time-scale factor
ξ
T
<
1 when solidification does
89
not take place. As long as the new characteristic time
t
∗
T,w
=
t
T,w
/ξ
T
does not exceed
t
sub
90
and
t
eva
, temperature displays a quasi-steady state and the problem kinetics is not affected.
91
This condition leads to
ξ
T
>
10
−
4
. In this work, we take
ξ
v
= 10
−
3
and
ξ
T
= 1 or
ξ
T
= 10
−
2
92
depending on whether or not solidification occurs in any parts of the domain.
93
2.2 Kinetic parameters constraints on grid resolution
94
If we assume equal diffusivity for all phases, then the relations between the Gibbs-Thomson
95
parameters (
β
j
,d
j
) and our model parameters (
M
j
,α
j
) (for
j
=
{
sol
,
sub
,
eva
}
) can be derived
96
from asymptotic analysis
1,3
and are valid under the following conditions:
97
1
. χε <
1
,
2
. εv
n,j
/D
k
<
1
,
3
. λ
j
u
j
<
1
,
4
. τ
j
v
n,j
/ε <
1
,
(12)
for
k
=
{
i,w,a,v
0
}
, where
ε
represents the interface width,
D
k
is the thermal or vapor diffu-
98
sivity coefficient,
v
n,j
is the normal velocity of the interface,
λ
j
and
τ
j
are parameters defined
99
in Section 1 in the Supporting Information, and
u
j
is the phase-change factor that multiplies
100
the terms
α
j
φ
2
l
φ
2
m
in the phase-field equations (for
{
l,m
}
=
{
i,w,a
}
; see Eqs. (23) and (24)
101
6
Table S1: Gibbs-Thomson parameters and
ε
limitations. Index
j
stands for solidification,
sublimation and evaporation. The
β
j
,
d
j
,
M
j
, and
α
j
values considered in this work are
listed in the bottom rows.
Solidification
Sublimation
Evaporation
β
j
(s/m)
∼
(10
,
100)
(2
,
200)
×
10
4
(2
,
200)
×
10
3
d
j
(m)
∼
5
×
10
−
10
∼
10
−
9
∼
5
×
10
−
10
u
j
(-)
∼
10
−
3
∼
10
−
8
∼
10
−
8
v
n,j
(m/s)
∼
(10
−
4
,
10
−
5
)
∼
(10
−
7
,
10
−
9
)
∼
(10
−
6
,
10
−
8
)
Cond. 1 (m)
ε <
∼
10
−
4
ε <
∼
10
−
4
ε <
∼
10
−
4
Cond. 2 (m)
ε <
∼
10
−
3
ε <
∼
10
ε <
∼
10
−
1
Cond. 3 (m)
ε <
∼
10
−
7
ε <
∼
10
−
7
ε <
∼
10
−
8
Cond. 4 (m)
ε <
∼
10
−
7
ε <
∼
10
−
7
ε <
∼
10
−
7
Parameter values considered in this work (lead to
ε <
10
−
6
)
β
j
(s/m)
125
1
.
4
×
10
5
1
.
53
×
10
4
d
j
(m)
4
×
10
−
9
10
−
7
6
×
10
−
8
M
j
(m/s)
2
.
12
×
10
−
5
4
.
33
×
10
−
7
1
.
34
×
10
−
6
α
j
(1/s)
7
.
96
×
10
4
1
.
3
×
10
7
6
.
69
×
10
7
in the main text). Condition 1 imposes an interface width smaller than the curvature radius,
102
condition 2 guarantees that thermal and vapor diffusion are faster than the interface mo-
103
tion, and conditions 3 and 4 ensure the preservation of a tanh-profile (i.e., quasi-equilibrium
104
profile) across the interface. We use these conditions to guide the choice of
ε
in our model;
105
however, we recognize that these conditions are only valid for equal diffusivity problems,
106
which is not true for the wet snow problem. These conditions impose a restriction on
ε
,
107
which represents the numerical interface width and is an adjustable model parameter.
108
To find a valid range for
ε
, we first estimate the Gibbs-Thomson parameters
β
j
and
d
j
109
from the literature
6–9
(see Table S1). Note that
β
j
may range several orders of magnitude
110
depending on ambient conditions and growth type.
8
We next estimate the typical values of
111
u
j
and
v
n,j
(Table S1), based on
T
and
ρ
v
distributions observed during the different phase
112
transitions. We use the Gibbs-Thomson conditions defined in Section 1 of the Supporting
113
Information to relate
u
j
and
v
n,j
, under the assumption that the curvature
χ
≈
0. In order
114
to express
λ
j
and
τ
j
as a function of
β
j
and
d
j
(see Eqs. (4), (9), and (10)), we assume
115
ρ
I
vs
/ρ
i
=
ρ
W
vs
/ρ
w
= 5
×
10
−
6
, which is the approximate value at the freezing point. For
116
7
condition 1, we consider typical ice curvatures
χ
∼
10
4
m
−
1
. With this information, we solve
117
for
ε
in Eq. (12). We indicate the less favorable valid range of
ε
in Table S1. We find that
118
condition 3 is the most restrictive, such that
ε <
10
−
8
m.
119
In order to speed up the simulations, we consider a value of
d
sol
,
d
eva
and
d
sub
one
120
or two orders of magnitude higher than its actual value (see bottom rows of Table S1).
121
This is because the
T
- and
ρ
v
-terms in the Gibbs-Thomson conditions (Eqs. (3) and (8))
122
dominate the interface motion compared to the curvature terms (
d
j
) for
χ
∼
10
4
m
−
1
in
123
most cases. Doing so, the
ε
constraint in Eq. (12) relaxes to
ε <
10
−
6
m, which allows us
124
to use a coarser spatial discretization, whereas the wet snow metamorphism kinetics is not
125
significantly affected. If compared against experiments, deviations in our modeled kinetics
126
from the observed kinetics could be attributed to this numerical treatment and, hence, the
127
actual values of
d
sol
,
d
eva
and
d
sub
must be considered.
128
2.3 Parameter values
129
We list the value of the parameters used in our model at the bottom of Table S1 and in
130
Table S3. For simplicity, we assume that the values of these parameters are constant. The
131
values of the parameters
k
j
and
g
j
in the saturated vapor pressures are listed in Table S2.
132
The parameter
ξ
T
is defined as
133
ξ
T
=
1
if Γ
iw
> ε,
10
−
2
if Γ
iw
≤
ε,
(13)
where Γ
iw
is the length (in 2D) of the ice-water interface and the parameter
ε
is close to
134
zero. In particular, we compute the interface as Γ
iw
=
R
Ω
φ
2
i
φ
2
w
d
x
and we take
ε
= 2
×
10
−
4
ε
135
in 1D and
ε
= 2
×
10
−
6
ε
√
S
Ω
in 2D, where
S
Ω
is the area of Ω.
136
8
Table S2: Saturated vapor pressure of ice and water. Interpolation parameter values.
10,11
j
0
1
2
3
k
j
−
0
.
5865
×
10
4
0
.
2224
×
10
2
0
.
1375
×
10
−
1
−
0
.
3403
×
10
−
4
g
j
−
0
.
2991
×
10
4
−
0
.
6017
×
10
4
0
.
1888
×
10
2
−
0
.
2836
×
10
−
1
j
4
5
6
7
k
j
0
.
2697
×
10
−
7
0
.
6918
-
-
g
j
0
.
1784
×
10
−
4
−
0
.
8415
×
10
−
9
0
.
4441
×
10
−
12
0
.
2859
×
10
1
Table S3: Parameter values. The values of
β
j
,
d
j
,
M
j
, and
α
j
(where
j
stands for sol, sub,
and eva) are listed in the bottom of Table S1. We take
ε
= 2
×
10
−
7
m for the simulations
shown in Section 4.3.3 in the main text.
Parameter
Value
Units
ε
5
×
10
−
7
m
Λ
1
N/m
σ
iw
0.033
N/m
σ
ia
0.109
N/m
σ
wa
0.076
N/m
T
melt
0
◦
C
P
a
101325
Pa
α
nucl
5
×
10
6
1/m
T
nucl
variable
◦
C
ρ
a
1.341
Kg/m
3
ρ
i
917
Kg/m
3
ρ
w
1000
Kg/m
3
c
p,a
1
.
044
×
10
3
J/(Kg
◦
C)
c
p,i
1
.
96
×
10
3
J/(Kg
◦
C)
c
p,w
4
.
2
×
10
3
J/(Kg
◦
C)
K
a
0.02
W/(m
◦
C)
K
i
2.29
W/(m
◦
C)
K
w
0.554
W/(m
◦
C)
L
sol
3
.
34
×
10
5
J/Kg
L
sub
2
.
83
×
10
6
J/Kg
D
v
0
2
.
178
×
10
−
5
m
2
/s
ξ
v
10
−
3
(-)
ξ
T
see Eq. (13)
(-)
9
2.4 Initial and boundary conditions
137
For the phase fields and the temperature we consider the boundary conditions
138
∇
φ
i
·
n
= 0
,
∇
φ
a
·
n
= 0
, T
=
T
B
(
x
,t
) on
∂
Ω
,
(14)
where
∂
Ω is the boundary of Ω and the boundary temperature
T
B
may take different values
139
on each boundary and may evolve in time. For the vapor concentration, we consider a fixed
140
concentration on the boundary or we impose no-flux through the boundary, depending on
141
each simulation, which are expressed as
142
ρ
v
=
ρ
v,B
or
∇
ρ
v
·
n
= 0 on
∂
Ω
,
(15)
where
ρ
v,B
is a constant value. The initial conditions of our problem may be written as
143
φ
l
(
x
,
0) = 0
.
5
1 + tanh
d
l
(
x
)
2
ε
, T
(
x
,
0) =
T
0
(
x
)
, ρ
v
(
x
,
0) =
ρ
I
vs
(
T
0
)
,
(16)
where
l
=
{
i,a
}
,
d
l
is the signed distance to the interface of phase
φ
l
and
T
0
(
x
) is a linear
144
function consistent with the temperature boundary condition
T
B
. The distance
d
l
depends
145
on the geometry considered in each simulation.
146
2.5 Spatial, time discretization and variable regularization
147
We use the Finite Element Method to solve our model. In particular, we use a uniform
148
mesh composed of bilinear (linear in 1D) basis functions. We compute the weak form by
149
multiplying the model equations (Eqs. (23)–(26) in the main text) by weighting functions,
150
integrating over Ω, and integrating by parts taking into account the boundary conditions
151
defined in Eqs. (14) and (15). To obtain the Galerkin form, we substitute the problem
152
unknowns and the weighting functions by the basis functions. The phase field variables
153
may take values
φ
l
<
0 and
φ
l
>
1 during transient states, with
l
=
{
i,w,a
}
. To avoid
154
10
singularities and numerical issues we regularize the functions
ρ
,
c
p
,
K
, and
ρ
SE
: We restrict
155
the value of
φ
l
between 0 and 1 only when we compute the functions
ρ
,
c
p
, and
K
. For
156
the regularization of
ρ
SE
we (1) implement an additional cubic interpolation locally around
157
φ
a
= 0 which smooths out the discontinuity at that point, and (2) extend
ρ
SE
off the ternary
158
diagram constant in the normal direction to the diagram sides.
159
References
160
(1) Karma, A.; Rappel, W.-J. Phase-field method for computationally efficient modeling
161
of solidification with arbitrary interface kinetics.
Physical review E
1996
,
53
, R3017.
162
(2) Gomez, H.; Bures, M.; Moure, A. A review on computational modelling of phase-
163
transition problems.
Philosophical Transactions of the Royal Society A
2019
,
377
,
164
20180203.
165
(3) Karma, A.; Rappel, W.-J. Quantitative phase-field modeling of dendritic growth in two
166
and three dimensions.
Physical review E
1998
,
57
, 4323.
167
(4) Karma, A. Phase-Field formulation for quantitative modeling of alloy solidification.
168
Phys. Rev. Lett.
2001
,
87
, 115701.
169
(5) T ́oth, G. I.; Pusztai, T.; Gr ́an ́asy, L. Consistent multiphase-field theory for interface
170
driven multidomain dynamics.
Phys. Rev. B
2015
,
92
, 184105.
171
(6) Kaempfer, T. U.; Plapp, M. Phase-field modeling of dry snow metamorphism.
Physical
172
Review E
2009
,
79
, 031502.
173
(7) Eames, I.; Marr, N.; Sabir, H. The evaporation coefficient of water: a review.
Interna-
174
tional Journal of Heat and Mass Transfer
1997
,
40
, 2963–2973.
175
(8) Libbrecht, K. G. Physical dynamics of ice crystal growth.
Annual Review of Materials
176
Research
2017
,
47
, 271–295.
177
11
(9) Krishnamachari, B.; McLean, J.; Cooper, B.; Sethna, J. Gibbs-Thomson formula for
178
small island sizes: Corrections for high vapor densities.
Physical Review B
1996
,
54
,
179
8899.
180
(10) Flatau, P. J.; Walko, R. L.; Cotton, W. R. Polynomial fits to saturation vapor pressure.
181
Journal of Applied Meteorology
1992
, 1507–1513.
182
(11) Wexler, A. Vapor pressure formulation for ice.
Journal of Research of the National
183
Bureau of Standards. Section A, Physics and Chemistry
1977
,
81
, 5.
184
12