# A phase-field model for wet snow metamorphism

Adrian Moure<EMAIL_ADDRESS>Xiaojing Fu<EMAIL_ADDRESS>Department of
Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, CA, United States

###### Abstract

The microstructure of snow determines its fundamental properties such as the
mechanical strength, reflectivity, or the thermo-hydraulic properties. Snow
undergoes continuous microstructural changes due to local gradients in
temperature, humidity or curvature, in a process known as snow metamorphism.
In this work, we focus on wet snow metamorphism, which occurs when temperature
is close to the melting point and involves phase transitions amongst liquid
water, water vapor, and solid ice. We propose a pore-scale phase-field model
that simultaneously captures the three relevant phase-change phenomena:
sublimation (deposition), evaporation (condensation), and melting
(solidification). The phase-field formulation allows one to track the
temperature evolution amongst the three phases and the water vapor
concentration in the air. Our three-phase model recovers the corresponding
two-phase transition model when one phase is not present in the system. 2D
simulations of the model unveils the impact of humidity and temperature on the
dynamics of wet snow metamorphism at the pore scale. We also explore the role
of liquid melt content in controlling the dynamics of snow metamorphism in
contrast to the dry regime, before percolation onsets. The model can be
readily extended to incorporate two-phase flow and may be the basis for
investigating other problems involving water phase transitions in a vapor-
solid-liquid system such as airplane icing or thermal spray coating.

###### keywords:

wet snow, firn, melt film, metamorphism

IR,NMR,UV

##  1 Introduction

Snow and firn are heterogenous porous material composed of two components,
water and air, distributed amongst three phases: solid ice, liquid water, and
air with water vapor. The microstructure of snow and firn dictates their
mechanical strength, reflectivity, and thermo-hydraulic properties, which
govern important processes such as snow avalanches, snowpack and glacial
hydrology, radar remote sensing, and the performance of snow vehicles 1.
However, these properties constantly evolve because snow is a
thermodynamically active material that undergoes continuous microstructural
changes caused by phase transitions between ice, liquid water, and vapor, in a
process known as snow metamorphism 2. Ultimately, snow metamorphism leads to a
denser snow composed of coarser and rounder ice grains. Based on whether
liquid water is present, snow metamorphism is considered in two types —dry and
wet. Colbeck and other authors published a series of pioneer papers in the
70’s and 80’s 1, 3, 2, 4, 5, 6, 7, 8 that derived the theoretical foundations
to understand the mechanisms driving each type of snow metamorphism. Dry snow
metamorphism 3, 9 results from sublimation and deposition and is driven by the
transport of water vapor between regions with different vapor pressure, i.e.,
regions with different temperature and/or grain curvature, according to the
Gibbs-Thomson condition 10. On the other side, wet snow metamorphism occurs
when the temperature is close to the melting point. The presence of liquid
water fundamentally alters the mechanisms of metamorphism by increasing the
thermodynamic activity of snow and strengthening the mechanical and thermal
connection between ice grains through capillary bridging 1, 3. Thus, compared
to dry metamorphism, wet snow metamorphism displays accelerated coarsening and
densification. Colbeck proposed two different regimes for wet metamorphism
depending on the liquid water content: the pendular and the funicular regimes
1. In the funicular regime, the liquid water content is high so that the water
completely surrounds the ice grains. In the pendular regime, the liquid water
content is low and most of the ice surface is in contact with air. Thus, the
pendular regime displays lower heat flow, grain growth, and snow
densification, and larger capillary forces and grain-to-grain mechanical
strength, compared to the funicular regime. In both cases, vapor kinetics
plays an important role 2. While Colbeck’s work provided mechanistic
understanding of snow metamorphism, it focused on theoretical snow
configurations under ideal environmental conditions. Subsequent experimental
works used X-ray microtomography (microCT) to study the evolution of the dry
snow microstructure in realistic scenarios under isothermal conditions 11, 12
and temperature gradients 13, 14. Many authors leveraged the microCT data to
perform a comprehensive analysis of local crystal growth laws 15, 16 and
validate mean field models 17, 18. The first experiments of wet metamorphism
focused on fully saturated snow 5. Improved experimental techniques enable
studies of wet metamorphism on samples with low liquid water content 19, 20,
21, 22. These works analyzed snow wetness and grain size and derived empirical
laws for the snow properties under wet conditions 23. However, more advanced
techniques are still being developed to robustly visualize liquid water within
snow samples, such as the use of hyperspectral imaging to map out liquid water
content in snow 24.

Mathematical models may complement experimental development to provide
insights into the wet snow metamorphism process. To this end, wet snow is
modeled as a three-phase system composed of liquid water, vapor and solid ice,
and thus can be grouped with other class problems such as snowflake growth and
water drop icing in a humidity-controlled atmosphere. In this broader group of
literature, some have proposed simplified 1D model to describe the role of
humidity on the evolution dynamics 25. More elaborated 2D continuum models
resorted to the phase-field method 26, 27 or the level-set method 28, 29,
which have been successful in reproducing the phase transitions between two
phases of water, i.e., solidification 30, 31, 32, sublimation 33, 34, 35, and
evaporation 36, 37, 38, 39. However, most of these continuum models are
limited to the two-phase regime, and do not simultaneously capture the
transitions amongst the three water phases; thus, these models do not address
the wet snow metamorphism problem.

Here, we propose a phase-field model that describes the non-equilibrium
evolution of the three-phase mixture composed of liquid water, water vapor and
solid ice. The modeling choice is inspired by recent applications of phase-
field methods in porous media problems with phase transitions, including the
modeling of microstructural evolution of porous materials and multicomponent
fluid mixtures 34, 40, 41, 42, 43, 44, phase transitions during two-phase flow
in porous media45, 46 and crystallization flow in porous media 47, 48, 49. The
three-phase model proposed in this work reduces to the classic two-phase model
in the absence of a third phase, and the resulting two-phase models readily
capture the Gibbs-Thomson conditions. The model allows us to analyze the pore-
scale dynamics of that dictates wet snow metamorphism. Currently, the model
does not account for granular compaction, melt film imbibition or melt
percolation. Rather, the primary focus of the model is to elucidate the role
of thermodynamic-driven phase transitions on microstructure evolution of a
quiescent wet snow mixture.

##  2 A phase-field model for wet snow mixture

The phase-field method 50, 51 is a mathematical technique that is based on a
free energy description of multi-phase mixtures, and is well-suited for
problems with moving interfaces. Within the phase-field framework, one defines
a dimensionless phase field variable, denoted $\phi_{l}\in[0,1]$, which
represents the volume fraction of the phase $l$ at any given point in the
domain. The phase variable smoothly transitions across the interface between
two distinct phases (Fig. 1). With the appropriate evolution equations, the
moving boundary problem is reformulated in a fixed domain, which avoids most
of the numerical issues caused by the moving interfaces.

Here, we propose a non-variational phase-field model that captures the phase
transitions between ice, liquid water, and water vapor, i.e., solidification,
sublimation, evaporation, and the opposite transitions. For simplicity, from
here on we denote the liquid water as water and the water vapor as vapor. We
assume that the air is at the atmospheric pressure. We assume that ice and
water have equal density for the phase evolution equations, which introduces
small mass conservation error. Although fluid flow is not considered in the
current model, future work that incorporates fluid flow and unequal density
may mitigate the mass conservation issue.

We note that in the spirit of numerical robustness and physical consistency, a
variational phase-field formulation 52 is often more desirable to ensure non-
negative entropy production and the convergence towards equilibrium 53, 54.
However, due to the close proximity in the density of the ice and liquid water
phase, we find that formulating a free energy for the wet snow system based on
the density field introduces numerical singularities in the formulation.

###  2.1 Model variables

Our model unknowns are the phase field variables for ice $\phi_{i}(\bm{x},t)$,
water $\phi_{w}(\bm{x},t)$, and air $\phi_{a}(\bm{x},t)$ phases (see Fig. 1),
as well as variables for the temperature $T(\bm{x},t)$ (units ∘C) and the
vapor density $\rho_{v}(\bm{x},t)$ (units kg/m3). All variables are defined
point-wise in the problem domain $\Omega$ (Fig. 1). Note that even though
temperature and vapor density are not phase field variables, here we leverage
the formulation to specify phase-dependent temperature evolution (see Section
2.2.2), and to localize the vapor dynamics onto the air phase (see Section
2.2.3).

Figure 1: Illustration of the key model variables
$\phi_{a},\phi_{w},\phi_{i},\text{ and }\rho_{v}$ as defined in the continuum
domain composed of wet snow. The 1D transect profile illustrates the diffusive
nature of the phase variables at the phase-phase boundaries. The model
variable temperature $T$ (not shown here) is also defined point-wise in the
entire domain.

###  2.2 Model equations

In this section, we will introduce the five partial differential equations
that describe the co-evolution of our model variables.

####  2.2.1 The phase evolution equations

The basic form of the phase evolution equations for the three phase field
variables are expressed as:

$\displaystyle\frac{\partial\phi_{i}}{\partial
t}=-\frac{M_{0}}{\Sigma_{i}}\frac{\delta\mathcal{F}^{\mathrm{tri}}(\overline{\phi})}{\delta\phi_{i}}-\alpha_{\mathrm{sol}}\phi_{i}^{2}\phi_{w}^{2}\frac{T-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}+\alpha_{\mathrm{sub}}\phi_{i}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{i}}+\theta_{\mathrm{nucl}}(\overline{\phi},T),$
(1) $\displaystyle\frac{\partial\phi_{w}}{\partial
t}=-\frac{M_{0}}{\Sigma_{w}}\frac{\delta\mathcal{F}^{\mathrm{tri}}(\overline{\phi})}{\delta\phi_{w}}+\alpha_{\mathrm{sol}}\phi_{i}^{2}\phi_{w}^{2}\frac{T-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}+\alpha_{\mathrm{eva}}\phi_{w}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{W}(T)}{\rho_{w}}-\theta_{\mathrm{nucl}}(\overline{\phi},T),$
(2) $\displaystyle\frac{\partial\phi_{a}}{\partial
t}=-\frac{M_{0}}{\Sigma_{a}}\frac{\delta\mathcal{F}^{\mathrm{tri}}(\overline{\phi})}{\delta\phi_{a}}-\alpha_{\mathrm{sub}}\phi_{i}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{i}}-\alpha_{\mathrm{eva}}\phi_{w}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{W}(T)}{\rho_{w}},$
(3)

where $\displaystyle\delta\cdot/\delta\cdot$ denotes the variational
derivative and we use the notation
$\overline{\phi}=\\{\phi_{i},\phi_{w},\phi_{a}\\}$. Note that Eqs.(1)-(3)
account for total mass conservation by assuming equal density of ice and water
($\rho_{i}=\rho_{w}$). The first terms in the right-hand side of these three
equations capture the Allen-Cahn kinetics of phase-field variables 55, where
$M_{0}$ is the phase-dependent mobility defined below in Sect. 2.3 and
$\mathcal{F}^{\mathrm{tri}}$ is the classical energy functional of a ternary
mixture 51:

$\mathcal{F}^{\mathrm{tri}}(\overline{\phi})=\int\left(\frac{3}{\varepsilon}F^{\mathrm{tri}}(\overline{\phi})+\frac{3}{2}\varepsilon\left[\Sigma_{i}(\nabla\phi_{i})^{2}+\Sigma_{w}(\nabla\phi_{w})^{2}+\Sigma_{a}(\nabla\phi_{a})^{2}\right]+\beta(1-\phi_{i}-\phi_{w}-\phi_{a})\right)\mathrm{\,d}\Omega.$
(4)

The parameter $\varepsilon$ represents the width of the phase-field interface
and $\beta$ is the Lagrange multiplier (see Section 3). The function
$F^{\mathrm{tri}}$ is the triple-well potential, defined as:

$F^{\mathrm{tri}}(\overline{\phi})=\frac{1}{2}\Sigma_{i}\phi_{i}^{2}(1-\phi_{i})^{2}+\frac{1}{2}\Sigma_{w}\phi_{w}^{2}(1-\phi_{w})^{2}+\frac{1}{2}\Sigma_{a}\phi_{a}^{2}(1-\phi_{a})^{2}+\Lambda\phi_{i}^{2}\phi_{w}^{2}\phi_{a}^{2},$
(5)

where the $\Lambda$-term improves the dynamic consistency of the model51 to
mitigate spurious effect along phase-phase boundaries. The parameters
$\Sigma_{i}$, $\Sigma_{w}$, and $\Sigma_{a}$ are related to the surface
tensions ($\sigma_{lm}$) between phases $l$ and $m$ such that

$\Sigma_{i}=\sigma_{iw}+\sigma_{ia}-\sigma_{wa},\quad\Sigma_{w}=\sigma_{iw}+\sigma_{wa}-\sigma_{ia},\quad\Sigma_{a}=\sigma_{wa}+\sigma_{ia}-\sigma_{iw}.$
(6)

We extend the Allen-Cahn kinetics with the $\alpha$-terms in Eqs. (1) to (3),
which account for the different phase transitions. In particular, the
$\alpha_{\text{sol}}$ term accounts for solidification or melting, the
$\alpha_{\text{sub}}$ term accounts for sublimation or deposition, and the
$\alpha_{\text{eva}}$ term accounts for evaporation or condensation. The
estimation of $\alpha$’s parameter values is described in the Supporting
Information (Section 1 and Table 1). The parameters $L_{\mathrm{sol}}$,
$c_{p,w}$, $\rho_{w}$ and $\rho_{i}$ are the solidification latent heat, water
specific heat capacity, water density, and ice density, respectively. We
follow Kaempfer and Plapp 34 and express the saturated vapor density of ice
$(\rho_{vs}^{I})$ and water $(\rho_{vs}^{W})$ as

$\rho_{vs}^{J}(T)=\rho_{a}0.62\frac{P_{vs}^{J}(T)}{P_{a}-P_{vs}^{J}(T)},$ (7)

where the index $J$ stands for $I$ and $W$, $\rho_{a}$ is the dry air density,
$P_{a}$ is the atmospheric pressure (units Pa), and $P_{vs}^{J}$ is the
saturated vapor pressure of phase $J$ (units Pa). The saturated vapor pressure
of ice and water are defined as 56, 57

$\displaystyle
P_{vs}^{I}(T)=\mathrm{exp}\left(\sum_{j=0}^{4}k_{j}(T+T_{0})^{j-1}+k_{5}\,\mathrm{ln}(T+T_{0})\right),$
(8) $\displaystyle
P_{vs}^{W}(T)=\mathrm{exp}\left(\sum_{j=0}^{6}g_{j}(T+T_{0})^{j-2}+g_{7}\,\mathrm{ln}(T+T_{0})\right),$
(9)

respectively, where $T_{0}=273.15$\text{\,}\mathrm{\SIUnitSymbolCelsius}$$,
and the values of the parameters $k_{j}$ and $g_{j}$ are listed in Table 2 in
the Supporting Information.

The $\theta_{\mathrm{nucl}}$ terms in Eqs. (1)–(2) account for ice and water
nucleation, which are necessary in order to trigger the nucleation of a third
phase when only two phases exist initially. In particular, here we account for
the nucleation of water when $T$ increases above $T_{\mathrm{melt}}$ in an
ice-air system and the nucleation of ice when $T$ decreases below the ice
nucleation temperature $T_{\mathrm{nucl}}$ in a water-air system. We propose a
phenomenological nucleation term

$\theta_{\mathrm{nucl}}(\overline{\phi},T)=\alpha_{\mathrm{nucl}}M_{0}\phi_{a}^{2}[\phi_{w}^{2}N_{i}(T)-\phi_{i}^{2}N_{w}(T)],$
(10)

where $\alpha_{\mathrm{nucl}}$ is a parameter controlling the strength of the
nucleation and the ice and water nucleation functions $N_{i}(T)$ and
$N_{w}(T)$, respectively, are expressed as

$N_{i}(T)=0.5-0.5\,\mathrm{tanh}[20(T-T_{\mathrm{nucl}}+0.1)],\quad
N_{w}(T)=0.5+0.5\,\mathrm{tanh}[20(T-T_{\mathrm{melt}}-0.1)].$ (11)

Note that $N_{i}$ becomes active when $T<T_{\mathrm{nucl}}$, which allows the
presence of undercooled water in the range
$T_{\mathrm{nucl}}<T<T_{\mathrm{melt}}$. According to Eq. (10), ice nucleation
occurs on the water-air interface when $T<T_{\mathrm{nucl}}$ (note the term
$\phi_{a}^{2}\phi_{w}^{2}$ multiplying $N_{i}$). Once ice nucleation advances,
the water-air interface splits into the water-ice and ice-air interfaces,
resulting in $\phi_{a}^{2}\phi_{w}^{2}=0$ and thus $\theta_{\mathrm{nucl}}$
stays zero even if $T<T_{\mathrm{nucl}}$. Water nucleation exhibits an
analogous behavior. Note that $\theta_{\mathrm{nucl}}$ only needs to be active
until the system enters into the spinodal decomposition region of the ternary
mixture.

####  2.2.2 Energy conservation equation

Here, we formulate conservation of energy in terms of temperature, expressed
as

$\rho(\overline{\phi})c_{p}(\overline{\phi})\frac{\partial T}{\partial
t}=\nabla\cdot\left[K(\overline{\phi})\nabla
T\right]+\rho(\overline{\phi})L_{\mathrm{sol}}\left(\frac{\partial\phi_{a}}{\partial
t}+\frac{\partial\phi_{i}}{\partial
t}\right)-\rho(\overline{\phi})L_{\mathrm{sub}}\frac{\partial\phi_{a}}{\partial
t}.$ (12)

Eq. (12) accounts for thermal diffusion and the latent heat released and/or
absorbed during the different phase transitions. $L_{\mathrm{sub}}$ is the
sublimation latent heat. The latent heat of evaporation is then captured
through $L_{\mathrm{eva}}=L_{\mathrm{sub}}-L_{\mathrm{sol}}$. The phase-
dependent density ($\rho$), specific heat capacity ($c_{p}$), and thermal
conductivity ($K$) are defined as:

$\rho(\overline{\phi})=\sum_{l=i,w,a}\phi_{l}\rho_{l};\;\;\;c_{p}(\overline{\phi})=\sum_{l=i,w,a}\phi_{l}c_{p,l};\;\;\;K(\overline{\phi})=\sum_{l=i,w,a}\phi_{l}K_{l},$
(13)

where $\rho_{l}$, $c_{p,l}$, and $K_{l}$ are the density, specific heat
capacity, and thermal conductivity, respectively, of phase $l$. We assume that
$\rho_{l}$, $c_{p,l}$, and $K_{l}$ are constant in time. Note that here we
violate our primary assumption that $\rho_{i}=\rho_{w}$ used in Eqs. (1)-(3),
which incurs mass conservation error; instead, we consider different density
for the ice and water phase for the energy equation in order to adopt the
correct heat capacity for each phase. Although doing so does not eliminate
errors in energy conservation due to errors in mass conservation, it improves
the thermal consistency of the model.

####  2.2.3 Mass conservation of vapor

The vapor mass conservation equation is expressed as

$\frac{\partial(\phi_{a}\rho_{v})}{\partial
t}=\nabla\cdot\left[\phi_{a}D_{v}(T)\nabla\rho_{v}\right]+\rho_{\mathrm{SE}}\frac{\partial\phi_{a}}{\partial
t},$ (14)

which accounts for vapor diffusion in the air and mass transfer of water
during sublimation and evaporation (and the opposite transitions). Here, we
formulate a phase-dependent density $\rho_{\mathrm{SE}}(\overline{\phi})$ such
that $\rho_{\mathrm{SE}}=\rho_{i}$ in case of sublimation and
$\rho_{\mathrm{SE}}=\rho_{w}$ in case of evaporation. The specific form for
$\rho_{\mathrm{SE}}(\overline{\phi})$ is discussed and provided in Eq. (17) of
Section 2.3. $D_{v}$ is the vapor diffusion coefficient, whose dependence on
$T$ (in $\text{\,}\mathrm{\SIUnitSymbolCelsius}$) follows the expression: 58

$D_{v}(T)=D_{v0}\left(\frac{T+T_{0}}{T_{0}}\right)^{1.81},$ (15)

where $D_{v0}$ is the vapor diffusion coefficient at the freezing point and
$T_{0}=273.15$.

###  2.3 Phase-dependent functions

The model equations introduced above rely on a few parameters that are phase-
dependent. For material property parameters such as $\rho,c_{p},K$, we simply
use phase-weighted averages as shown in Eq. (13). However, parameters such as
the mobility function $M_{0}$ in Eqs. (1)–(3) and $\rho_{\mathrm{SE}}$ in Eq.
(14) are related to phase change dynamics at phase-phase boundaries or triple
junction regions, and thus require more careful mathematical formulation. Here
we resort to the ternary diagram to illustrate these modeling choices.

A ternary diagram is the common tool to represent three-phase systems. In a
ternary diagram, the vertex $l$ represents the point $\phi_{l}=1$, while the
opposite side of that vertex represents the points $\phi_{l}=0$ (for
$l=\\{i,w,a\\}$). The interior of the diagram represents points where the
three phases coexist. Thus, any valid $\\{\phi_{i},\phi_{w},\phi_{a}\\}$
configuration is represented by a point in the ternary diagram. In particular,
evaporation takes place on the side $\phi_{i}=0$, sublimation on the side
$\phi_{w}=0$, and solidification on the side $\phi_{a}=0$.

At the boundaries of the ternary diagram (i.e. $\phi_{a}=0$, $\phi_{w}=0$, and
$\phi_{i}=0$), the mobility function $M_{0}$ must take the values
$M_{\mathrm{sol}}$, $M_{\mathrm{sub}}$, and $M_{\mathrm{eva}}$ that correspond
to the kinetics of solidification, sublimation and evaporation, respectively.
The exact expressions for these parameters are provided in Eqs. (36), (1.2),
and (1.3) of the Supporting Information. Within the interior region of the
ternary diagram, however, $M_{0}$ needs to be interpolated based on
$M_{\mathrm{sol}}$, $M_{\mathrm{sub}}$, and $M_{\mathrm{eva}}$. A common
choice is the cubic interpolation, which presents issues for two reasons59,
60: (1) the triple junction is a region, not a point, where the three phases
coexist; (2) the values of $M_{\mathrm{sol}}$, $M_{\mathrm{sub}}$, and
$M_{\mathrm{eva}}$ range several orders of magnitude in this region (see Table
1 in the Supporting Information). Thus, a cubic interpolation of $M_{0}$ would
result in large gradients in mobility within the interior of the ternary
diagram, leading to severe numerical challenges. In addition, a non-unique
value of $M_{0}$ within the triple junction is not physically consistent and
result in unrealistic behavior of the triple junction.

To avoid these two issues, we propose a piecewise constant interpolation of
$M_{0}$, similar to that proposed in Miyoshi and Takaki 59, where $M_{0}$
takes a constant value near each of the three sides and in the interior region
of the ternary diagram:

$M_{0}(\overline{\phi})=\left\\{\begin{array}[]{ll}M_{\mathrm{sol}}&\mathrm{if}\;\phi_{a}<=\phi_{i},\;\;\phi_{a}<=\phi_{w},\;\;\mathrm{and}\;\;\phi_{a}<0.01,\\\
M_{\mathrm{eva}}&\mathrm{if}\;\phi_{i}<\phi_{a},\;\;\phi_{i}<=\phi_{w},\;\;\mathrm{and}\;\;\phi_{i}<0.01,\\\
M_{\mathrm{sub}}&\mathrm{if}\;\phi_{w}<\phi_{a},\;\;\phi_{w}<\phi_{i},\;\;\mathrm{and}\;\;\phi_{w}<0.01,\\\
M_{\mathrm{av}}&\mathrm{otherwise},\end{array}\right.$ (16)

where we take $M_{\mathrm{av}}$ as the geometric mean (i.e.,
$M_{\mathrm{av}}=\sqrt[3]{M_{\mathrm{sub}}M_{\mathrm{sol}}M_{\mathrm{eva}}}$).
We find that Eq. (16), compared to a cubic interpolation, improves the
numerical efficiency and allows us to simulate triple junctions while
respecting the actual kinetics of the different phase transitions.

The density $\rho_{\mathrm{SE}}$ must take the value $\rho_{i}$ in case of
sublimation and $\rho_{w}$ in case of evaporation (i.e., $\phi_{w}=0$ and
$\phi_{i}=0$, respectively, in the ternary diagram). The value of
$\rho_{\mathrm{SE}}$ along the side $\phi_{a}=0$ is irrelevant because
$\rho_{\mathrm{SE}}\frac{\partial\phi_{a}}{\partial t}$ (see Eq. (14)) is zero
on that side. We propose a cubic interpolation between the sides $\phi_{w}=0$
and $\phi_{i}=0$ in the ternary diagram, such that the derivatives of
$\rho_{\mathrm{SE}}(\overline{\phi})$ are zero on the sides $\phi_{w}=0$ and
$\phi_{i}=0$. We define the function $\rho_{\mathrm{SE}}$ as

$\rho_{\mathrm{SE}}(\overline{\phi})=\left\\{\begin{array}[]{ll}\frac{\rho_{i}+\rho_{w}}{2}+\frac{\rho_{w}-\rho_{i}}{4}x_{\mathrm{SE}}(3-x_{\mathrm{SE}}^{2})&\mathrm{if}\;\phi_{a}\neq
1,\\\
\frac{\rho_{i}+\rho_{w}}{2}&\mathrm{if}\;\phi_{a}=1,\end{array}\right.\quad\mathrm{with}\quad
x_{\mathrm{SE}}=\frac{1-\phi_{a}-2\phi_{i}}{1-\phi_{a}}.$ (17)

Equation (17) displays a discontinuity at $\phi_{a}=1$. From a physical and
numerical point of view, this discontinuity does not represent an issue
because the term $\rho_{\mathrm{SE}}\frac{\partial\phi_{a}}{\partial t}$ is
zero if $\phi_{a}=1$. Nevertheless, we implement a regularization of Eq. (17)
(see Section 2.5 in the Supporting Information).

###  2.4 The Gibbs-Thomson condition

The Gibbs-Thomson effect describes the deviations in the equilibrium chemical
potential at a two-phase boundary due to the curvature of the interface. Such
effect readily applies to the three types of interface considered in wet snow
and influences both the kinetics and equilibrium conditions experienced by
curved interfaces. In Section 1 of the Supporting Information, we show that
our model is equivalent to the two-phase models for solidification,
sublimation, and evaporation when only two phases are present, and the
resulting reduced model equations recover the Gibbs-Thomson for two-phase
boundaries under certain conditions. In particular, our model reproduces,
under certain parameter regimes, the kinetics defined by the Gibbs-Thomson
condition for ice-liquid interfaces:

$\frac{T-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}=-d_{\mathrm{sol}}\chi-\beta_{\mathrm{sol}}v_{n},$
(18)

where $T$ is the interface temperature, $\beta_{\mathrm{sol}}$ is the kinetic
attachment coefficient, $d_{\mathrm{sol}}$ is the capillary length, $\chi$ is
the curvature of the interface (positive for spherical ice grains), and
$v_{n}$ is the normal velocity of the interface (positive for ice growth).
Note that, at thermodynamic equilibrium (i.e., $v_{n}=0$), Eq. (18) simplifies
to the classic Gibbs-Thomson equation:

$T^{\mathrm{equil}}=T_{\mathrm{melt}}-\frac{L_{\mathrm{sol}}d_{\mathrm{sol}}}{c_{p,w}}\chi,$
(19)

which states that the actual equilibrium temperature for solidification is
lower than $T_{\mathrm{melt}}=0^{\circ}$C for spherical ice grains. Similarly,
our model captures, under certain parameter regimes, the kinetics defined by
the Gibbs-Thomson condition for air-ice interfaces:

$\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{vs}^{I}(T)}=d_{\mathrm{sub}}\chi+\beta_{\mathrm{sub}}v_{n},$
(20)

and for air-liquid interfaces:

$\frac{\rho_{v}-\rho_{vs}^{W}(T)}{\rho_{vs}^{W}(T)}=d_{\mathrm{eva}}\chi+\beta_{\mathrm{eva}}v_{n},$
(21)

where $\beta_{\mathrm{sub}}$ and $d_{\mathrm{sub}}$ are the Gibbs-Thomson
coefficients for sublimation and $\beta_{\mathrm{eva}}$ and $d_{\mathrm{eva}}$
are the Gibbs-Thomson coefficients for evaporation. The model’s versatility in
capturing the Gibbs-Thomson conditions for all three types of two-phase
boundaries rely on the phase-dependent formulation of $M_{0}$ as discussed in
Section 2.3.

Note that because we do not assume equal $\Sigma_{\alpha}^{\prime}s$, spurious
phase (nonphysical third phase) may appear along phase-phase boundaries,
leading to inaccuracies in the equilibrium interface solution. We mitigate
this issue using the $\Lambda$-term in the free energy formulation (Eq. (5)).
In addition, we note that the Gibbs-Thomson condition is only recovered in
phase-field models under certain parameter conditions 34 or mathematical
constructions through asymptotic analysis 61, 62, 53. Our current formulation
does not entail the mathematical additions needed to recover the Gibbs-Thomson
condition for all interfaces for all parameter regimes; however, under certain
parameter regimes, the model recovers these conditions. More details are
provided in Section 3.2 and Supplementary Information.

##  3 Numerical implementation

We first apply the phase constraint $\phi_{w}=1-\phi_{i}-\phi_{a}$ to Eqs.
(1)–(3), which allows us to compute the Lagrange multiplier $\beta$ in Eq. (4)
and eliminate one phase evolution equation:

$\beta=\frac{1}{\frac{1}{\Sigma_{a}}+\frac{1}{\Sigma_{w}}+\frac{1}{\Sigma_{i}}}\frac{3}{\varepsilon}\left(\frac{1}{\Sigma_{i}}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{i}}+\frac{1}{\Sigma_{w}}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{w}}+\frac{1}{\Sigma_{a}}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{a}}\right)$ (22)

The resulting system of four coupled partial differential equations are:

$\displaystyle\frac{\partial\phi_{i}}{\partial t}$ $\displaystyle\quad=$
$\displaystyle-\frac{3M_{0}(\overline{\phi})}{\varepsilon\Sigma_{T}}\left[(\Sigma_{w}+\Sigma_{a})\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{i}}-\Sigma_{a}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{w}}-\Sigma_{w}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{a}}\right]+3M_{0}(\overline{\phi})\varepsilon\nabla^{2}\phi_{i}$
$\displaystyle-\alpha_{\mathrm{sol}}\phi_{i}^{2}\phi_{w}^{2}\frac{T-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}+\alpha_{\mathrm{sub}}\phi_{i}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{i}}$
$\displaystyle+\alpha_{\mathrm{nucl}}M_{0}(\overline{\phi})\phi_{a}^{2}[\phi_{w}^{2}N_{i}(T)-\phi_{i}^{2}N_{w}(T)],$
(23) $\displaystyle\frac{\partial\phi_{a}}{\partial t}$ $\displaystyle\quad=$
$\displaystyle-\frac{3M_{0}(\overline{\phi})}{\varepsilon\Sigma_{T}}\left[-\Sigma_{w}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{i}}-\Sigma_{i}\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{w}}+(\Sigma_{i}+\Sigma_{w})\frac{\partial
F^{\mathrm{tri}}}{\partial\phi_{a}}\right]+3M_{0}(\overline{\phi})\varepsilon\nabla^{2}\phi_{a}$
$\displaystyle-\alpha_{\mathrm{sub}}\phi_{i}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{i}}-\alpha_{\mathrm{eva}}\phi_{w}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{W}(T)}{\rho_{w}},$
(24)
$\displaystyle\frac{1}{\xi_{T}}\rho(\overline{\phi})c_{p}(\overline{\phi})\frac{\partial
T}{\partial t}$ $\displaystyle\quad=$
$\displaystyle\nabla\cdot\left[K(\overline{\phi})\nabla
T\right]+\rho(\overline{\phi})L_{\mathrm{sol}}\left(\frac{\partial\phi_{a}}{\partial
t}+\frac{\partial\phi_{i}}{\partial
t}\right)-\rho(\overline{\phi})L_{\mathrm{sub}}\frac{\partial\phi_{a}}{\partial
t},$ (25)
$\displaystyle\frac{1}{\xi_{v}}\frac{\partial(\phi_{a}\rho_{v})}{\partial t}$
$\displaystyle\quad=$
$\displaystyle\nabla\cdot\left[\phi_{a}D_{v}(T)\nabla\rho_{v}\right]+\rho_{\mathrm{SE}}(\overline{\phi})\frac{\partial\phi_{a}}{\partial
t},$ (26)

where $F^{\mathrm{tri}}$, $\rho_{\mathrm{SE}}$ and $M_{0}$ are defined in Eqs.
(5), (17), and (16), respectively. The saturated vapor densities
$\rho^{I}_{vs}$ and $\rho^{W}_{vs}$ are defined in Eqs. (7)–(9). The parameter
$\Sigma_{T}=\Sigma_{i}\Sigma_{w}+\Sigma_{i}\Sigma_{a}+\Sigma_{w}\Sigma_{a}$,
where $\Sigma_{l}$’s are defined in Eq. (6). The nucleation functions $N_{i}$
and $N_{w}$ are defined in Eq. (11). The phase-dependent functions $\rho$,
$c_{p}$, and $K$ are defined through Eq. (13) and $D_{v}$ follows Eq. (15).
The time scaling parameters $\xi_{T}$ and $\xi_{v}$ are introduced in the next
section (Section 3.1). The rest of the model parameters are listed in Tables 1
and 3 in the Supporting Information.

###  3.1 Temporal scaling

Wet snow metamorphism involves multiple coupled processes at different time
scales. Assuming a characteristic length scale of $L=10^{-4}\,$m, the
characteristic times for different processes can range from $\sim
10^{-4}\,\mathrm{s}$ for vapor diffusion to $\sim 10^{4}\,\mathrm{s}$ for
sublimation (see Section 2.1 in the Supporting Information for detailed
calculations). A monolithic numerical solver requires the use of time steps
that resolve the fastest process (i.e., vapor diffusion at $\sim
10^{-4}\,\mathrm{s}$ ), which would significantly increase the computational
cost. Here, we leverage the procedure explained in Kaempfer and Plapp 34 to
speed up the simulations. This approach leverages the fact that $T$ and
$\rho_{v}$ are quasi-steady compared to the phase transition kinetics. The
procedure consists of multiplying the right-hand side of the $T$ and
$\rho_{v}$ equations (Eqs. (12) and (14), respectively) by a time-scale factor
($\xi_{T}$ and $\xi_{v}$ respectively), so that the process of vapor and
thermal diffusion is numerically slowed down by orders of magnitude while
still maintaining a quasi-steady state in $T$ and $\rho_{v}$. This approach
enables the use of larger time steps without introducing noticeable errors in
the simulations. In this work, we take $\xi_{v}=10^{-3}$ and $\xi_{T}=1$ or
$\xi_{T}=10^{-2}$ depending on whether or not solidification occurs in any
parts of the domain. More details on this procedure are described in Section
2.1 of the Supporting Information.

###  3.2 Kinetic parameters constraints

In order to accurately capture the interface kinetics set by the Gibbs-Thomson
conditions, certain model parameters need to be constrained by the interface
physics. To this end, the relations between the Gibbs-Thomson parameters
$(\beta_{j},d_{j})$ and our model parameters $(M_{j},\alpha_{j})$ (for
$j=\\{\mathrm{sol},\mathrm{sub},\mathrm{eva}\\}$; see Section 1) can be
formulated from asymptotic analysis 32, 61 but are only valid under certain
conditions. In particular, we note that these asymptotic analysis do not
readily apply to the wet snow problem as they are originally developed for
compositional problems with symmetric 32, 61 and asymmetric diffusivities 62,
63. Solutions using anti-trapping current have been devised for the thermal
problems, but only with symmetric diffusivities 53. Thus, in order to robustly
capture the Gibbs-Thomson condition under all parameter regimes for the wet
snow problem, further mathematical additions remain to be developed.

We note, however, in the earlier and similar work on dry snow metamorphism by
Kaempfer and Plapp 34, it is shown that if we impose
$\varepsilon<D\beta^{\prime}$, where $D$ are the different thermal or vapor
diffusivities and $\beta^{\prime}$ is the scaled kinetic coefficient, then no
anti-trapping current is needed to correctly recover the sharp-interface
limit. In particular, for the wet snow problem, this corresponds to the
following sets of constraints:

$\displaystyle\varepsilon<\beta_{\text{sol}}\frac{K_{w}}{\rho_{w}c_{p,w}}\approx
10^{-5}\text{m};\;(\text{for water-ice interface})$ (27)
$\displaystyle\varepsilon<\beta_{\text{sub}}\frac{\rho_{vs}}{\rho_{i}}\frac{K_{a}}{\rho_{a}c_{p,a}}\approx
10^{-7}\text{m};\;(\text{for air-ice interface})$ (28)
$\displaystyle\varepsilon<\beta_{\text{eva}}\frac{\rho_{vs}}{\rho_{i}}\frac{K_{w}}{\rho_{w}c_{p,w}}\approx
10^{-8}\text{m}.\;(\text{for air-water interface})$ (29)

These conditions impose a restriction on $\varepsilon$, which represents the
diffusive interface width at phase-phase boundaries. The value of
$\varepsilon$ dictates the spatial discretization of our model, which must be
fine enough to resolve the interface width with sufficient numerical elements.
Given that the most restrictive conditions of the above requires
$\varepsilon<10^{-8}\,\mathrm{m}$, this severely restricts the domain size we
can simulate. In order to speed up the simulations, we consider a value
$\varepsilon<10^{-6}\,\mathrm{m}$, which allows us to use a coarser spatial
discretization. As a result, our simulation are less accurate in capturing
interface kinetics along the air-ice and air-water interfaces. However,
considering the phase change kinetics involving the air phase is slow compared
to the solidification/melting processes along the water-ice interface, we
suspect the overall error on wet snow processes is small.

###  3.3 Numerical methods

We perform both 1D and 2D simulations of our model equations using the finite
element method. In particular, we use a uniform mesh composed of bilinear
(linear in 1D) basis functions. For the time integration, we use a semi-
implicit algorithm based on the generalized-$\alpha$ method 64, 65, where we
treat implicitly all the terms except the function $M_{0}$ and the parameter
$\xi_{T}$, which are constant during each time step. We use the Newton-Raphson
method with adaptive time stepping scheme to solve the resulting nonlinear
system. To avoid singularities and numerical issues we regularize the
functions $\rho$, $c_{p}$, $K$, and $\rho_{\mathrm{SE}}$ (see Section 2.5 in
the Supporting Information).

We impose no-flux boundary conditions for the phase variables and Dirichlet
boundary condition for the temperature. The vapor boundary condition is either
Dirichlet or no-flux depending on specific case studies. In most simulations,
we initialize the vapor density with the saturated vapor density for ice at
the prescribed temperature.

##  4 Results

###  4.1 Directional solidification and nucleation of a third phase

In this section, we demonstrate the ability of our model to capture
directional solidification as imposed by a temperature gradient, and the
nucleation of a third phase in a system that is initially composed of two
phases only. In particular, we present three 1D simulations in Fig. 2 that
show (a) the complete freezing of water in an initial ice-water-air system
(Fig. 2A), (b) water nucleation in an initial ice-air system (Fig. 2B), and
(c) ice nucleation in an initial water-air system (Fig. 2C). To perform these
1D simulations, we impose fixed temperatures and no-flux boundary condition
for the vapor phase on the two ends of the domain, which measures
$1\text{\,}\mathrm{mm}$ in length and is discretized with 2000 elements.

Figure 2: Time evolution (from top to bottom) of (A) directional water
freezing towards the vapor phase, (B) water nucleation during melting of an
initially dry system, and (C) ice nucleation during an initially ice-free
system. The shaded regions in the main panels represent the different phases
(light blue, dark blue, and unshaded for ice, water, and air, respectively).
The red solid line indicates the temperature (left axis) and the dashed blue
line indicates the vapor density in the air (right axis). The insets show the
diffusive profiles of $\phi_{i}$ (light blue), $\phi_{w}$ (dark blue), and
$\phi_{a}$ (black) around the interface where phase change is occuring.

In the freezing simulation (Fig. 2A), the initially imposed temperature
gradient results in the directional solidification of water towards the air
phase. At each time step, the temperature distribution appears piecewise
linear, which confirms the quasi-steady behavior of $T$ as discussed in
Section 3.1. Although the $\rho_{v}$ dynamics are also quasi-steady, its
distribution is not always linear (see $t=27$ and $27.5\text{\,}\mathrm{s}$ in
Fig. 2A). This nonlinear response is primarily triggered when a new ice-air
interface is established upon complete freezing
($\sim$27.5\text{\,}\mathrm{s}$$), and a new equilibrium is being established
to satisfy the saturated vapor density in the presence of ice, prescribed by
Eq. (8).

The nucleation simulations demonstrate significantly more nonlinear dynamics.
In the water nucleation simulation (Fig. 2B), we impose
$0\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and
$6\text{\,}\mathrm{\SIUnitSymbolCelsius}$ on the left and right boundaries of
the domain respectively to induce melting at the initial ice-air interface. In
the ice nucleation simulation (Fig. 2C), the water-air system is initially at
$-2\text{\,}\mathrm{\SIUnitSymbolCelsius}$ uniformly and we consider the ice
nucleation temperature
$T_{\mathrm{nucl}}=$-1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. Ice nucleation
onsets, preferentially towards the water phase, well before the temperature on
the water-air interface increases to $T_{\mathrm{nucl}}$ (Fig. 2C,
$t=$0.2\text{\,}\mathrm{s}$$). At $t=$0.3\text{\,}\mathrm{s}$$, when
$T>T_{\mathrm{nucl}}$ on the interface, the ice phase has fully developed.
Both types of nucleation described here are highly energetic events, leading
to rapid changes in temperature and saturated vapor density at the nucleation
region. As a result, the $T$ and $\rho_{v}$ profiles are no longer linear, in
contrast to the directional freezing problem (Fig. 2A). We note that here,
nucleation only occurs at the phase-phase interface due to its mathematical
formulation in Eq. (10).

We note that the phase field variables maintain the tanh-profile (quasi-
equilibrium profile) during the freezing simulation (insets of Fig. 2A), but
such profiles are temporarily lost during the nucleation simulations (insets
Figs. 2B, C). In the latter cases, the loss of the tanh-profile is due to
changes in mobility $M_{0}$ (Eq. (16)) during the dynamic nucleation events,
and could result in small errors in curvature calculations and thus influence
the kinetics of phase change through the Gibbs-Thomson condition.
Nevertheless, the tanh-profile is recovered shortly after nucleation (e.g., at
$t=$1\text{\,}\mathrm{s}$$ in Figs. 2B, C); thus the undesired loss of the
tanh-profile is temporary (occurs only during nucleation) and for this reason,
we assume its impact on the overall kinetics of the problem is negligible.

###  4.2 Simulation of dry snow metamorphism in 2D

Under thermal conditions when liquid water phase is absent ($\phi_{w}$=0), our
model reduces to the two-phase model for dry snow metamorphism first proposed
in Kaempfer and Plapp 34:

$\displaystyle\frac{\partial\phi_{i}}{\partial t}$
$\displaystyle=-3M_{0}\left(\frac{1}{\varepsilon}\phi_{i}(1-\phi_{i})(1-2\phi_{i})-\varepsilon\nabla^{2}\phi_{i}\right)+\alpha_{\mathrm{sub}}\phi_{i}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{i}}$
(30)
$\displaystyle\frac{1}{\xi_{T}}\rho(\overline{\phi})c_{p}(\overline{\phi})\frac{\partial
T}{\partial t}$ $\displaystyle=\nabla\cdot\left[K(\overline{\phi})\nabla
T\right]-\rho(\overline{\phi})L_{\mathrm{sub}}\frac{\partial\phi_{a}}{\partial
t},$ (31)
$\displaystyle\frac{1}{\xi_{v}}\frac{\partial(\phi_{a}\rho_{v})}{\partial t}$
$\displaystyle=\nabla\cdot\left[\phi_{a}D_{v}(T)\nabla\rho_{v}\right]+\rho_{\mathrm{SE}}\frac{\partial\phi_{a}}{\partial
t},$ (32)

with $\partial\phi_{a}/\partial t=-\partial\phi_{i}/\partial t$. Here, we show
a 2D example of the dry metamorphism model using a simple snow geometry
composed of four circular ice grains in contact (see $t_{0}$ in Fig. 3A). The
domain measures $0.2\times$0.2\text{\,}\mathrm{m}\mathrm{m}^{2}$$ and is
meshed with $400\times 400$ elements. We impose a fixed temperature
$T_{B}=$-5\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ and no-flux of $\rho_{v}$
on all boundaries of the domain. Over the 90 minute time span of the
simulation, we plot the time evolution of the ice phase $\phi_{i}$ and vapor
concentration $\rho_{v}$ (Fig. 3A), along with the $T$ (Fig. 3B) and the total
ice interface length $S_{i}$ (Figs. 3C and D). Also known as coarsening, dry
snow metamorphism under isothermal condition is a curvature-driven process,
which minimizes the total interfacial length and the curvature of the system.
The process penalizes regions of high curvature and is driven by vapor mass
transfer between regions of different curvatures 3. Our simulation readily
illustrates the above mechanism. We observe that the ice grains quickly sinter
(before $t_{1}=$40\text{\,}\mathrm{s}$$) and the coarsening process gradually
slows down, as measured by the decay in slope in Figs. 3C and D. Accompanying
these geometric changes are the dynamics in $T$ distribution due to latent
heat generation in regions of higher curvature (see Fig. 3B). The simulation
also captures the entrapment of an air bubble during the coarsening process.
At steady state (e.g., $t=$3\text{\,}\mathrm{h}$$ in Fig. 3D), the vapor
concentration inside the bubble is set by its curvature according to the
Gibbs-Thomson condition (Eq. (20)).

Figure 3: Dry snow metamorphism. Time evolution of (A) the ice geometry and
vapor concentration and (B) the temperature distribution at
$t_{0}=$0\text{\,}\mathrm{s}$$, $t_{1}=$40\text{\,}\mathrm{s}$$,
$t_{2}=$100\text{\,}\mathrm{s}$$, and $t_{3}=$240\text{\,}\mathrm{s}$$. The
blue solid line in (B) is the isoline $\phi_{i}=0.5$. Time evolution of the
total ice perimeter $S_{i}$ for short (C) and longer (D) times. The ice
geometry at $t=$30\text{\,}\mathrm{min}$$, $1\text{\,}\mathrm{h}$, and
$3\text{\,}\mathrm{h}$ is included in (D).

###  4.3 Simulations of wet snow dynamics in 2D

In this section, we perform a suite of numerical simulations of our full model
to investigate wet snow dynamics as it undergoes phase transition. Due to high
computational cost of the model, here we only focus on small sized domains
($\sim$mm) in 2D.

####  4.3.1 Influence of vapor concentration

In this part, we study the impact of vapor concentration on the kinetics of
melt generation when the initially dry material is subjected to melting. To
achieve this, we consider the same four-grain geometry as the initial dry snow
simulation in Fig. 3, we impose a boundary temperature of
$T_{B}=$1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ and run two simulations with
different vapor concentrations on the boundary, denoted $\rho_{v,B}$. In
particular, we impose an undersaturated vapor density,
$\rho_{v,B}=0.8\rho^{I}_{vs}($1\text{\,}\mathrm{\SIUnitSymbolCelsius}$)$, in
the first simulation (Fig. 4A) and a saturated vapor density,
$\rho_{v,B}=\rho^{I}_{vs}($1\text{\,}\mathrm{\SIUnitSymbolCelsius}$)$, in the
second one (Fig. 4B).

We quantify these two simulations by plotting the time evolution of the
spatially averaged temperature $\overline{T}$, volume of ice phase $V_{i}$ and
water phase $V_{w}$ in Figs. 5A and C. The results illustrate that, when the
initial vapor phase is undersaturated, melt nucleation is delayed and the
kinetics of melt generation is slower (Fig. 5C). We explain this from the
perspectives of thermal balance. Because both sublimation and melt generation
requires energy input, an unsaturated vapor phase will reduce the amount of
energy available for melting. This is evident when we observe that the average
temperature decrease is stronger in the unsaturated vapor case (Fig. 5A,
dashed line). This lower temperature also causes a delay in water nucleation
(Figs. 4A, Figs. 5C). We also compare the final snapshots of the two
simulations, taken at different times ($t=$45\text{\,}\mathrm{s}$$ for the
first and $t=$16\text{\,}\mathrm{s}$$ for the second simulation) but both
correspond to a state when ice no longer exists (Fig. 5B). It is evident that
there is significantly less melt water in the first simulation, which is
caused by the significant sublimation that claims part of the ice in order to
saturate the vapor phase.

Figure 4: Impact of vapor concentration on melting of snow. Time evolution of
the ice and water phases and vapor concentration for boundary vapor
concentration of (A) $\rho_{v,B}=0.8\rho_{vs}^{I}$ and (B)
$\rho_{v,B}=\rho_{vs}^{I}$ when snow melts with
$T_{B}=$1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. The initial snow geometry
is plotted in Fig. 3.  Figure 5: More detailed analysis of the simulations
shown in Figure 4. (A) Time evolution of the spatial-average temperature
$(\overline{T})$ when $\rho_{v,B}=0.8\rho_{vs}^{I}$ (dashed line) and
$\rho_{v,B}=\rho_{vs}^{I}$ (solid line). (B) Ice and water phases at times
$t_{5}=$16\text{\,}\mathrm{s}$$ and $t_{6}=$45\text{\,}\mathrm{s}$$ when
$\rho_{v,B}=0.8\rho_{vs}^{I}$ (top) and at time
$t_{5}=$16\text{\,}\mathrm{s}$$ when $\rho_{v,B}=\rho_{vs}^{I}$ (bottom). (C)
Time evolution of the ice ($V_{i}$, left) and water ($V_{w}$, right) volume
when $\rho_{v,B}=0.8\rho_{vs}^{I}$ (dashed lines) and
$\rho_{v,B}=\rho_{vs}^{I}$ (solid lines).

Finally, we want to remark that while water nucleation occurs in many parts of
the ice interface, melting into the ice phase only advances in some regions of
the interface because the process is limited by the thermal energy influx
(Figs. 4 and 5). In particular, we observe that melting usually starts at the
grain contact regions because larger curvatures induce higher temperatures, as
already shown in the case with dry snow (Fig. 3B).

####  4.3.2 Influence of temperature

In this part, we study the impact of temperature on wet snow evolution. Using
the same initial dry snow geometry (Fig. 3), here we fix the vapor
concentration on the boundary such that
$\rho_{v,B}=\rho^{I}_{vs}($1\text{\,}\mathrm{\SIUnitSymbolCelsius}$)$ and
consider two different boundary temperature $T_{B}$’s, namely,
$T_{B}=$1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (Fig. 6A) and
$T_{B}=$2\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (Fig. 6B). We further
quantify these two simulations with the time evolution of $V_{i}$ and $V_{w}$
(Fig. 6C) and the ice-water interface length $S_{iw}$ (Fig. 6D). As expected,
we observe that melting is faster at a higher boundary temperature (Fig. 6C).
In particular, our results illustrate that a higher temperature accelerates
melting by promoting melt nucleation and thus increases the amount of ice-
water interface available for melt advancing (Fig. 6D). In particular, we
observe that a film of nucleated water completely surrounds the ice phase at
the higher temperature (Fig. 6B), while water film only appears at some
regions of the ice interface at the lower temperature(Fig. 6A). After
nucleation, we observe that melting only advances at the three grain contact
locations for low $T_{B}$, while melting progresses in a few more locations,
some away from grain contacts, at high $T_{B}$ (Fig. 6B).

Finally, upon complete melting, the two simulations display a small difference
in $V_{w}$ (Fig. 6C, at $t=$10\text{\,}\mathrm{s}$$). This difference is
caused by the imposed $T$ and $\rho_{v}$ boundary conditions. Since we impose
$\rho_{v,B}=\rho^{I}_{vs}($1\text{\,}\mathrm{\SIUnitSymbolCelsius}$)$ for both
simulations, vapor concentration on the boundary is out of equilibrium for the
case of $T_{B}=$2\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. This induces a
vapor outflux through the boundary which is reflected in final value of
$V_{w}$ via sublimation/evaporation.

Figure 6: Impact of temperature on melting of snow. (A,B) Time evolution of
the ice and water phases and vapor concentration (upper half of each panel)
and temperature distribution (bottom half of each panel) for boundary
temperature (A) $T_{b}=$1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ and (B)
$T_{b}=$2\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. The initial snow geometry
is plotted in Fig. 3. (C) Time evolution of the ice ($V_{i}$, light blue) and
water ($V_{w}$, dark blue) volume for
$T_{b}=$1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (solid lines) and
$T_{b}=$2\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (dashed lines). (D) Time
evolution of the ice-water interface length $S_{iw}$ for
$T_{b}=$1\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (solid line) and
$T_{b}=$2\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (dashed line).

####  4.3.3 Influence of liquid water content

Figure 7: Influence of LWC on quasi-isothermal snow metamorphism, where the
boundary temperature $T_{B}$ is fixed to $T_{\mathrm{melt}}$. (A) Time
evolution of the ice and water phases for initial LWC (LWC0) of 0 (top row),
0.15 (center row), and 0.46 (bottom row; fully wet). (B) Time evolution of the
normalized interface length of the air phase $(\overline{S}_{a})$ for
LWC${}_{0}=0$ (black) and LWC${}_{0}=0.15$ (red) and the ice phase
$(\overline{S}_{i})$ for LWC${}_{0}=0.46$ (blue). The interface length is
normalized by its initial value. (C) Time evolution of the ice volume
$(V_{i})$ for LWC${}_{0}=0$ (black), 0.15 (red), and 0.46 (blue).

In this section, we study the influence of the liquid water content (LWC) on
wet snow metamorphism. We consider a simplified snow geometry composed of
circular ice grains in a domain of
$0.2\times$0.2\text{\,}\mathrm{m}\mathrm{m}^{2}$$ meshed with $800\times 800$
elements (Fig. 7A, $t=0$s). We run three simulations with different initial
LWC (LWC0) that correspond to dry (LWC${}_{0}=0$), partially wet
(LWC${}_{0}=0.15$), and fully wet (LWC${}_{0}=0.46$) snow. The initial ice
phase distribution is the same for all the three cases. For the partially wet
snow, we impose a water film of $2.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$
thickness uniformly surrounding the ice grains. For the fully wet snow, we
fill the pore space entirely with the liquid phase. We adopt the parameter
value $\varepsilon=2\times 10^{-7}\,$m. To keep the system in quasi-isothermal
conditions, we impose a fixed temperature $T_{B}=T_{\mathrm{melt}}$ and no-
flux in terms of vapor density $\rho_{v}$ on all boundaries.

The simulations results (Figs. 7 and 8) illustrate the influence of LWC0 on
the detailed dynamics of metamorphism. In particular, when comparing against
dry snow, we find that partially wet snow evolves towards larger and fewer ice
clusters faster (e.g., LWC${}_{0}=0$ and LWC${}_{0}=0.15$ at
$t=$10\text{\,}\mathrm{s}$$ in Fig. 7A). While this accelerated coarsening in
partially wet snow is visually apparent, plot of the normalized interface
length $\overline{S}_{a}$ does not readily capture this difference (Fig. 7B,
black and red lines). Note that $\overline{S}_{a}$ in the partially wet snow
accounts for the interface of the combined phase ice and water
$(\phi_{i}+\phi_{w})$. Meanwhile, the coarsening speed of the fully wet snow
is much faster, as captured visually in Fig. 7A and also in terms of the
normalized ice interface length $\overline{S}_{i}$ (Fig. 7B, blue line).

Figure 8: Influence of LWC on quasi-isothermal snow metamorphism, where the
boundary temperature $T_{B}$ is fixed to $T_{\mathrm{melt}}$. Time evolution
of the temperature distribution for initial LWC (LWC0) of 0 (top row), 0.15
(center row), and 0.46 (bottom row; fully wet). Light and dark blue lines
represent the ice and water interfaces, respectively.

Wet snow also experiences a broader range of dynamic temperatures than dry
snow during metamorphism (Fig. 8). This may be explained by the fact that in
wet snow, the thermal fluctuations are dominated by the latent heat
released/absorbed during freezing/melting process, which is orders of
magnitude faster than the process of deposition/sublimation in dry snow (see
Table 1 in the Supporting Information). The temperature distribution also
reflects the Gibbs-Thomson effect: because the equilibrium temperature
($T^{\mathrm{equil}}$, Eq. (19)) decreases with increasing ice grain
curvature, smaller ice grains require colder temperature
($T^{\mathrm{equil}}<T_{\mathrm{melt}}$) to be stable (Fig. 8). Because we
impose a temperature $T_{\mathrm{melt}}$ at the boundaries, which is higher
than the equilibrium temperature for curved grains, we find that all ice
grains eventually melt away in the case of fully saturated snow (Figs. 7A, C,
LWC$=0.46$). Note that here, if we initialize the simulation with _concave_
grain shapes, this will lead to $T^{\mathrm{equil}}>T_{\mathrm{melt}}$ and
complete freezing of the water in the case of fully wet snow.

A key characteristic of dry snow metamorphism is the reduction in
microstructure surface area while maintaining an almost constant ice volume
(Fig. 7B-C, black curves). Interestingly, we find that for the partially wet
snow case investigated here, the volume of the air phase is also roughly
constant throughout the simulation (not shown here). We attribute this to the
fact that water storage capacity of the air is limited (the initial vapor
concentration is already close to the saturated vapor density). Nevertheless,
we do observe a small mass exchange between the ice and water phases for
partially wet snow, as evidenced by the changes in $V_{i}$ for
LWC${}_{0}=0.15$ in Fig. 7C.

##  5 Numerical challenges and model improvements

We have so far presented a series of 1D and 2D simulations of wet snow in
domains of limited sizes ($\sim 0.2$mm). The reason we do not explore larger
problems in 3D is due to the high computational cost associated with these
simulations. This limits our ability to compare with experimental studies,
which are always in 3D and done on samples that are $\sim$cm in size (e.g., in
Colbeck 5). Here we summarize these computational challenges as well as points
for model improvements in future work.

In order to accurately capture the kinetics defined by the Gibbs-Thomson
conditions, our model requires that the numerical interface width parameter
$\varepsilon<10^{-8}\,$m (see Section 3.2). The consequence of this
requirement is the need for a rather fine numerical mesh that resolves such
interface, making it computationally expensive to simulate larger problems. An
adaptive mesh refinement algorithm, rather than the uniform mesh used in this
work, would allow one to run mm- and cm-scale simulations in reasonable
computational times. Additionally, mathematical additions to our model based
on an asymptotic analysis for the thermal problem with asymmetric
diffusivities could further improve robustness of our model at the interface.

Another cause for high computational cost is the need to resolve the kinetics
of different processes. In this work, we resort to a monolithic scheme, where
our time step size is limited by the fastest process. We have alleviated this
issue partially using the time-scaling strategy explained in Section 3.1.
Another way to speed up time integration is to neglect the phase transitions
related to the vapor phase during melting/freezing scenarios. The errors
introduced with this approach would be small as long as melting/freezing
occurs for short time intervals compared to the total time of interest.

Finally, our model is currently not suited for studying the dynamics of the
triple junction —the region where the air, liquid and ice phases contact— when
the system is near the triple point. At the triple point, defined here as
$\\{T,\rho_{v}\\}=\\{T_{\mathrm{melt}},\rho^{I}_{vs}(T_{\mathrm{melt}})\\}$,
the triple junction is expected to display an equilibrium configuration
defined by the surface tensions in Eq. (6). In our model, an initially non-
equilibrium triple junction configuration will evolve towards equilibrium
through the co-movement of the water-air, air-ice and water-ice interfaces.
However, because each interface moves at different speed, as defined by the
phase-dependent mobility (Eq. (16)), a true equilibrium may not be possible.
Instead, the model produces oscillations around the triple junction that
impede the achievement of the equilibrium configuration 59, 60. Far away from
the triple point, phase change kinetics around the triple junction dominate
over the surface tension kinetics and these oscillations do not appear.

For future work, we will focus on more realistic representation of the
nucleation process that includes randomized Gaussian noise 66, 67, which would
produce a more realistic behavior of the wet snow metamorphism problem. The
current model also does not consider the density difference between the solid
ice and liquid water phase, which would introduce small mass conservation
error (the model is volume conserved). Our model might be expanded based on
Huang et al. 68, Zhang et al. 27, Hagiwara et al. 26 to account for mass
conservation during melting/freezing, but further research is needed.

##  6 Conclusions

We propose a non-variational phase-field model for wet snow metamorphism. The
model accounts for the ice, water, and air phases, along with the temperature
and vapor dynamics, and captures the actual kinetics of solidification,
sublimation, and evaporation (and the opposite transitions). To the best of
our knowledge, this is the first phase-field model that simultaneously
reproduces the different phase transitions of water. The model results unveil
the intimate coupling amongst the various transport and phase change processes
involved in wet snow metamorphism. Our results show that local humidity
conditions affect snow melting and melt refreezing. The results also reveal
the differences between dry and wet snow metamorphism when water flow is not
considered: partially and fully wet snow experience larger thermal
fluctuations due to the dominant freeze/melt processes and the coarsening rate
is higher for fully wet snow compared to dry snow as expected. However, our
observations are based on numerically generated partial wet snow where a
liquid film of constant thickness is imposed over all ice grain surfaces. We
recognize that this way of initializing the simulation may not represent
actual liquid film distribution in wet snow and thus may result in unrealistic
dynamics. In addition, because we do not consider thin film flow and other
capillary-driven processes of the liquid phase, our model may be missing key
fluid mechanical processes that influence how wet snow evolves.

The current results focus on 1D and 2D problems, due to the high computational
cost of the model equations. Future modification, extension and application of
this model in 3D may be used to quantitatively analyze the evolution of the
snow pore structure during wet snow metamorphism, which dictates the thermo-
mechanical and hydraulic properties of the snow at larger scales. The proposed
modeling framework has proven useful to study wet snow metamorphism and may be
the basis to investigate other mm- or cm-scale problems involving water phase
transitions such as, e.g., the freezing of a water droplet on a surface 69,
70, 71, 25, the human-induced thermo-mechanical changes of icy planetary
surfaces, or water spray cooling.

The authors acknowledge the partial support from the Resnick Sustainability
Institute at California Institute of Technology, the National Science
Foundation under Grant No. EAR-2243631 and the ACS Petroleum Research Fund
Doctoral New Investigator Grant No. 66867-DNI9. The authors also acknowledge
insightful discussions with Dr. Quirine Krol.

## References

  * Colbeck 1973 Colbeck, S. C. Theory of metamorphism of wet snow. 1973,
  * Colbeck 1996 Colbeck, S. C. The basic ideas behind snow metamorphism. Proceedings of the International Snow Science Workshop. 1996
  * Colbeck 1980 Colbeck, S. C. Thermodynamics of snow metamorphism due to variations in curvature. _Journal of Glaciology_ 1980, _26_ , 291–301
  * Colbeck 1982 Colbeck, S. C. An overview of seasonal snow metamorphism. _Reviews of Geophysics_ 1982, _20_ , 45–61
  * Colbeck 1986 Colbeck, S. C. Statistics of coarsening in water-saturated snow. _Acta Metallurgica_ 1986, _34_ , 347–352
  * Colbeck 1997 Colbeck, S. C. A Review of Sintering in Seasonal Snow. 1997,
  * Raymond and Tusima 1979 Raymond, C. F.; Tusima, K. Grain coarsening of water-saturated snow. _Journal of Glaciology_ 1979, _22_ , 83–105
  * Wakahama 1975 Wakahama, G. The role of meltwater in densification processes of snow and firn. _International Association of Hydrological Sciences Publication_ 1975, _114_ , 66–72
  * Gundlach et al. 2018 Gundlach, B.; Ratte, J.; Blum, J.; Oesert, J.; Gorb, S. Sintering and sublimation of micrometre-sized water-ice particles: The formation of surface crusts on icy Solar System bodies. _Monthly Notices of the Royal Astronomical Society_ 2018, _479_ , 5272–5287
  * Johnson 1965 Johnson, C. A. Generalization of the Gibbs-Thomson equation. _Surface Science_ 1965, _3_ , 429–444
  * Löwe et al. 2011 Löwe, H.; Spiegel, J.; Schneebeli, M. Interfacial and structural relaxations of snow under isothermal conditions. _Journal of Glaciology_ 2011, _57_ , 499–510
  * Kaempfer and Schneebeli 2007 Kaempfer, T. U.; Schneebeli, M. Observation of isothermal metamorphism of new snow and interpretation as a sintering process. _Journal of Geophysical Research: Atmospheres_ 2007, _112_
  * Pinzer et al. 2012 Pinzer, B.; Schneebeli, M.; Kaempfer, T. Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography. _The Cryosphere Discussions_ 2012, _6_ , 1673–1714
  * Wang and Baker 2014 Wang, X.; Baker, I. Evolution of the specific surface area of snow during high-temperature gradient metamorphism. _Journal of Geophysical Research: Atmospheres_ 2014, _119_ , 13–690
  * Fife et al. 2014 Fife, J.; Gibbs, J.; Gulsoy, E.; Park, C.-L.; Thornton, K.; Voorhees, P. The dynamics of interfaces during coarsening in solid–liquid systems. _Acta Materialia_ 2014, _70_ , 66–78
  * Krol and Löwe 2016 Krol, Q.; Löwe, H. Analysis of local ice crystal growth in snow. _Journal of Glaciology_ 2016, _62_ , 378–390
  * Legagneux and Domine 2005 Legagneux, L.; Domine, F. A mean field model of the decrease of the specific surface area of dry snow during isothermal metamorphism. _Journal of Geophysical Research: Earth Surface_ 2005, _110_
  * Krol and Löwe 2018 Krol, Q.; Löwe, H. Upscaling ice crystal growth dynamics in snow: Rigorous modeling and comparison to 4D X-ray tomography data. _Acta Materialia_ 2018, _151_ , 478–487
  * Brun 1989 Brun, E. Investigation on wet-snow metamorphism in respect of liquid-water content. _Annals of Glaciology_ 1989, _13_ , 22–26
  * Ebner et al. 2020 Ebner, P. P.; Coulin, A.; Borner, J.; Wolfsperger, F.; Hohl, M.; Schneebeli, M. Liquid-water content and water distribution of wet snow using electrical monitoring. _The Cryosphere Discussions_ 2020, 1–34
  * Stein et al. 1997 Stein, J.; Laberge, G.; Lévesque, D. Monitoring the dry density and the liquid water content of snow using time domain reflectometry (TDR). _Cold Regions Science and Technology_ 1997, _25_ , 123–136
  * Coléou and Lesaffre 1998 Coléou, C.; Lesaffre, B. Irreducible water saturation in snow: experimental results in a cold laboratory. _Annals of glaciology_ 1998, _26_ , 64–68
  * Brun et al. 1992 Brun, E.; David, P.; Sudul, M.; Brunot, G. A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting. _Journal of Glaciology_ 1992, _38_ , 13–22
  * Donahue et al. 2022 Donahue, C.; Skiles, S. M.; Hammonds, K. Mapping liquid water content in snow at the millimeter scale: an intercomparison of mixed-phase optical property models using hyperspectral imaging and in situ measurements. _The Cryosphere_ 2022, _16_ , 43–59
  * Sebilleau et al. 2021 Sebilleau, J.; Ablonet, E.; Tordjeman, P.; Legendre, D. Air humidity effects on water-drop icing. _Physical Review E_ 2021, _104_ , L032802
  * Hagiwara et al. 2017 Hagiwara, Y.; Ishikawa, S.; Kimura, R.; Toyohara, K. Ice growth and interface oscillation of water droplets impinged on a cooling surface. _Journal of Crystal Growth_ 2017, _468_ , 46–53
  * Zhang et al. 2022 Zhang, W.; Shahmardi, A.; Choi, K.-s.; Tammisola, O.; Brandt, L.; Mao, X. A phase-field method for three-phase flows with icing. _Journal of Computational Physics_ 2022, _458_ , 111104
  * Vahab et al. 2016 Vahab, M.; Pei, C.; Hussaini, M. Y.; Sussman, M.; Lian, Y. An adaptive coupled level set and moment-of-fluid method for simulating droplet impact and solidification on solid surfaces with application to aircraft icing. 54th AIAA Aerospace Sciences Meeting. 2016; p 1340
  * Vahab et al. 2022 Vahab, M.; Murphy, D.; Shoele, K. Fluid dynamics of frozen precipitation at the air–water interface. _Journal of Fluid Mechanics_ 2022, _933_ , A36
  * Cristini and Lowengrub 2004 Cristini, V.; Lowengrub, J. Three-dimensional crystal growth—II: nonlinear simulation and control of the Mullins–Sekerka instability. _Journal of Crystal Growth_ 2004, _266_ , 552–567
  * Wang et al. 1993 Wang, S.-L.; Sekerka, R.; Wheeler, A.; Murray, B.; Coriell, S.; Braun, R.; McFadden, G. Thermodynamically-consistent phase-field models for solidification. _Physica D: Nonlinear Phenomena_ 1993, _69_ , 189–200
  * Karma and Rappel 1996 Karma, A.; Rappel, W.-J. Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics. _Physical review E_ 1996, _53_ , R3017
  * Bouvet et al. 2022 Bouvet, L.; Calonne, N.; Flin, F.; Geindreau, C. Snow Equi-Temperature Metamorphism Described by a Phase-Field Model Applicable on Micro-Tomographic Images: Prediction of Microstructural and Transport Properties. _Journal of Advances in Modeling Earth Systems_ 2022, _14_ , e2022MS002998
  * Kaempfer and Plapp 2009 Kaempfer, T. U.; Plapp, M. Phase-field modeling of dry snow metamorphism. _Physical Review E_ 2009, _79_ , 031502
  * Reitzle et al. 2019 Reitzle, M.; Ruberto, S.; Stierle, R.; Gross, J.; Janzen, T.; Weigand, B. Direct numerical simulation of sublimating ice particles. _International Journal of Thermal Sciences_ 2019, _145_ , 105953
  * Fei et al. 2022 Fei, L.; Qin, F.; Wang, G.; Luo, K. H.; Derome, D.; Carmeliet, J. Droplet evaporation in finite-size systems: Theoretical analysis and mesoscopic modeling. _Physical Review E_ 2022, _105_ , 025101
  * Schweigler et al. 2017 Schweigler, K.; Said, M. B.; Seifritz, S.; Selzer, M.; Nestler, B. Experimental and numerical investigation of drop evaporation depending on the shape of the liquid/gas interface. _International Journal of Heat and Mass Transfer_ 2017, _105_ , 655–663
  * Ronsin et al. 2020 Ronsin, O. J.; Jang, D.; Egelhaaf, H.-J.; Brabec, C. J.; Harting, J. A phase-field model for the evaporation of thin film mixtures. _Physical Chemistry Chemical Physics_ 2020, _22_ , 6638–6652
  * Wang et al. 2021 Wang, Z.; Zheng, X.; Chryssostomidis, C.; Karniadakis, G. E. A phase-field method for boiling heat transfer. _Journal of Computational Physics_ 2021, _435_ , 110239
  * Fu et al. 2016 Fu, X.; Cueto-Felgueroso, L.; Juanes, R. Thermodynamic coarsening arrested by viscous fingering in partially miscible binary mixtures. _Physical Review E_ 2016, _94_ , 033111
  * Shimizu and Tanaka 2017 Shimizu, R.; Tanaka, H. Impact of complex topology of porous media on phase separation of binary mixtures. _Science Advances_ 2017, _3_ , eaap9570
  * Mao et al. 2019 Mao, S.; Kuldinow, D.; Haataja, M. P.; Košmrlj, A. Phase behavior and morphology of multicomponent liquid mixtures. _Soft Matter_ 2019, _15_ , 1297–1311
  * Bhopalam et al. 2022 Bhopalam, S. R.; Bueno, J.; Gomez, H. Elasto-capillary fluid–structure interaction with compound droplets. _Computer Methods in Applied Mechanics and Engineering_ 2022, _400_ , 115507
  * Yue 2020 Yue, P. Thermodynamically consistent phase-field modelling of contact angle hysteresis. _Journal of Fluid Mechanics_ 2020, _899_ , A15
  * Fu et al. 2017 Fu, X.; Cueto-Felgueroso, L.; Juanes, R. Viscous fingering with partially miscible fluids. _Physical Review Fluids_ 2017, _2_ , 104001
  * Cueto-Felgueroso et al. 2018 Cueto-Felgueroso, L.; Fu, X.; Juanes, R. Pore-scale modeling of phase change in porous media. _Physical Review Fluids_ 2018, _3_ , 084302
  * Fu et al. 2018 Fu, X.; Cueto-Felgueroso, L.; Juanes, R. Nonequilibrium Thermodynamics of Hydrate Growth on a Gas-Liquid Interface. _Physical Review Letters_ 2018, _120_ , 144501
  * Fu et al. 2020 Fu, X.; Jimenez-Martinez, J.; Nguyen, T. P.; Carey, J. W.; Viswanathan, H.; Cueto-Felgueroso, L.; Juanes, R. Crustal fingering facilitates free-gas methane migration through the hydrate stability zone. _Proceedings of the National Academy of Sciences_ 2020, _117_ , 31660–31664
  * Zhang et al. 2022 Zhang, J.; Yuan, J.; Yuan, H. Numerical and experimental investigation of phase change characteristics and desalination effect of seawater freezing desalination in narrow channels. _Desalination_ 2022, _544_ , 116145
  * Gomez et al. 2019 Gomez, H.; Bures, M.; Moure, A. A review on computational modelling of phase-transition problems. _Philosophical Transactions of the Royal Society A_ 2019, _377_ , 20180203
  * Boyer and Lapuerta 2006 Boyer, F.; Lapuerta, C. Study of a three component Cahn-Hilliard flow model. _ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique_ 2006, _40_ , 653–687
  * Steinbach 2009 Steinbach, I. Phase-field models in materials science. _Modelling and Simulation in Materials Science and Engineering_ 2009, _17_ , 073001
  * Folch and Plapp 2005 Folch, R.; Plapp, M. Quantitative phase-field modeling of two-phase growth. _Phys. Rev. E_ 2005, _72_ , 011602
  * Tóth et al. 2015 Tóth, G. I.; Pusztai, T.; Gránásy, L. Consistent multiphase-field theory for interface driven multidomain dynamics. _Phys. Rev. B_ 2015, _92_ , 184105
  * Allen and Cahn 1979 Allen, S. M.; Cahn, J. W. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. _Acta Metallurgica_ 1979, _27_ , 1085–1095
  * Flatau et al. 1992 Flatau, P. J.; Walko, R. L.; Cotton, W. R. Polynomial fits to saturation vapor pressure. _Journal of Applied Meteorology_ 1992, 1507–1513
  * Wexler 1977 Wexler, A. Vapor pressure formulation for ice. _Journal of Research of the National Bureau of Standards. Section A, Physics and Chemistry_ 1977, _81_ , 5
  * Massman 1998 Massman, W. A review of the molecular diffusivities of H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2 in air, O2 and N2 near STP. _Atmospheric Environment_ 1998, _32_ , 1111–1127
  * Miyoshi and Takaki 2016 Miyoshi, E.; Takaki, T. Validation of a novel higher-order multi-phase-field model for grain-growth simulations using anisotropic grain-boundary properties. _Computational Materials Science_ 2016, _112_ , 44–51
  * Moelans 2022 Moelans, N. New phase-field model for polycrystalline systems with anisotropic grain boundary properties. _Materials & Design_ 2022, _217_ , 110592
  * Karma and Rappel 1998 Karma, A.; Rappel, W.-J. Quantitative phase-field modeling of dendritic growth in two and three dimensions. _Physical review E_ 1998, _57_ , 4323
  * Karma 2001 Karma, A. Phase-Field formulation for quantitative modeling of alloy solidification. _Phys. Rev. Lett._ 2001, _87_ , 115701
  * Ohno et al. 2016 Ohno, M.; Takaki, T.; Shibuta, Y. Variational formulation and numerical accuracy of a quantitative phase-field model for binary alloy solidification with two-sided diffusion. _Phys. Rev. E_ 2016, _93_ , 012802
  * Chung and Hulbert 1993 Chung, J.; Hulbert, G. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method. 1993,
  * Jansen et al. 2000 Jansen, K. E.; Whiting, C. H.; Hulbert, G. M. A generalized-$\alpha$ method for integrating the filtered Navier–Stokes equations with a stabilized finite element method. _Computer Methods in Applied Mechanics and Engineering_ 2000, _190_ , 305–319
  * Karma and Rappel 1999 Karma, A.; Rappel, W.-J. Phase-field model of dendritic sidebranching with thermal noise. _Physical review E_ 1999, _60_ , 3614
  * Cates and Nardini 2023 Cates, M.; Nardini, C. Classical Nucleation Theory for Active Fluid Phase Separation. _Physical Review Letters_ 2023, _130_ , 098203
  * Huang et al. 2022 Huang, Z.; Lin, G.; Ardekani, A. M. A consistent and conservative Phase-Field model for thermo-gas-liquid-solid flows including liquid-solid phase change. _Journal of Computational Physics_ 2022, _449_ , 110795
  * Snoeijer and Brunet 2012 Snoeijer, J. H.; Brunet, P. Pointy ice-drops: How water freezes into a singular shape. _American Journal of Physics_ 2012, _80_ , 764–771
  * Zeng et al. 2022 Zeng, H.; Lyu, S.; Legendre, D.; Sun, C. Influence of gravity on the freezing dynamics of drops on a solid surface. _Physical Review Fluids_ 2022, _7_ , 103605
  * Seguy et al. 2023 Seguy, L.; Protière, S.; Huerre, A. Role of geometry and adhesion in droplet freezing dynamics. _Physical Review Fluids_ 2023, _8_ , 033601
  * Eames et al. 1997 Eames, I.; Marr, N.; Sabir, H. The evaporation coefficient of water: a review. _International Journal of Heat and Mass Transfer_ 1997, _40_ , 2963–2973
  * Libbrecht 2017 Libbrecht, K. G. Physical dynamics of ice crystal growth. _Annual Review of Materials Research_ 2017, _47_ , 271–295
  * Krishnamachari et al. 1996 Krishnamachari, B.; McLean, J.; Cooper, B.; Sethna, J. Gibbs-Thomson formula for small island sizes: Corrections for high vapor densities. _Physical Review B_ 1996, _54_ , 8899

##  1 Equivalence to models of two-phase transition

In this section we show that our wet snow metamorphism model is equivalent to
the two-phase models for solidification, sublimation, and evaporation when
only two phases are present in the system. Only for this section, we assume
that $\theta_{\mathrm{nucl}}=0$ (see Eqs. (1)–(3) in the main text). To prove
the equivalency, we first need to show that a two-phase system remains as a
two-phase system and the non-existing phase does not appear in $\Omega$. It is
trivial to show that the phase-fields equations in our model, Eqs. (1)–(3) in
the main text, have the property $\phi_{l}(\bm{x},t)=0\;\forall t$ if
$\phi_{l}(\bm{x},0)=0$ and $\theta_{\mathrm{nucl}}=0$, which proves this first
condition. The second condition constitutes the actual equivalence, which we
analyze for each phase transition in the following paragraphs.

###  1.1 Solidification

Solidification and melting occurs when $\phi_{a}=0$. We can rewrite the
$\phi_{i}$ evolution equation (Eq. (1) in the main text) imposing $\phi_{a}=0$
as

$\frac{\partial\phi_{i}}{\partial
t}=-3M_{0}\left(\frac{1}{\varepsilon}\phi_{i}(1-\phi_{i})(1-2\phi_{i})-\varepsilon\nabla^{2}\phi_{i}\right)-\alpha_{\mathrm{sol}}\phi_{i}^{2}\phi_{w}^{2}\frac{T-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}.$
(33)

The Lagrange multiplier in our formulation guarantees that
$\displaystyle\phi_{i}+\phi_{w}=0$, which implies that
$\displaystyle\frac{\partial\phi_{w}}{\partial
t}=-\frac{\partial\phi_{i}}{\partial t}$. Imposing $\phi_{a}=0$, the
temperature evolution equation (Eq. (12) in the main text) is expressed as:

$\rho(\overline{\phi})c_{p}(\overline{\phi})\frac{\partial T}{\partial
t}=\nabla\cdot\left[K(\overline{\phi})\nabla
T\right]+\rho(\overline{\phi})L_{\mathrm{sol}}\frac{\partial\phi_{i}}{\partial
t}.$ (34)

We can disregard the vapor density equation since there is no air phase.
Equations (33) and (34) constitute a two-phase model for solidification (and
melting) known as the generalized Stefan problem 32, 50. Under certain
parameter regimes, this model reproduces the kinetics defined by the Gibbs-
Thomson condition for the ice-liquid water interface, which reads

$\frac{T-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}=-d_{\mathrm{sol}}\chi-\beta_{\mathrm{sol}}v_{n},$
(35)

where $T$ is the interface temperature, $\beta_{\mathrm{sol}}$ is the kinetic
attachment coefficient, $d_{\mathrm{sol}}$ is the capillary length, $\chi$ is
the curvature of the interface (positive for spherical ice grains), and
$v_{n}$ is the normal velocity of the interface (positive for ice growth).
Karma and Rappel 32, 61 showed that the two-phase model captures the Gibbs-
Thomson condition kinetics when (i) $\varepsilon\rightarrow 0$ and (ii) the
model parameters $M_{0}$ and $\alpha_{\mathrm{sol}}$ follow the relations

$M_{0}=M_{\mathrm{sol}}=\frac{\varepsilon}{3\tau_{\mathrm{sol}}},\quad\alpha_{\mathrm{sol}}=\frac{\lambda_{\mathrm{sol}}}{\tau_{\mathrm{sol}}},\quad
d_{\mathrm{sol}}=a_{1}\frac{\varepsilon}{\lambda_{\mathrm{sol}}},\quad\beta_{\mathrm{sol}}=a_{1}\left(\frac{\tau_{\mathrm{sol}}}{\varepsilon\lambda_{\mathrm{sol}}}-a_{2}\frac{\varepsilon}{D^{*}_{iw}}\right),$
(36)

where $a_{1}\approx 5$, $a_{2}\approx 0.1581$ 111The $a_{1}$ and $a_{2}$
values are computed by following the same procedure as in Karma and Rappel 61,
applied to Eqs. (33) and (34)., and $D^{*}_{iw}\approx(D_{i}+D_{w})/2$
represents thermal diffusivity, such that $D_{j}=K_{j}/(\rho_{j}c_{p,j})$ for
$j=\\{i,w\\}$. However, we note that the above relations assume equal thermal
diffusivity around the freezing/melting interface, which is not true for ice
and water in our case. Asymptotic analysis assuming asymmetric diffusivities
have been conducted for compositional problems 62, 54, but the same type of
analysis for the thermal problems remain to be developed. Nevertheless,
Kaempfer and Plapp 34 has shown that the phase-field solidification model
recovers the sharp-interface solution if
$\displaystyle\varepsilon~{}<~{}\beta_{\text{sol}}{D^{*}_{iw}}\approx
10^{-5}$m, which is true in the simulations presented in this paper.

###  1.2 Sublimation

Sublimation and deposition occurs when $\phi_{w}=0$. After imposing
$\phi_{w}=0$, Eqs. (1), (12), and (14) in the main text are expressed as

$\displaystyle\frac{\partial\phi_{i}}{\partial
t}=-3M_{0}\left(\frac{1}{\varepsilon}\phi_{i}(1-\phi_{i})(1-2\phi_{i})-\varepsilon\nabla^{2}\phi_{i}\right)+\alpha_{\mathrm{sub}}\phi_{i}^{2}\phi_{a}^{2}\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{i}}$
(37) $\displaystyle\rho(\overline{\phi})c_{p}(\overline{\phi})\frac{\partial
T}{\partial t}=\nabla\cdot\left[K(\overline{\phi})\nabla
T\right]-\rho(\overline{\phi})L_{\mathrm{sub}}\frac{\partial\phi_{a}}{\partial
t},$ (38) $\displaystyle\frac{\partial(\phi_{a}\rho_{v})}{\partial
t}=\nabla\cdot\left[\phi_{a}D_{v}(T)\nabla\rho_{v}\right]+\rho_{\mathrm{SE}}\frac{\partial\phi_{a}}{\partial
t},$ (39)

with $\partial\phi_{a}/\partial t=-\partial\phi_{i}/\partial t$. These
equations constitute the two-phase sublimation model 34 and capture the
kinetics defined by the Gibbs-Thomson condition

$\frac{\rho_{v}-\rho_{vs}^{I}(T)}{\rho_{vs}^{I}(T)}=d_{\mathrm{sub}}\chi+\beta_{\mathrm{sub}}v_{n},$
(40)

as long as the parameters $M_{0}$, $\alpha_{\mathrm{sub}}$, and
$\rho_{\mathrm{SE}}$ follow the relations 34:

$\displaystyle
M_{0}=M_{\mathrm{sub}}=\frac{\varepsilon}{3\tau_{\mathrm{sub}}},\quad\alpha_{\mathrm{sub}}=\frac{\lambda_{\mathrm{sub}}}{\tau_{\mathrm{sub}}},\quad\rho_{\mathrm{SE}}=\rho_{i},$
$\displaystyle
d_{\mathrm{sub}}\frac{\rho_{vs}^{I}}{\rho_{i}}=a_{1}\frac{\varepsilon}{\lambda_{\mathrm{sub}}},\quad\beta_{\mathrm{sub}}\frac{\rho_{vs}^{I}}{\rho_{i}}=a_{1}\left(\frac{\tau_{\mathrm{sub}}}{\varepsilon\lambda_{\mathrm{sub}}}-a_{2}\frac{\varepsilon}{D^{*}_{ia}}-a_{2}\frac{\varepsilon}{D_{v0}}\right),$
(41)

where $\beta_{\mathrm{sub}}$ and $d_{\mathrm{sub}}$ are the Gibbs-Thomson
coefficients for sublimation and $D^{*}_{ia}$ represents the average thermal
diffusivity of ice and air. Again, we note that the above relations assume
equal thermal diffusivity around the interface, which is not true for ice and
air in our case. Asymptotic analysis assuming for this exact problem remains
to be developed. Meanwhile, Kaempfer and Plapp 34 has shown that the above
phase-field model recovers the sharp-interface solution if
$\displaystyle\varepsilon~{}<~{}\beta_{\text{sub}}\frac{\rho_{vs}^{I}}{\rho_{i}}{D^{*}_{ia}}$
and
$\displaystyle\varepsilon~{}<~{}\beta_{\text{sub}}\frac{\rho_{vs}^{I}}{\rho_{i}}{D_{v0}}$.
These two conditions yield the requirement that
$\displaystyle\varepsilon~{}<10^{-7}$m, which we do not satisfy in the
simulations presented. Thus, our numerical results do not strictly recover the
Gibbs-Thomson condition for air-ice interfaces.

###  1.3 Evaporation

The derivation for the evaporation case $(\phi_{i}=0)$ is analogous to the
sublimation case. The final two-phase model for evaporation (without flow) is
analogous to the 2-phase sublimation model (Eqs. (37)–(39)), where the
evaporation latent heat is expressed as
$(L_{\mathrm{eva}}=L_{\mathrm{sub}}-L_{\mathrm{sol}})$. Here, we simply list
the final relations between the model parameters $M_{0}$,
$\alpha_{\mathrm{eva}}$, $\rho_{\mathrm{SE}}$, and the Gibbs-Thomson
coefficients $\beta_{\mathrm{eva}}$ and $d_{\mathrm{eva}}$, which read

$\displaystyle
M_{0}=M_{\mathrm{eva}}=\frac{\varepsilon}{3\tau_{\mathrm{eva}}},\quad\alpha_{\mathrm{eva}}=\frac{\lambda_{\mathrm{eva}}}{\tau_{\mathrm{eva}}},\quad\rho_{\mathrm{SE}}=\rho_{w},$
$\displaystyle
d_{\mathrm{eva}}\frac{\rho_{vs}^{W}}{\rho_{w}}=a_{1}\frac{\varepsilon}{\lambda_{\mathrm{eva}}},\quad\beta_{\mathrm{eva}}\frac{\rho_{vs}^{W}}{\rho_{w}}=a_{1}\left(\frac{\tau_{\mathrm{eva}}}{\varepsilon\lambda_{\mathrm{eva}}}-a_{2}\frac{\varepsilon}{D^{*}_{wa}}-a_{2}\frac{\varepsilon}{D_{v0}}\right),$
(42)

where $D^{*}_{wa}$ represents the average thermal diffusivity of water and
air. We note again that the above relations assume equal thermal diffusivity
around the interface, which is not true for water and air in our case.
Asymptotic analysis assuming for this exact problem remains to be developed.
Meanwhile, Kaempfer and Plapp 34 has shown that the above phase-field model
recovers the sharp-interface solution if
$\displaystyle\varepsilon~{}<~{}\beta_{\text{eva}}\frac{\rho_{vs}^{W}}{\rho_{w}}{D^{*}_{wa}}$
and
$\displaystyle\varepsilon~{}<~{}\beta_{\text{eva}}\frac{\rho_{vs}^{I}}{\rho_{w}}{D_{v0}}$.
These two conditions yield the requirement that
$\displaystyle\varepsilon~{}<10^{-8}$m, which we do not satisfy in the
simulations presented. Thus, our numerical results do not strictly recover the
Gibbs-Thomson condition for air-water interfaces.

##  2 Additional details for numerical implementation

###  2.1 Determining temporal scaling coefficients

Wet snow metamorphism involves processes at different time scales. The
characteristic times for solidification $(t_{\mathrm{sol}})$, evaporation
$(t_{\mathrm{eva}})$, sublimation $(t_{\mathrm{sub}})$, thermal diffusion
$(t_{T,i},\,t_{T,w},\,\mathrm{and}\,t_{T,a})$, and vapor diffusion $(t_{v})$,
assuming a characteristic length scale of $L=10^{-4}\,$m, are

$\displaystyle t_{\mathrm{sol}}\equiv L/v_{n,\mathrm{sol}}\sim
1\,\mathrm{s},\;t_{\mathrm{sub}}\equiv L/v_{n,\mathrm{sub}}\sim
10^{4}\,\mathrm{s},\;t_{\mathrm{eva}}\equiv L/v_{n,\mathrm{eva}}\sim
10^{3}\,\mathrm{s},\;t_{T,i}\equiv L^{2}/D_{i}\sim 10^{-2}\,\mathrm{s},$
$\displaystyle t_{T,w}\equiv L^{2}/D_{w}\sim
10^{-1}\,\mathrm{s},\;t_{T,a}\equiv L^{2}/D_{a}\sim
10^{-3}\,\mathrm{s},\;t_{v}\equiv L^{2}/D_{v0}\sim 10^{-4}\,\mathrm{s},$ (43)

where thermal diffusivities are calculated as $D_{l}=K_{l}/(\rho_{l}c_{p,l})$,
for phase $l$, and we consider typical interface velocities $v_{n,j}$, for
$j=\\{\mathrm{sol},\mathrm{sub},\mathrm{eva}\\}$ (see Table 1). These values
show that sublimation and evaporation are several orders of magnitude slower
than solidification and any diffusion process. A monolithic numerical solver
requires the use of time steps that capture the lower characteristic time
(i.e., $t_{v}$), which would significantly increase the computational times.

Here, we leverage the procedure explained in Kaempfer and Plapp 34 to speed up
the simulations. This approach leverages the fact that $T$ and $\rho_{v}$
exhibit a quasi-steady state compared to the phase transition kinetics. The
procedure consists of multiplying the right-hand side of the $T$ and
$\rho_{v}$ equations (Eqs. (12) and (14) in the main text, respectively) by a
time-scale factor, so that the characteristic times for vapor and thermal
diffusion increase some orders of magnitude while $T$ and $\rho_{v}$ keep
displaying a quasi-steady state. This approach enables the use of larger time
steps without introducing noticeable errors in the simulations. We multiply
the right-hand side of the $\rho_{v}$ equation by the time-scale factor
$\xi_{v}<1$, which provokes an increase of the vapor diffusion characteristic
time to $t^{*}_{v}=t_{v}/\xi_{v}$. As long as $t^{*}_{v}$ does not exceed
$t_{\mathrm{sol}}$, $t_{\mathrm{sub}}$, and $t_{\mathrm{eva}}$, vapor
concentration displays a quasi-steady state and the overall problem kinetics
is not affected. This condition leads to $\xi_{v}>10^{-4}$. For the $T$
equation we consider the largest thermal diffusion characteristic time,
$t_{T,w}$. In Eq. (2.1) we observe that $t_{T,w}\sim t_{\mathrm{sol}}$, while
$t_{T,w}<t_{\mathrm{sub}}$ and $t_{T,w}<t_{\mathrm{eva}}$. Thus, $T$ displays
a quasi-steady state when solidification does not occur. We follow the same
procedure and multiply the right-hand side of the $T$ equation by the time-
scale factor $\xi_{T}<1$ when solidification does not take place. As long as
the new characteristic time $t^{*}_{T,w}=t_{T,w}/\xi_{T}$ does not exceed
$t_{\mathrm{sub}}$ and $t_{\mathrm{eva}}$, temperature displays a quasi-steady
state and the problem kinetics is not affected. This condition leads to
$\xi_{T}>10^{-4}$. In this work, we take $\xi_{v}=10^{-3}$ and $\xi_{T}=1$ or
$\xi_{T}=10^{-2}$ depending on whether or not solidification occurs in any
parts of the domain.

###  2.2 Kinetic parameters constraints on grid resolution

If we assume equal diffusivity for all phases, then the relations between the
Gibbs-Thomson parameters $(\beta_{j},d_{j})$ and our model parameters
$(M_{j},\alpha_{j})$ (for $j=\\{\mathrm{sol},\mathrm{sub},\mathrm{eva}\\}$)
can be derived from asymptotic analysis 32, 61 and are valid under the
following conditions:

$1.\;\chi\varepsilon<1,\quad\quad 2.\;\varepsilon v_{n,j}/D_{k}<1,\quad\quad
3.\;\lambda_{j}u_{j}<1,\quad\quad 4.\;\tau_{j}v_{n,j}/\varepsilon<1,$ (44)

for $k=\\{i,w,a,v0\\}$, where $\varepsilon$ represents the interface width,
$D_{k}$ is the thermal or vapor diffusivity coefficient, $v_{n,j}$ is the
normal velocity of the interface, $\lambda_{j}$ and $\tau_{j}$ are parameters
defined in Section 1 in the Supporting Information, and $u_{j}$ is the phase-
change factor that multiplies the terms $\alpha_{j}\phi_{l}^{2}\phi_{m}^{2}$
in the phase-field equations (for $\\{l,m\\}=\\{i,w,a\\}$; see Eqs. (23) and
(24) in the main text). Condition 1 imposes an interface width smaller than
the curvature radius, condition 2 guarantees that thermal and vapor diffusion
are faster than the interface motion, and conditions 3 and 4 ensure the
preservation of a tanh-profile (i.e., quasi-equilibrium profile) across the
interface. We use these conditions to guide the choice of $\varepsilon$ in our
model; however, we recognize that these conditions are only valid for equal
diffusivity problems, which is not true for the wet snow problem. These
conditions impose a restriction on $\varepsilon$, which represents the
numerical interface width and is an adjustable model parameter.

Table 1: Gibbs-Thomson parameters and $\varepsilon$ limitations. Index $j$ stands for solidification, sublimation and evaporation. The $\beta_{j}$, $d_{j}$, $M_{j}$, and $\alpha_{j}$ values considered in this work are listed in the bottom rows.  | Solidification | Sublimation | Evaporation
---|---|---|---
$\beta_{j}$ (s/m) | $\sim(10,100)$ | $(2,200)\times 10^{4}$ | $(2,200)\times 10^{3}$
$d_{j}$ (m) | $\sim 5\times 10^{-10}$ | $\sim 10^{-9}$ | $\sim 5\times 10^{-10}$
$u_{j}$ (-) | $\sim 10^{-3}$ | $\sim 10^{-8}$ | $\sim 10^{-8}$
$v_{n,j}$ (m/s) | $\sim(10^{-4},10^{-5})$ | $\sim(10^{-7},10^{-9})$ | $\sim(10^{-6},10^{-8})$
Cond. 1 (m) | $\varepsilon<\,\sim 10^{-4}$ | $\varepsilon<\,\sim 10^{-4}$ | $\varepsilon<\,\sim 10^{-4}$
Cond. 2 (m) | $\varepsilon<\,\sim 10^{-3}$ | $\varepsilon<\,\sim 10$ | $\varepsilon<\,\sim 10^{-1}$
Cond. 3 (m) | $\varepsilon<\,\sim 10^{-7}$ | $\varepsilon<\,\sim 10^{-7}$ | $\varepsilon<\,\sim 10^{-8}$
Cond. 4 (m) | $\varepsilon<\,\sim 10^{-7}$ | $\varepsilon<\,\sim 10^{-7}$ | $\varepsilon<\,\sim 10^{-7}$
Parameter values considered in this work (lead to $\varepsilon<10^{-6}$)
$\beta_{j}$ (s/m) | $125$ | $1.4\times 10^{5}$ | $1.53\times 10^{4}$
$d_{j}$ (m) | $4\times 10^{-9}$ | $10^{-7}$ | $6\times 10^{-8}$
$M_{j}$ (m/s) | $2.12\times 10^{-5}$ | $4.33\times 10^{-7}$ | $1.34\times 10^{-6}$
$\alpha_{j}$ (1/s) | $7.96\times 10^{4}$ | $1.3\times 10^{7}$ | $6.69\times 10^{7}$

To find a valid range for $\varepsilon$, we first estimate the Gibbs-Thomson
parameters $\beta_{j}$ and $d_{j}$ from the literature 72, 73, 34, 74 (see
Table 1). Note that $\beta_{j}$ may range several orders of magnitude
depending on ambient conditions and growth type 73. We next estimate the
typical values of $u_{j}$ and $v_{n,j}$ (Table 1), based on $T$ and $\rho_{v}$
distributions observed during the different phase transitions. We use the
Gibbs-Thomson conditions defined in Section 1 of the Supporting Information to
relate $u_{j}$ and $v_{n,j}$, under the assumption that the curvature
$\chi\approx 0$. In order to express $\lambda_{j}$ and $\tau_{j}$ as a
function of $\beta_{j}$ and $d_{j}$ (see Eqs. (36), (1.2), and (1.3)), we
assume $\rho^{I}_{vs}/\rho_{i}=\rho^{W}_{vs}/\rho_{w}=5\times 10^{-6}$, which
is the approximate value at the freezing point. For condition 1, we consider
typical ice curvatures $\chi\sim 10^{4}\,\mathrm{m}^{-1}$. With this
information, we solve for $\varepsilon$ in Eq. (44). We indicate the less
favorable valid range of $\varepsilon$ in Table 1. We find that condition 3 is
the most restrictive, such that $\varepsilon<10^{-8}\,\mathrm{m}$.

In order to speed up the simulations, we consider a value of
$d_{\mathrm{sol}}$, $d_{\mathrm{eva}}$ and $d_{\mathrm{sub}}$ one or two
orders of magnitude higher than its actual value (see bottom rows of Table 1).
This is because the $T$\- and $\rho_{v}$-terms in the Gibbs-Thomson conditions
(Eqs. (35) and (40)) dominate the interface motion compared to the curvature
terms ($d_{j}$) for $\chi\sim 10^{4}\,\mathrm{m}^{-1}$ in most cases. Doing
so, the $\varepsilon$ constraint in Eq. (44) relaxes to
$\varepsilon<10^{-6}\,\mathrm{m}$, which allows us to use a coarser spatial
discretization, whereas the wet snow metamorphism kinetics is not
significantly affected. If compared against experiments, deviations in our
modeled kinetics from the observed kinetics could be attributed to this
numerical treatment and, hence, the actual values of $d_{\mathrm{sol}}$,
$d_{\mathrm{eva}}$ and $d_{\mathrm{sub}}$ must be considered.

###  2.3 Parameter values

We list the value of the parameters used in our model at the bottom of Table 1
and in Table 3. For simplicity, we assume that the values of these parameters
are constant. The values of the parameters $k_{j}$ and $g_{j}$ in the
saturated vapor pressures are listed in Table 2. The parameter $\xi_{T}$ is
defined as

$\xi_{T}=\left\\{\begin{array}[]{ll}1&\mathrm{if}\quad\Gamma_{iw}>\epsilon,\\\
10^{-2}&\mathrm{if}\quad\Gamma_{iw}\leq\epsilon,\end{array}\right.$ (45)

where $\Gamma_{iw}$ is the length (in 2D) of the ice-water interface and the
parameter $\epsilon$ is close to zero. In particular, we compute the interface
as $\Gamma_{iw}=\int_{\Omega}\phi_{i}^{2}\phi_{w}^{2}\mathrm{d}x$ and we take
$\epsilon=2\times 10^{-4}\varepsilon$ in 1D and $\epsilon=2\times
10^{-6}\varepsilon\sqrt{S_{\Omega}}$ in 2D, where $S_{\Omega}$ is the area of
$\Omega$.

Table 2: Saturated vapor pressure of ice and water. Interpolation parameter values 56, 57. $j$ | 0 | 1 | 2 | 3
---|---|---|---|---
$k_{j}$ | $-0.5865\times 10^{4}$ | $0.2224\times 10^{2}$ | $0.1375\times 10^{-1}$ | $-0.3403\times 10^{-4}$
$g_{j}$ | $-0.2991\times 10^{4}$ | $-0.6017\times 10^{4}$ | $0.1888\times 10^{2}$ | $-0.2836\times 10^{-1}$
$j$ | 4 | 5 | 6 | 7
$k_{j}$ | $0.2697\times 10^{-7}$ | $0.6918$ | - | -
$g_{j}$ | $0.1784\times 10^{-4}$ | $-0.8415\times 10^{-9}$ | $0.4441\times 10^{-12}$ | $0.2859\times 10^{1}$
Table 3: Parameter values. The values of $\beta_{j}$, $d_{j}$, $M_{j}$, and $\alpha_{j}$ (where $j$ stands for sol, sub, and eva) are listed in the bottom of Table 1. We take $\varepsilon=2\times 10^{-7}\,$m for the simulations shown in Section 4.3.3 in the main text.  Parameter | Value | Units
---|---|---
$\varepsilon$ | $5\times 10^{-7}$ | m
$\Lambda$ | 1 | N/m
$\sigma_{iw}$ | 0.033 | N/m
$\sigma_{ia}$ | 0.109 | N/m
$\sigma_{wa}$ | 0.076 | N/m
$T_{\mathrm{melt}}$ | 0 |  ∘C
$P_{a}$ | 101325 | Pa
$\alpha_{\mathrm{nucl}}$ | $5\times 10^{6}$ | 1/m
$T_{\mathrm{nucl}}$ | variable |  ∘C
$\rho_{a}$ | 1.341 | Kg/m3
$\rho_{i}$ | 917 | Kg/m3
$\rho_{w}$ | 1000 | Kg/m3
$c_{p,a}$ | $1.044\times 10^{3}$ | J/(Kg$\,{}^{\circ}$C)
$c_{p,i}$ | $1.96\times 10^{3}$ | J/(Kg$\,{}^{\circ}$C)
$c_{p,w}$ | $4.2\times 10^{3}$ | J/(Kg$\,{}^{\circ}$C)
$K_{a}$ | 0.02 | W/(m$\,{}^{\circ}$C)
$K_{i}$ | 2.29 | W/(m$\,{}^{\circ}$C)
$K_{w}$ | 0.554 | W/(m$\,{}^{\circ}$C)
$L_{\mathrm{sol}}$ | $3.34\times 10^{5}$ | J/Kg
$L_{\mathrm{sub}}$ | $2.83\times 10^{6}$ | J/Kg
$D_{v0}$ | $2.178\times 10^{-5}$ | m2/s
$\xi_{v}$ | $10^{-3}$ | (-)
$\xi_{T}$ | see Eq. (45) | (-)

###  2.4 Initial and boundary conditions

For the phase fields and the temperature we consider the boundary conditions

$\nabla\phi_{i}\cdot\bm{n}=0,\quad\nabla\phi_{a}\cdot\bm{n}=0,\quad
T=T_{B}(\bm{x},t)\quad\mathrm{on}\;\partial\Omega,$ (46)

where $\partial\Omega$ is the boundary of $\Omega$ and the boundary
temperature $T_{B}$ may take different values on each boundary and may evolve
in time. For the vapor concentration, we consider a fixed concentration on the
boundary or we impose no-flux through the boundary, depending on each
simulation, which are expressed as

$\rho_{v}=\rho_{v,B}\quad\mathrm{or}\quad\nabla\rho_{v}\cdot\bm{n}=0\quad\mathrm{on}\;\partial\Omega,$
(47)

where $\rho_{v,B}$ is a constant value. The initial conditions of our problem
may be written as

$\phi_{l}(\bm{x},0)=0.5\left(1+\mathrm{tanh}\left(\frac{d_{l}(\bm{x})}{2\varepsilon}\right)\right),\quad
T(\bm{x},0)=T_{0}(\bm{x}),\quad\rho_{v}(\bm{x},0)=\rho^{I}_{vs}(T_{0}),$ (48)

where $l=\\{i,a\\}$, $d_{l}$ is the signed distance to the interface of phase
$\phi_{l}$ and $T_{0}(\bm{x})$ is a linear function consistent with the
temperature boundary condition $T_{B}$. The distance $d_{l}$ depends on the
geometry considered in each simulation.

###  2.5 Spatial, time discretization and variable regularization

We use the Finite Element Method to solve our model. In particular, we use a
uniform mesh composed of bilinear (linear in 1D) basis functions. We compute
the weak form by multiplying the model equations (Eqs. (23)–(24) in the main
text) by weighting functions, integrating over $\Omega$, and integrating by
parts taking into account the boundary conditions defined in Eqs. (46) and
(47). To obtain the Galerkin form, we substitute the problem unknowns and the
weighting functions by the basis functions. The phase field variables may take
values $\phi_{l}<0$ and $\phi_{l}>1$ during transient states, with
$l=\\{i,w,a\\}$. To avoid singularities and numerical issues we regularize the
functions $\rho$, $c_{p}$, $K$, and $\rho_{\mathrm{SE}}$: We restrict the
value of $\phi_{l}$ between 0 and 1 only when we compute the functions $\rho$,
$c_{p}$, and $K$. For the regularization of $\rho_{\mathrm{SE}}$ we (1)
implement an additional cubic interpolation locally around $\phi_{a}=0$ which
smooths out the discontinuity at that point, and (2) extend
$\rho_{\mathrm{SE}}$ off the ternary diagram constant in the normal direction
to the diagram sides.