# Freezing of sessile droplet and frost halo formation

Sivanandan Kavuri1, George Karapetsas2, Chander Shekhar Sharma3 and Kirti
Chandra<EMAIL_ADDRESS>1Department of Chemical Engineering,
Indian Institute of Technology Hyderabad, Sangareddy, 502 284, Telangana,
India
2Department of Chemical Engineering, Aristotle University of Thessaloniki,
Thessaloniki 54124, Greece
3Department of Mechanical Engineering, Indian Institute of Technology Ropar,
140001 Rupnagar, India

###### Abstract

The freezing of a sessile droplet unveils fascinating physics, characterised
by the emergence of a frost halo on the underlying substrate, the progression
of the liquid-ice interface, and the formation of a cusp-like morphology at
the tip of the droplet. We investigate the freezing of a volatile sessile
droplet, focusing on the frost halo formation, which has not been
theoretically explored. The formation of the frost halo is associated with the
inherent evaporation process in the early freezing stages. We observe a
negative evaporation flux enveloping the droplet in the initial stages, which
indicates that vapour produced during freezing condenses on the substrate
close to the contact line, forming a frost halo. The condensate accumulation
triggers re-evaporation, resulting in a temporal shift of the frost halo
region away from the contact line. Eventually, it disappears due to the
diffusive nature of the water vapour far away from the droplet. We found that
increasing the relative humidity increases the lifetime of the frost halo due
to a substantial reduction in evaporation that prolonged the presence of net
condensate on the substrate. Increasing liquid volatility increases the
evaporation flux, and condensation occurs closer to the droplet, as a higher
amount of vapour is in the periphery of the droplet. We also found that
decreasing the thermal conductivity of the substrate increases the total
freezing time. The slower freezing process is accompanied by increased
vaporized liquid, resulting in condensation with its concentration reaching
supersaturation.

Freezing, sessile droplet, frost halo, evaporation, lubrication approximation

##  I Introduction

Droplet freezing is ubiquitous in many practical applications, such as icing
in aircraft and wind turbines, spray, food and pharmaceutical engineering, and
natural phenomena, such as raindrop formation [1, 2, 3, 4]. Thus, several
researchers examined the intricate physics underlying the freezing behaviour
of a sessile droplet on supercooled substrates [5, 6]. However, despite
significant advancements, there are still considerable gaps in our fundamental
understanding of the freezing phenomenon and frost halo formation around a
sessile droplet during the freezing process because of the coupling between
the heat and mass transfer processes.

The freezing of a sessile droplet occurs in two distinct phases. The early
recalescence phase involves the rapid development of an ice-crystal scaffold
and the temperature rise in the remaining liquid [7, 6]. The energy released
as a liquid transforms into ice causes an increase in temperature.
Subsequently, the second, much slower phase follows, which causes the
remaining liquid to solidify through isothermal freezing. At the same time,
more heat is released as the freezing front (liquid-ice interface) propagates
towards the apex of the sessile droplet. The condensation of the vapor
produced by the evaporation induced by the release of latent heat during this
freezing process creates a frost halo around the frozen droplet [6].
Eventually, the apex of the frozen droplet exhibits an inflection point and a
cusp-like morphology [8, 9] as a result of the volume expansion [10] caused by
the density difference between the ice and water phases. The different stages,
namely the early recalescence, the intermediate stage when the freezing front
moves towards the apex of the droplet and the final stage associated with the
formation of the cusp at the tip of the droplet, during freezing of a water
droplet on a hydrophilic Polymethyl methacrylate (PMMA) substrate obtained
experimentally, are illustrated in figure 1.

The earlier experimental study by Hu and Jin [7] investigated the icing of a
water droplet by examining the temperature field and demonstrated that the
drop freezes from its bottom. Enríquez  _et al._ [11], Marin  _et al._ [12]
employed shadowgraphy to investigate the freezing of a water droplet on a
copper plate immersed in a cooling bath of ethylene-glycol, ethanol and dry
ice. They examined the temporal evolution of the freezing front and cusp
formation at the tip of the droplet. Peng  _et al._ [13], Pan  _et al._ [14]
studied the effect of wettability on the freezing behaviour of a water
droplet. Peng  _et al._ [13] considered a superhydrophobic aluminum substrate,
whereas Pan  _et al._ [14] used micropatterned stainless steel surfaces to
regulate wettability. A few researchers also examined the effect of substrate
curvature on the freezing behaviour of droplets [15, 16, 17, 18]. It is to be
noted that the aforementioned studies and several others [19, 6, 20, 21, 22,
23, 24, 25] focused on the onset of freezing, freezing time and cusp formation
at the tip of the droplet. They emphasised the influence of various factors
(such as the wettability and roughness of the substrate and environmental
conditions) on freezing dynamics. Jung  _et al._ [6] investigated the
formation of the frost halo and its subsequent freezing in the vicinity of the
droplet by conducting experiments in controlled humidity conditions. Jung  _et
al._ [6] considered three different substrates, namely PMMA, polished titanium
and copper substrates, which have heat conductivities ranging by three orders
of magnitude but exhibit similar wettability. By carefully balancing the heat
diffusion and vapour transfer, they observed frost halo formation on the
substrate surrounding the droplet. A few studies also reported condensation
frosting under saturation conditions around a freezing drop produced due to
evaporation from the drop [26, 27, 28]. Further, the consequence of the frost
halo formation during the freezing of multiple droplets was studied by Graeber
_et al._ [29]. They demonstrated that the vapour produced by the freezing
process could alter the saturation conditions around neighbouring drops. This
vapour can condense and freeze when the surrounding air is cooler than the
saturation temperature. As a result, it forms tiny ice nuclei that
significantly enhance the freezing dynamics. It has been recognised that drops
do not freeze in isolation, and inter-droplet interactions, triggered by
evaporation and halo formation, play a significant role in their frost
spreading across the substrate [26]. Moreover, the change in the nucleation
mechanism due to the frost formation alters the icephobicity of a substrate.

Figure 1: Different stages observed during the freezing of a water droplet on
a hydrophilic Polymethyl methacrylate (PMMA) substrate. Courtesy: Vineet Kumar
- some preliminary experiments conducted in our group.

A few researchers also studied the propagation of the freezing front, volume
expansion and cusp formation at the tip of the droplet by conducting numerical
simulations of the Navier-stokes equations in the framework of front-tracking
[30, 31], couple level-set and volume-of-fluid [32] and lattice Boltzmann
methods by incorporating the many-body dissipative particle dynamics with the
energy conservation [33, 34]. Tembely  _et al._ [35] investigated the impact
of a supercooled droplet on a cold substrate using a volume-of-fluid method.
Ajaev and Davis [36] conducted theoretical and experimental studies to examine
the effect of tri-junction conditions corresponding to the fixed contact
angle, fixed contact line, fixed growth angle and dynamic growth angle on the
evolution of three interfaces. They showed that the dynamic growth angle
condition could only predict the cusp formation observed in their experiments.
Zadravzil  _et al._ [9] and Tembely and Dolatabadi [37] investigated the
freezing dynamics of a two-dimensional sessile droplet by employing the
lubrication approximation. They treated the shear stress singularity at the
solid-liquid interface using a precursor layer technique and disjoining
pressure formulation [38]. Zadravzil  _et al._ [9] investigated the effect of
a cooler porous substrate on the freezing of a sessile droplet. They analysed
the spreading and imbibition that occurred along with the solidification
inside the pores, resulting in the occlusion and contact line arrest of the
droplet. Recently, Tembely and Dolatabadi [37] considered the curvature
expression of the droplet-air and liquid-ice interfaces, disjoining pressure,
the Gibbs-Thomson effect, natural convection, and substrate thermal and
physicochemical properties, and investigated the impact of these parameters on
the temporal evolution of the freezing front and cusp formation. It is to be
noted that none of these numerical studies has considered evaporation and the
recalescence phase during the freezing of sessile droplets. However, a few
previous lubrication-based studies have examined only evaporation without
considering freezing in sessile droplets on heated substrates (e.g., see Refs.
[39, 40, 41, 42]). Stan  _et al._ [43] observed that a frozen drop suspended
in a vacuum migrates as a consequence of symmetry breaking driven by a
localised increase in evaporation due to the heat release during the
recalescence phase of freezing. In the experimental study conducted by
Sebilleau  _et al._ [44], the freezing of a water drop on a copper substrate
was investigated under the conditions of a substrate temperature of
$-15^{\circ}$C and an ambient temperature ranging from $15^{\circ}$C to
$20^{\circ}$C. Their findings revealed a substantial impact of relative
humidity on the freezing front propagation rates, while wettability did not
significantly influence the process.

To summarize, the freezing behaviour, the propagation of the freezing front
and the cusp formation at the apex of the frozen sessile droplet have been
studied experimentally, numerically and theoretically. A sessile droplet
undergoes frost halo formation around the droplet at the early recalescence
phase. Although the frost halo was observed experimentally by Jung  _et al._
[6], it has not been theoretically investigated. To the best of our knowledge,
the present study is the first attempt to model the formation of the frost
halo using a lubrication approximation by incorporating evaporation and
condensation of the vapor produced due to the release of the latent heat
during freezing. We also estimate the temporal evolution of the width of the
frost halo after the recalescence phase for different dimensionless
parameters.

The rest of this paper is organized as follows. In §II, we present the problem
formulation, governing equations, scalings and the numerical procedure adopted
in the framework of lubrication approximation. After establishing the
mechanism of the frost halo formation, a detailed parametric study is
conducted in §III. Finally, concluding remarks are provided in §IV.

##  II Problem formulation

The freezing phenomenon of a sessile liquid droplet placed on a cold,
horizontal solid substrate is investigated numerically using the lubrication
approximation. Figure 2 depicts a schematic diagram of a two-dimensional (2D)
droplet at a particular instant, $t$, during freezing, along with various
parameters used in our modelling. The bottom of the substrate with thickness
$(d_{w})$, thermal conductivity $(\lambda_{w})$ and specific heat $(C_{pw})$
is maintained at a constant temperature, $T_{c}$. We assume that, at $t=0$,
the early recalescence phase has already occurred with the liquid reaching the
melting temperature. At this point, the slower solidification step driven by
the heat conduction ensues, and a very thin layer of ice of uniform thickness,
$S_{\infty}$, has formed along the solid substrate. The height and half-width
of the droplet are denoted by $H$ and $L$, respectively. The liquid is assumed
to be incompressible and Newtonian, with constant density $(\rho_{l})$,
specific heat capacity $(C_{pl})$, thermal conductivity $(\lambda_{l})$ and
viscosity $(\mu_{l})$. The surface tension of the liquid-gas interface
$(\gamma_{lg})$, density of frozen solid phase $(\rho_{s})$, specific heat
capacity $(C_{ps})$ and thermal conductivity $(\lambda_{s})$ are assumed to be
constant. A Cartesian coordinate system $(x,z)$ with its origin at the centre
of the droplet on the solid substrate is employed in our study as shown in
figure 2. Here, $z=s(x,t)$ and $z=h(x,t)$ represent the liquid-ice and liquid-
gas interfaces, respectively. In the present work, we consider the drop to be
very thin $(H\ll L)$. Thus, the aspect ratio of the droplet, $\epsilon=H/L$,
is assumed to be very small. This assumption permits us to use the lubrication
theory, employed below, to derive a set of evolution equations that govern the
freezing dynamics of the sessile droplet considering the liquid, ice and gas
phases. However, earlier studies demonstrated the validity of the lubrication
model for contact angles up to $60^{\circ}$ [45, 37].

###  II.1 Dimensional governing equations

####  II.1.1 Liquid phase

The dynamics in the liquid phase is governed by the mass, momentum and energy
conservation equations, which are given by

$\displaystyle\rho_{l}\left(\frac{\partial{\rm{\bf v}}}{\partial t}+{\rm{\bf
v}}\cdot\nabla{\rm{\bf v}}\right)$ $\displaystyle=$ $\displaystyle-\nabla
p+\mu_{l}\nabla^{2}{\rm{\bf v}},$ (1) $\displaystyle\nabla\cdot{\rm{\bf v}}$
$\displaystyle=$ $\displaystyle 0,$ (2)
$\displaystyle\rho_{l}C_{pl}\left(\frac{\partial T_{l}}{\partial t}+{\rm{\bf
v}}\cdot\nabla T_{l}\right)$ $\displaystyle=$
$\displaystyle\lambda_{l}\nabla^{2}T_{l},$ (3)

where ${\rm{\bf v}}$, $p$ and $T_{l}$ denote the velocity, pressure and
temperature in the liquid phase, respectively; $\nabla$ represents the
gradient operator. At the free surface ($z=h(x,t)$), the liquid velocity,
${\rm{\bf v}}$, and the velocity of the liquid-gas interface, ${\rm{\bf
v}}_{lg}$, are related as

${\rm{\bf v}}={\rm{\bf v}}_{lg}+(J_{v}/\rho_{l})\;{\rm\bf n}_{l},$ (4)

where $J_{v}$ denotes the evaporative flux and ${\rm\bf n}_{l}$ is the outward
unit normal on the interface. However, the tangential components of both
velocities, ${\rm{\bf v}}_{\tau}={\rm{\bf v}}-({\rm{\bf v}}\cdot{\rm\bf
n}_{l}){\rm\bf n}_{l}={\rm{\bf v}}_{lg}-({\rm{\bf v}}_{lg}\cdot{\rm\bf
n}_{l}){\rm\bf n}_{l}$, are the same. Moreover, at $z=h(x,t)$, the velocity
field satisfies the local mass, force and energy balance in the liquid and gas
phases [42]. Thus,

$\rho_{l}({\rm{\bf v}}-{\rm{\bf v}}_{lg})\cdot{\rm\bf n}_{l}=\rho_{g}({\rm{\bf
v}}_{g}-{\rm{\bf v}}_{lg})\cdot{\rm\bf n}_{l},$ (5)

where $\rho_{g}$ and ${\rm{\bf v}}_{g}$ denote the density and velocity field
of the gas phase, respectively. By ignoring the effect of vapour recoil and
considering viscous effects from the gas phase to be negligible, the
tangential and normal to the free surface components of the force balance at
the liquid-gas interface are given by

$\displaystyle-p+{\rm\bf n}_{l}\cdot\tau\cdot{\rm\bf n}_{l}$ $\displaystyle=$
$\displaystyle-p_{g}+\gamma_{lg}\kappa_{lg}+\Pi,$ (6) $\displaystyle{\rm\bf
n}_{l}\cdot\tau\cdot{\rm\bf t}_{l}$ $\displaystyle=$ $\displaystyle 0.$ (7)

Here, $\tau=\mu_{l}\left(\nabla{\rm{\bf v}}+\nabla{\rm{\bf v}}^{T}\right)$
denotes the extra stress tensor, $\kappa_{lg}=-\nabla_{s,l}\cdot{\rm\bf
n}_{l}$ denotes the mean curvature of the free surface, and
$\nabla_{s,l}=(I-{\rm\bf n}_{l}{\rm\bf n}_{l})\cdot\nabla$ represents the
surface gradient operator. The disjoining pressure ($\Pi$) that accounts for
the van der Waals interaction is defined as

$\Pi=A\left[\left(\frac{B}{h-s}\right)^{n}-\left(\frac{B}{h-s}\right)^{m}\right],$
(8)

where $A=A_{Ham}/B^{n}\geq 0$ is a constant that describes the magnitude of
the energy of the intermolecular interactions between the liquid-gas and
liquid-ice interfaces, and $B$ denotes the precursor layer thickness. Here,
$n>m>1$ and $A_{Ham}$ denotes the dimensional Hamaker constant.

(a)

(b)

Figure 2: (a) Schematic diagram of a sessile droplet undergoing freezing on a
solid substrate. Here, $S_{\infty}$ and $H_{\infty}$ are the thickness of an
initial microscopic ice layer and the thickness of the precursor layer,
respectively. (b) Schematic representation of the left ($m_{cl}$) and right
($m_{cr}$) edges of the condensate mass and the thickness of the halo
($d_{halo}=m_{cr}-m_{cl}$) in the vicinity of the two-dimensional (2D)
droplet.

The heat balance over the liquid-gas interface is given by

$J_{v}L_{v}+\lambda_{l}\nabla T_{l}\cdot{\rm\bf n}_{l}-\lambda_{g}\nabla
T_{g}\cdot{\rm\bf n}_{l}-h_{tc}(T_{g}-T_{lg})=0,$ (9)

where $\lambda_{g}$ and $T_{g}$ denote the thermal conductivity and
temperature of the gas phase, respectively. Here, $T_{lg}$ indicates the
temperature at the liquid-gas interface, $h_{tc}$ represents the convective
heat transfer coefficient, $L_{v}$ is the specific internal latent heat of
vaporization.

At the liquid-ice interface $(z=s(x,t))$, the velocity is given by

${\rm{\bf v}}={\rm{\bf v}}_{ls}-(J_{s}/\rho_{l})\;{\rm\bf n}_{s},$ (10)

where $J_{s}$ denotes the freezing mass flux and ${\rm{\bf v}}_{ls}$ the
velocity of the liquid-ice interface. As we impose the no-slip condition at
the liquid-ice interface, the tangential component of the velocity is given by

${\rm{\bf v}}\cdot{\rm\bf t}_{s}=0.$ (11)

Here, ${\rm\bf n}_{s}$ and ${\rm\bf t}_{s}$ are the outward unit normal and
unit tangential vectors on the liquid-ice interface, respectively.

####  II.1.2 Solid (ice) phase

The energy conservation equation in the solid (ice) phase is given by

$\rho_{s}C_{ps}\frac{\partial T_{s}}{\partial t}=\lambda_{s}\nabla^{2}T_{s},$
(12)

where $T_{s}$ denotes the temperature in the solid phase.

At the solid substrate ($z=0$), we impose continuity of temperature

$T_{s}=T_{w},$ (13)

where $T_{w}$ is the temperature of the substrate at $z=0$.

At the freezing front ($z=s(x,t)$), the boundary condition for the temperature
is expressed as

$T_{s}=T_{l}=T_{f}.$ (14)

We also assume that equilibrium temperature at the freezing front, $T_{f}$, is
the same as the melting temperature, $T_{m}$.

At $z=s(x,t)$, the conservation of mass and energy between the liquid and
solid phases leads to

$J_{s}=\rho_{l}({\rm{\bf v}}_{ls}-{\rm{\bf v}})\cdot{\rm\bf
n}_{s}=\rho_{s}({\rm{\bf v}}_{ls}-{\rm{\bf v}}_{s})\cdot{\rm\bf n}_{s},$ (15)
$\rho_{s}({\rm{\bf v}}_{ls}-{\rm{\bf v}}_{s})\cdot{\rm\bf
n}_{s}H_{s}-\lambda_{s}\nabla T_{s}\cdot{\rm\bf n}_{s}=\rho_{l}({\rm{\bf
v}}_{ls}-{\rm{\bf v}})\cdot{\rm\bf n}_{s}H_{l}-\lambda_{l}\nabla
T_{l}\cdot{\rm\bf n}_{s},$ (16)

where ${\rm{\bf v}}_{s}$ denotes the velocity of the ice phase; $H_{s}$ and
$H_{l}$ denote the enthalpy of the ice and liquid, respectively. By combining
eqs. (15) and (16) and assuming that ${\rm{\bf v}}_{s}=0$, we get

$J_{s}\Delta H_{sl}-\lambda_{s}\nabla T_{s}\cdot{\rm\bf
n}_{s}+\lambda_{l}\nabla T_{l}\cdot{\rm\bf n}_{s}=0,$ (17)

where $\Delta H_{sl}=H_{s}-H_{l}$ denotes the enthalpy jump at the liquid-ice
interface. Considering that $H_{s}=C_{ps}(T_{f}-T_{m})+L_{f}(T_{m})$ and
$H_{l}=C_{pl}(T_{f}-T_{m})$, we evaluate

$\Delta H_{sl}=(C_{ps}-C_{pl})(T_{f}-T_{m})+L_{f},$ (18)

where $L_{f}$ denotes the latent heat of fusion considering the melting
temperature, $T_{m}$ as the reference temperature. As will be shown below, eq.
(17) can be used to evaluate the position of the freezing front, $s(x,t)$.

Table 1: Typical values of the physical parameters considered in our simulations. These properties are for water-air system as considered by Jung  _et al._ [6]. Property | Notation | Value
---|---|---
Density of the liquid phase | $\rho_{l}$ | 1000 Kg m-3
Density of the frozen solid phase | $\rho_{s}$ | 900 Kg m-3
Density of the gas phase | $\rho_{g}$ |  $1.65\times 10^{-3}$ Kg m-3
Viscosity of the liquid phase | $\mu_{l}$ |  $10^{-3}$ Pa$\cdot$s
Viscosity of the gas phase | $\mu_{g}$ |  $1.81\times 10^{-5}$ Pa$\cdot$s
Melting temperature | $T_{m}$ | 273 K
Ambient temperature | $T_{g}$ | 258.5 K
Temperature at the bottom of the substrate | $T_{c}$ | 258.5 K
Thickness of the substrate | $d_{w}$ |  $3\times 10^{-3}$ m
Thermal conductivity of the substrate | $\lambda_{w}$ | 0.1 - 500 W m-1 K-1
Thermal conductivity of the liquid phase | $\lambda_{l}$ | 0.57 W m-1 K-1
Thermal conductivity of the frozen solid phase | $\lambda_{s}$ | 2.21 W m-1 K-1
Thermal conductivity of the gas phase | $\lambda_{g}$ | 0.02 W m-1 K-1
Specific heat capacity of the substrate | $C_{pw}$ | 300 - 3000 J Kg-1K-1
Specific heat capacity of the liquid phase | $C_{pl}$ | 4220 J Kg-1K-1
Specific heat capacity of the frozen solid phase | $C_{ps}$ | 2050 J Kg-1K-1
Specific heat capacity of the gas phase | $C_{pg}$ | 4220 J Kg-1K-1
Surface tension of the liquid-gas interface at $T_{m}$ | $\gamma_{lgo}$ | 0.07 N m-1
Surface tension of the liquid-ice interface | $\gamma_{ls}$ | 0.02 N m-1
Latent heat of evaporation | $L_{v}$ |  $2.45\times 10^{6}$ J Kg-1
Latent heat of fusion | $L_{f}$ |  $3.35\times 10^{5}$ J Kg-1
Diffusion coefficient of the vapour in the gas phase | $D_{m}$ |  $1.89\times 10^{-5}$ m2/s
Universal gas constant | $R_{g}$ | 8.314 J K-1 mol-1
Accommodation coefficient | $\alpha$ | $\approx 1$
Relative humidity | ${RH}$ | $0-1$
Heat transfer coefficient | $h_{tc}$ | 10 - 7000 W m-2 K-1

####  II.1.3 Gas phase

In order to account for situations when a droplet freezes in an environment
with varying humidity, we consider that the gas phase consists of both air and
vapour. Based on the assumption that $\lambda_{g}\ll\lambda_{l}$ and thus we
neglect thermal effects in the gas phase and thus consider that the gas has a
constant bulk temperature, $T_{g}$. The gas phase velocity (${\rm{\bf
v}}_{g}$) is assumed to be varied linearly between the liquid-gas interface
and far away from the droplet, such that ${\rm{\bf v}}_{g}\cdot{\rm\bf
t}_{l}={\rm{\bf v}}_{lg}\cdot{\rm\bf t}_{l}$ at $z=h$ and ${\rm{\bf v}}_{g}=0$
at $z=d_{g}$; the latter assumption is valid in the case where the gas phase
can be considered to be quiescent far from the droplet. The concentration
balance in the gas phase is given by

$\frac{\partial\rho_{v}}{\partial t}+{\rm{\bf
v}}_{g}\cdot\nabla\rho_{v}=D_{m}\nabla^{2}\rho_{v},$ (19)

where $\rho_{v}$ is the vapor concentration in the gas phase and $D_{m}$
represents the diffusion coefficient of the vapour in the gas phase. We
consider the case of a well-mixed gas phase where proper ventilation maintains
constant vapour concentration at some distance from the droplet $(\rho_{vi})$.
Thus, at ${z=d_{g}}$, $\rho_{v}=\rho_{vi}$.

The liquid in the droplet and the vapour in the gas phase are coupled as a
consequence of the evaporation and condensation at the liquid-gas interface.
In the gas phase, the vapour mass flux is related to the departure from the
uniform vapour density, which is given by

$J_{v}=-D_{m}\left({\rm\bf
n}_{l}\cdot\nabla\rho_{v}\right)\;\mbox{at}\;z=h(x,t).$ (20)

Additionally, the kinetic theory leads to a linear constitutive relation
between the mass and the departure from equilibrium at the interface, which is
known as the Hertz–Knudsen relationship [46, 47, 42]. This is given by

$J_{v}=\alpha\left(\frac{R_{g}T_{m}}{2\pi
M}\right)^{1/2}\left(\rho_{ve}(T_{lg})-\rho_{v}\right),$ (21)

where $R_{g}$ is the universal gas constant, $M$ is the molecular weight and
$\alpha$ is the accommodation coefficient (close to unity). The equilibrium
vapour concentration, $\rho_{ve}(T_{lg})$, can be obtained by employing a
linear temperature-dependent equation of state:

$\rho_{ve}(T_{lg})=\rho_{ve}(T_{m})\left[1+\frac{M(p-p_{g})}{\rho_{l}R_{g}T_{m}}+\frac{ML_{v}}{R_{g}T_{m}}\left(\frac{T_{lg}}{T_{m}}-1\right)\right].$
(22)

Here, $L_{v}$ is the latent heat of evaporation. The combination of eqs. (20)
and (21) leads to

$D_{m}\left({\rm\bf
n}_{l}\cdot\nabla\rho_{v}\right)=-\alpha\left(\frac{R_{g}T_{m}}{2\pi
M}\right)^{1/2}\left(\rho_{ve}(T_{lg})-\rho_{v}\right),$ (23)

which can be used to evaluate the local interfacial vapour concentration,
$\rho_{v}$. We use eq. (23) as the boundary condition at $z=h(x,t)$ to solve
eq. (19).

####  II.1.4 Solid substrate

The energy conservation equation in the solid substrate is given by

$\rho_{w}C_{pw}\frac{\partial T_{w}}{\partial t}=\lambda_{w}\nabla^{2}T_{w},$
(24)

where $\rho_{w}$ represents the density of the substrate. The above equation
is subjected to the continuity of thermal flux at the ice-substrate interface
($z=0$):

$\lambda_{s}\frac{\partial T_{s}}{\partial z}=\lambda_{w}\frac{\partial
T_{w}}{\partial z},$ (25)

and at the bottom of the substrate, ${z=-d_{w}}$,

$T_{w}=T_{c}.$ (26)

The properties of fluids and range of physical conditions considered in our
study are listed in Table 1.

###  II.2 Scaling and evolution equations

The governing equations and boundary conditions are nondimensionalized using
the following scalings (wherein tilde denotes the dimensionless variable):

$\begin{gathered}(x,z,h,d_{w},d_{g})=L(\tilde{x},\epsilon\tilde{z},\epsilon\tilde{h},\epsilon\tilde{D_{w}},\epsilon\tilde{D_{g}}),~{}t=\frac{L}{U}\tilde{t},~{}(u,w)=U(\tilde{u},\epsilon\tilde{w}),\\\
(u_{g},w_{g})=U(\tilde{u_{g}},\epsilon\tilde{w_{g}}),~{}(p,\Pi)=\frac{\mu_{l}UL}{H^{2}}(\tilde{p},\tilde{\Pi}),~{}T_{i}=\Delta
T\;\tilde{T}_{i}+T_{c}\;(i=l,s,w),\\\
\gamma_{i}=\gamma_{lgo}\tilde{\gamma}_{i}\;(i=lg,ls),~{}J_{s}=\epsilon\rho_{s}U\;\tilde{J}_{s},~{}J_{v}=\epsilon\rho_{ve}(T_{m})U\tilde{J}_{v},~{}\rho_{v}=\rho_{ve}(T_{m})\;\tilde{\rho}_{v},\\\
\nabla=\frac{1}{L}\tilde{\nabla},~{}\nabla_{s,i}=\frac{1}{L}\tilde{\nabla}_{s,i}\;(i=l,s),\end{gathered}$
(27)

where $\Delta T=T_{m}-T_{c}$, $\widetilde{\nabla}={\rm{\bf
e}}_{x}\widetilde{\partial}_{x}+{\rm{\bf
e}}_{z}\epsilon^{-1}\widetilde{\partial}_{z}$,
$\widetilde{\nabla}_{s,i}=(I-{\rm\bf n}{\rm\bf
n})\cdot\widetilde{\nabla}\;(i=l,s)$ and $U=\epsilon^{3}\gamma_{lgo}/\mu_{l}$
is the velocity scale. Henceforth, the tilde notation is suppressed, and
${\partial/\partial x}$, ${\partial/\partial z}$ and ${\partial/\partial t}$
are represented by the subscripts $x$, $z$ and $t$, respectively. By employing
these scalings, the governing equations are rendered dimensionless, and the
dimensionless groups that emerge are presented in Table 2.

Table 2: Dimensionless groups used in the present study. Dimensionless group | Description
---|---
$\epsilon=H/L$ | Droplet aspect ratio
$D_{w}=d_{w}/H$ | Scaled wall thickness
$D_{g}=d_{g}/H$ | Scaled height of the gas phase domain
$T_{v}=(T_{g}-T_{c})/(T_{m}-T_{c})$ | Scaled gas temperature
$D_{v}=\rho_{ve}(T_{m})/\rho_{l}$ | Density ratio (vapour-liquid)
$D_{s}=\rho_{s}/\rho_{l}$ | Density ratio (solid-liquid)
$\Lambda_{s}=\lambda_{s}/\lambda_{l}$ |  Thermal conductivity ratio (ice-liquid)
$\Lambda_{w}=\lambda_{w}/\lambda_{l}$ |  Thermal conductivity ratio (wall-liquid)
$\chi=\epsilon\rho_{ve}(T_{m})UL_{v}H/(\lambda_{l}\Delta T)$ |  Scaled latent heat of vaporization
$Bi=h_{tc}H/\lambda_{l}$ | Biot number
$A_{n}=HA/\gamma_{lgo}$ | Scaled Hamaker constant
$Ste=\lambda_{l}\Delta T/(\epsilon\rho_{s}UL_{f}H)$ | Stefan number
$Pe_{v}={\epsilon UH/D_{m}}$ | Péclet number
$V_{\rho,r}=\rho_{vi}/\rho_{ve}(T_{m})$ | Vapour density ratio
$K=\epsilon U/\alpha\sqrt{2\pi M/\left(R_{g}T_{m}\right)}$ | Dimensionless number associated with equilibrium condition at the interface
$\Delta=\mu_{l}ULM/\left(H^{2}\rho_{l}R_{g}T_{m}\right)$ | Dimensionless number for pressure effect on evaporation
$\Psi=L_{v}M\Delta T/\left(R_{g}T_{m}^{2}\right)$ | Dimensionless number associated with evaporative cooling

Incorporating the lubrication approximation ($\epsilon\ll 1$) and following
the typical procedure (given in detail in the Appendix §A), we obtain the
following evolution equations for the height of the liquid-gas and the liquid-
ice interface:

$\partial_{t}h-\partial_{t}s=-\partial_{x}q_{l}-D_{v}J_{v}-D_{s}J_{s},$ (28)
$\partial_{t}s=Ste\left(\frac{\Lambda_{w}T_{f}}{D_{w}+s\Lambda_{w}/\Lambda_{s}}+\chi
J_{v}-Bi(T_{v}-T_{l}\big{|}_{h})\right),$ (29)

where

$q_{l}=\frac{\partial_{x}p}{2}\left(\frac{h^{3}}{3}-s^{2}h+\frac{2s^{3}}{3}\right)-h\partial_{x}p\left(\frac{h^{2}}{2}-sh+\frac{s^{2}}{2}\right).$
(30)

Moreover, the temperature profiles in the liquid, ice and solid substrate are
given by

$\displaystyle T_{l}$ $\displaystyle=$ $\displaystyle\left(-\chi
J_{v}+Bi(T_{v}-T_{l}\big{|}_{h})\right)(z-s)+T_{f},$ (31) $\displaystyle
T_{s}$ $\displaystyle=$
$\displaystyle\frac{T_{f}}{D_{w}+s\Lambda_{w}/\Lambda_{s}}\left(D_{w}+z\frac{\Lambda_{w}}{\Lambda_{s}}\right),$
(32) $\displaystyle T_{w}$ $\displaystyle=$
$\displaystyle\frac{T_{f}}{D_{w}+s\;\Lambda_{w}/\Lambda_{s}}(z+D_{w}).$ (33)

####  II.2.1 Gas phase - Kármán-Pohlhausen approximation

Special care is needed in order to evaluate the vapour concentration profile
in the gas phase in order to retain the advection terms. To this end, we
employ the Kármán-Pohlhausen integral approximation and define the integrated
form of $\rho_{v}$, which is given by

$\int_{h}^{D_{g}}\rho_{v}dz=f.$ (34)

In order to be able to evaluate this integral, we need to prescribe the form
of $\rho_{v}$ as a function of the vertical coordinate. To this end, we assume
that $\rho_{v}$ can be approximated by a polynomial of the form

$\rho_{v}=c_{1}z^{2}+c_{2}z+c_{3}.$ (35)

By substituting the corresponding polynomial in eq. (34) and applying the
appropriate boundary conditions at the liquid-gas interface and far from the
droplet (see eqs. (A.18) and (A.20) in the Appendix §A), it is possible to
evaluate the polynomial constants and eventually derive the following
expressions for the constants $c_{1}$, $c_{2}$ and $c_{3}$.

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle\frac{f-\frac{Pe_{v}J_{v}}{2}(D_{g}-h)^{2}-V_{\rho,r}(D_{g}-h)}{\frac{2}{3}(h-D_{g})^{3}},$
(36) $\displaystyle c_{2}$ $\displaystyle=$ $\displaystyle-
Pe_{v}J_{v}-2h\left[\frac{f-\frac{Pe_{v}J_{v}}{2}(D_{g}-h)^{2}-V_{\rho,r}(D_{g}-h)}{\frac{2}{3}(h-D_{g})^{3}}\right],$
(37) $\displaystyle c_{3}$ $\displaystyle=$ $\displaystyle
V_{\rho,r}-{D_{g}}^{2}\left[\frac{f-\frac{Pe_{v}J_{v}}{2}(D_{g}-h)^{2}-V_{\rho,r}(D_{g}-h)}{\frac{2}{3}(h-D_{g})^{3}}\right]+Pe_{v}J_{v}D_{g}$
(38)
$\displaystyle+2hD_{g}\left[\frac{f-\frac{Pe_{v}J_{v}}{2}(D_{g}-h)^{2}-V_{\rho,r}(D_{g}-h)}{\frac{2}{3}(h-D_{g})^{3}}\right],$

where $V_{\rho,r}=\rho_{vi}/\rho_{ve}(T_{m})$ denotes that vapour density
ratio, which is related with the relative humidity ($RH$) as

$RH=\frac{V_{\rho,r}}{1+\Psi(T_{v}-1)}.$ (39)

By integrating the vapour concentration balance (eq. A.17) and using the
boundary conditions, we get the following integrated form of the concentration
equation which can eventually be solved with respect to $f$:

$\displaystyle Pe_{v}\left[\frac{\partial f}{\partial t}+\frac{\partial
g_{v}}{\partial x}-\rho_{v}\big{|}_{h}(D_{v}J_{v})\right]$ $\displaystyle=$
$\displaystyle\frac{\partial\rho_{v}}{\partial
z}\big{|}_{D_{g}}-\frac{\partial\rho_{v}}{\partial z}\big{|}_{h}$
$\displaystyle+$ $\displaystyle\epsilon^{2}\left[\frac{\partial}{\partial
x}\left(\frac{\partial f}{\partial
x}+\rho_{v}\big{|}_{h}h_{x}\right)+\frac{\partial\rho_{v}}{\partial
x}\big{|}_{h}h_{x}\right].$ (40)

In the above equation $g_{v}$ is defined as

$\int_{h}^{D_{g}}u_{g}\rho_{v}dz=g_{v}.$ (41)

To evaluate eq. (41), we need to prescribe the $x$-component of the velocity
profile in the gas phase, which can be approximated by considering the
following linear profile

$u_{g}=az+b.$ (42)

The constants $a$ and $b$ can be evaluated by simply considering that, at the
liquid-gas interface ($z=h(x,t)$), the velocity of the gas is equal to the
velocity of liquid

$u_{g}=u,~{}{\rm and}~{}w_{g}=w,$ (43)

and at the far-field ($z=D_{g}$) the gas phase is quiescent,

$u_{g}=0,~{}{\rm and}~{}w_{g}=0.$ (44)

Thus,

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{\frac{\partial_{x}p}{2}(h^{2}-s^{2})-h\partial_{x}p(h-s)}{(h-{D_{g}})},$
(45) $\displaystyle b$ $\displaystyle=$
$\displaystyle-{D_{g}}\left[\frac{\frac{\partial_{x}p}{2}(h^{2}-s^{2})-h\partial_{x}p(h-s)}{(h-{D_{g}})}\right].$
(46)

###  II.3 Initial and boundary conditions - Numerical procedure and validation

In our modelling, the droplet is deposited on a thin precursor layer that
resides on top of a thin layer of ice. The precursor layer in our study can be
related to the quasi-liquid layer, typically characterized by a thickness
ranging from 1-100 nm, as discussed by Nagata  _et al._ [48], Slater and
Michaelides [49]. While the quasi-liquid layer is a recognized phenomenon in
supersaturated conditions, its existence in undersaturated conditions remains
unestablished. Nevertheless, we note that a similar approach has successfully
been used in earlier theoretical studies of droplet freezing (see [9, 37]) to
alleviate the contact line singularity. The following initial conditions are
imposed on the domain.

$\displaystyle
h(x,t=0)=\max\left(h_{\infty}+s_{\infty}+1-x^{2},h_{\infty}+s_{\infty}\right),$
(47) $\displaystyle f(x,t=0)=V_{\rho,r}\left({D_{g}}-h(x,t=0)\right).$ (48)

The dimensionless equilibrium precursor layer thickness
($h_{\infty}=(h-s_{\infty})=H_{\infty}/{H}$) far from the droplet can be
estimated by considering that in this region, the fluid is flat with zero mean
curvature and sufficiently thin such that the attractive van der Waals forces
suppress evaporation. Therefore, by taking the constitutive equation for the
evaporation flux and setting it to zero, the dimensionless equilibrium
precursor thickness $h_{\infty}$ can be evaluated by solving the following
nonlinear equation

$1-{V_{\rho,r}}+\Delta\left(-\epsilon^{-2}A_{n}\left[\left(\frac{\beta}{h_{\infty}}\right)^{n}-\left(\frac{\beta}{h_{\infty}}\right)^{m}\right]\right)+\Psi\left(\frac{1+Bi({T_{v}}h_{\infty})}{1+Bi(h_{\infty})}-1\right)=0.$
(49)

In all our simulations, we have taken $\beta=0.01$. The initial thickness of
the thin ice layer is taken to be $s_{\infty}=10^{-3}$, but we have verified
that our findings remain unchanged when considering, e.g. one or two orders of
magnitude smaller values of $s_{\infty}$. To prevent the precursor layer from
freezing, we adopt a similar approach to that of Zadravzil  _et al._ [9] and
introduce the thickness-dependent Stefan number $(Ste(x))$, which is given by
the following expression.

$Ste(x)=\frac{1}{2}(1+\tanh[4\times 10^{3}((h-s)-1.4\beta)])Ste.$ (50)

Finally, we impose the following set of boundary conditions

$\displaystyle h_{x}(0,t)=h_{xxx}(0,t)$ $\displaystyle=$ $\displaystyle
h_{x}(x_{\infty},t)=h_{xxx}(x_{\infty},t)=0,$ (51) $\displaystyle
s_{x}(0,t)=s_{xxx}(0,t)$ $\displaystyle=$ $\displaystyle
s_{x}(x_{\infty},t)=s_{xxx}(x_{\infty},t)=0,$ (52) $\displaystyle
h(x_{\infty},t)-s(x_{\infty},t)$ $\displaystyle=$ $\displaystyle
h_{\infty},~{}s(x_{\infty},t)=s_{\infty},$ (53) $\displaystyle
f_{x}(0,t)=0,~{}f(x_{\infty},t)$ $\displaystyle=$
$\displaystyle{V_{\rho,r}}\left({D_{g}}-s(x_{\infty},t)-h_{\infty}\right).$
(54)

where $x_{\infty}$ represents the end of the domain.

The Galerkin Finite Element Method is utilized to discretize all dimensionless
governing equations. To mitigate boundary effects, our simulations incorporate
a computational domain that is sufficiently large, extending six times the
dimensionless initial radius. The optimum size of this computational domain is
determined by investigating the effect of its size. A grid convergence test
has been conducted for a representative parameter set considered in our study.
We found that the results obtained with 2001, 4801 and 8001 grid points are
indistinguishable (figure 15 in Appendix §B). Consequently, we selected 4801
grid points for the remaining simulations. The solutions are obtained using
the Newton-Raphson scheme, which is applied after the simulation has
progressed in time using an implicit Euler method and an adaptive time step.

####  II.3.1 Validation

In order to validate our solver, we simulate a case that was considered by
Zadravzil  _et al._ [9] in figure 3. This result mimics figure 15 of Zadravzil
_et al._ [9]. It shows the evolution of the freezing front, $s$ (dot-dashed
line) and the shape of the droplet, $h$ (solid line), placed on a cold
substrate for $Ste=0.04$, $T_{v}=0.5$, $A_{n}=6.25$,
$D_{s}=\Lambda_{S}=\Lambda_{W}=V_{\rho,r}=K=Bi=1$ and
$D_{w}=D_{v}=\Delta=\Psi=\chi=Pe_{v}=0$. Under this limit, we recover the
formulation of Zadravzil  _et al._ [9] for a non-porous substrate without
gravity effect. As depicted in figure 3, the ice-front and the morphology of
the liquid-gas interface in our study exhibit qualitative similarities to
those observed in Zadravzil  _et al._ [9].

Figure 3: Evolution of the freezing front, $s$ (dot-dashed lines) and shape of
the droplet, $h$ (solid line) placed on a cold substrate. The rest of the
dimensionless parameters are $\epsilon=0.2$, $Ste=0.04$, $T_{v}=0.5$,
$A_{n}=6.25$, $D_{s}=\Lambda_{S}=\Lambda_{W}=V_{\rho,r}=K=Bi=1$, and
$D_{w}=D_{v}=\Delta=\Psi=\chi=Pe_{v}=0$. This set of parameters is similar to
that of Zadravzil  _et al._ [9]. The value of the dimensionless total freezing
time $(t_{f})$ of the droplet is 15. Here, the evaporation is neglected.

An important feature of the freezing process, typically observed in
experiments, is the formation of a cusp at the drop apex when freezing ends.
Marin  _et al._ [12] experimentally observed a tip angle of $\sim
139^{\circ}$. They also reported that the tip angle remains independent of
substrate temperature, wettability, and solidification rate. In figure 16 of
Appendix §C, we consider a case for a typical set of parameters considered in
the present study ($Ste=1.22\times 10^{-3}$, $T_{v}=1$, $A_{n}=6.25$,
$D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.89$, $\Lambda_{W}=698$, $\chi=0.01$,
$K=8\times 10^{-4}$, $Bi=0.16$, $D_{w}=15$, $V_{\rho,r}=0.70$, $\epsilon=0.2$,
$D_{v}=1.65\times 10^{-6}$, $\Delta=10^{-4}$, $\Psi=0.94$ and $Pe_{v}=1$). At
this point, it is important to comment on the limitations of our numerical
model. As mentioned above, to maintain the thickness of the precursor layer
far from the drop and ensure accurate predictions for condensation and
evaporation, we employ a penalty function, see eq. (50). However, this
approach introduces a constraint on the freezing rate as the liquid layer
thickness approaches that of the precursor layer. Consequently, our model is
expected to fail at the end of the freezing process and thus cannot entirely
capture the formation of a cusp at the conclusion of freezing. With this
limitation in mind, and also considering that the cusp is expected to form
only at the very end of the freezing process, we can infer the tip angle by
the shape of the ice-gas interface in the vicinity of the tip. Thus, as shown
in figure 16 of Appendix §C, we provide an estimation for the tip angle by
drawing the tangents at the interface near the cusp and evaluate the tip angle
to be $\sim 145^{\circ}$, a prediction very close to the experimentally
reported value.

##  III Results and Discussion

###  III.1 Comparison with earlier studies - effect of evaporation

(a) (b)

(c) (d)

Figure 4: (a) Variation of (a) the evaporation flux $(J_{v})$ and (b) the
total mass of the condensate $(m_{c})$ along the substrate for different
values of $t/t_{f}$. Temporal variation of (c) the width of the frost halo
$(d_{halo})$ and (d) the left $(m_{cl})$ and right $(m_{cr})$ ends of the
condensation halo over the substrate. The maximum width of the halo is also
marked in panels (c) and (d). The rest of the dimensionless parameters are
$D_{w}=7.5$, $\epsilon=0.2$, $Ste=1.7\times 10^{-4}$, $T_{v}=0.2$,
$A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.82$, $\Lambda_{W}=191$,
$K=8\times 10^{-4}$, $Bi=2.29$ and $D_{v}=4.85\times 10^{-6}$. The value of
the dimensionless total freezing time $(t_{f})$ of the droplet is 1220,
respectively.

We begin our investigation by examining the consequences of including
evaporation in our model, which was neglected by the earlier theoretical
investigations [37, 9]. To this end, we consider a case with the following set
of parameters: $D_{w}=7.5$, $\epsilon=0.2$, $Ste=1.7\times 10^{-4}$,
$T_{v}=0.2$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.82$,
$\Lambda_{W}=191$, $K=8\times 10^{-4}$, $Bi=2.29$ and $D_{v}=4.85\times
10^{-6}$. This set of parameters corresponds to a system with similar physical
properties to that examined by Tembely and Dolatabadi [37]. While we observe
(not shown) that the freezing front and droplet shape with and without
evaporation obtained by setting $V_{\rho,r}=1.0$, $\chi=0$, $\Delta=0$,
$\Psi=0$ and $Pe_{v}=0$ (as considered by Tembely and Dolatabadi [37]) are
qualitatively the same, the droplet freezes faster when evaporation is
accounted for in the model. We also found that in the evaporation situation,
as opposed to when it is neglected, the freezing front propagates towards the
apex of the sessile droplet more quickly. For the set of parameters considered
in figure 4, the values of the dimensionless total freezing time $(t_{f})$ of
the droplet without and with evaporation are 1349 and 1220, respectively. This
behaviour can be attributed to evaporative cooling and the effect of latent
heat of evaporation mitigating the effect of latent heat of fusion in the
droplet. It is also to be noted that the model incorporating evaporation
exhibits a reduction in the expansion volume, albeit small, due to the
difference in densities between the ice and water phases, as compared to a
model without evaporation [37]. Note that our results cannot be directly
compared to Tembely and Dolatabadi [37] since our simulations correspond to a
thin ($\epsilon=0.2$) two-dimensional droplet, taking into account the
limitations of lubrication approximation, in contrast to a thick
($\epsilon=0.7$) axisymmetric droplet considered by Tembely and Dolatabadi
[37]. Another important difference is the fact that Tembely and Dolatabadi
[37] introduced a modified latent heat into their formulation. As reported
therein, it differs by almost a factor of three, with the normal latent heat,
significantly affecting the results of their simulations. It was reported that
this modification is necessary to achieve a good agreement with experiments.
Nevertheless, it is important to note that the introduction of the modified
latent heat is, in fact, incorrect and stems from an inconsistent use by these
authors of the reference temperature in the evaluation of the enthalpy jump at
the interface. As explained in §II.1.2, the enthalpy jump is given by eq.
(18). Considering that the local temperature at the liquid-solid interface is
very close to the melting temperature of ice, i.e. $T_{f}\approx T_{m}$ (even
when considering the Gibbs-Thomson effect as was done by Tembely and
Dolatabadi [37]), the effect of change in heat capacity between the two phases
in eq. (18) should actually be negligible. Additionally, Tembely and
Dolatabadi [37] account differently for the heat loss due to convection in the
gas phase, considering the dependence of the local gas temperature on the
droplet height. The latter effect has been neglected in our formulation. Both
these effects significantly impact the local temperature profile inside the
droplet, generating a super-cooling of the droplet, which affects the
predicted shape of the liquid-gas and liquid-ice interfaces, as will be
discussed in detail in the parametric study that follows. To account for the
super-cooling of the droplet, in figure 4, we consider a lower ambient
temperature $(T_{v}=0.2)$ than that of Tembely and Dolatabadi [37]. The effect
of $T_{v}$ is also investigated in §III.5.

The most distinct difference between the situations with and without
evaporation occurs during the early stages of freezing. It can be seen in
figure 4(a) that the droplet exhibits negative values of $J_{v}$ at the early
times in the region around the droplet. This indicates that the vapour
generated during the freezing process condenses on the substrate in the
vicinity of the contact line, forming a frost halo surrounding the freezing
drop. To determine the size and temporal evolution of the condensation halo,
we introduce the parameter $d_{halo}$, which represents the width of the frost
halo formed around the droplet. This width corresponds to the region where the
net condensate is accumulated, as determined numerically by integrating the
evaporation flux over time. Figure 4(b) displays the mass of the net
condensate near the droplet on the substrate at different instants. In this
figure, the area where ($m_{c}$) becomes negative indicates the presence of
the net condensate at that particular time. Moreover, figure 4(b) reveals that
initially, the condensate accumulates near the contact line ($x\approx 1$),
and as time progresses, this region moves away from the contact line.
Considering the diffusive nature of water vapor, we expect minimal
condensation far away from the droplet, which is also supported by our
simulations shown in figures 4(a) and 4(b), where the magnitude of the
negative evaporation flux and net condensate decreases as we move away from
the contact line. The interplay between condensation and evaporation processes
influences the location and width of the frost halo. The condensation halo
expands due to net condensation and contracts when net evaporation occurs. As
the mass of the condensate increases in a particular region of the substrate,
the condensate starts to re-evaporate, causing a shift in the area where the
evaporation flux is negative. Consequently, the region where the net
condensate is present changes as time progresses. In order to better
comprehend this phenomenon, we have plotted the temporal variation of
($d_{halo}$) in figure 4(c). The temporal variation of the leftmost ($m_{cl}$)
and rightmost ($m_{cr}$) ends of the net condensate are shown in figure 4(d).
It is evident from figure 4(c) that $d_{halo}$ initially increases, reaches a
maximum and then decreases before reaching the end of its lifetime. Moreover,
as shown by figure 4(d), the location of the frost halo also gradually shifts
away from the contact line over time, as shown by the evolution of the
leftmost ($m_{cl}$) and rightmost ($m_{cr}$) ends of frost halo. ($m_{cl}$)
and ($m_{cr}$) are initially close to each other, then move apart as
$d_{halo}$ grows, and finally come closer to each other towards the end of the
lifetime of the halo. It can be seen in figure 4(d) that, for the set of
parameters considered, the leftmost ($m_{cl}$) and rightmost ($m_{cr}$) ends
of the halo approach each other at ($x\approx 4.3$) just before eventually
disappearing. These results confirm that the droplet exhibits a frost halo
during the early stages of the freezing process $(0<t/t_{f}<0.014)$.

In the discussion that follows, we conduct a parametric investigation to
comprehend the underlying physics and effects of various parameters on the
formation of frost halo. Since the formation of the frost halo has been
studied in detail in the experimental work of Jung  _et al._ [6], we will
henceforth focus our attention on the range of parameters considered by these
authors and consider a ‘base’ case with the typical values of $Ste=1.22\times
10^{-3}$, $T_{v}=0$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.89$,
$\Lambda_{W}=0.33$, $\chi=0.01$, $K=8\times 10^{-4}$, $Bi=0.16$, $D_{w}=15$,
$D_{v}=1.65\times 10^{-6}$, $\Delta=10^{-4}$, $\Psi=0.94$ and $Pe_{v}=1$.
These parameters correspond to the freezing of a slender supercooled water
droplet at $-14.5^{\circ}\textrm{C}$ on a PMMA substrate and are kept
unchanged in the subsequent figures unless noted otherwise. Under the limit of
lubrication approximation, we consider a thin droplet with an equilibrium
contact angle of $14.3^{\circ}$. A similar system was examined in figure S5(A)
of Jung  _et al._ [6], albeit for a more hydrophobic system.

###  III.2 Effect of $RH$

In this section, we investigate the effect of relative humidity on the
freezing behaviour of the droplet by varying the value of the vapour density
ratio $V_{\rho,r}$. The relationship between the vapour density ratio
$V_{\rho,r}$ and relative humidity is given in eq. (39). It can be noted that
evaporation is only anticipated when the gas phase contains unsaturated humid
air. Figure 5 shows the profile of the vapour concentration along the liquid-
gas interface at different values of $t/t_{f}$ for $V_{\rho,r}=0.6$ and $0.85$
to examine the effect of vapour density ratio. The predicted total freezing
times $(t_{f})$ for $V_{\rho,r}=0.6$ and $0.85$ are 4700 and 5121,
respectively. We note that at very early times ($t/t_{f}=10^{-5}$), since the
droplet environment is not saturated, water evaporates with vapour
concentrating above the droplet; far from the droplet vapour concentration
remains equal to the vapour density ratio of the ambient gas phase. A local
maximum of the vapour concentration arises at the contact line since volatile
droplets evaporate preferentially. Later, vapour diffuses ahead of the
droplet, coming in contact with the supercooled precursor film ahead of the
droplet and condenses there, as shown below.

(a) (b)

Figure 5: Variation of the vapour concentration $(\rho_{v}\big{|}_{h})$ along
the substrate at different values of $t/t_{f}$ for (a) $V_{\rho,r}=0.60$ and
(b) $V_{\rho,r}=0.85$. The values of the rest of the dimensionless parameters
are $Ste=1.22\times 10^{-3}$, $T_{v}=0$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$,
$\Lambda_{S}=3.89$, $\Lambda_{W}=0.33$, $\chi=0.01$, $K=8\times 10^{-4}$,
$Bi=0.16$, $D_{w}=15$, $\epsilon=0.2$, $D_{v}=1.65\times 10^{-6}$,
$\Delta=10^{-4}$, $\Psi=0.94$ and $Pe_{v}=1$. The dotted line represents the
maximum vapour concentration attainable at the liquid-gas interface.

(a) (b)

(c) (d)

Figure 6: Variation of the evaporative flux $(J_{v})$ along the substrate for
(a) $V_{\rho,r}=0.35$, (b) $V_{\rho,r}=0.60$, (c) $V_{\rho,r}=0.85$ and (d)
$V_{\rho,r}=0.95$. The rest of the dimensionless parameters are the same as
figure 5. The dotted line is plotted to distinguish the positive and negative
values of the evaporation flux. The values of $t_{f}$ for $V_{\rho,r}=0.35$,
0.60, 0.85 and 0.95 are $t_{f}=4473$, 4700, 5121 and 5550, respectively.

(a) (b)

Figure 7: Temporal variation of (a) the width of the frost halo $(d_{halo})$
and (b) the left $(m_{cl})$ and right $(m_{cr})$ ends of the condensation halo
over the substrate for different values of the vapour density ratio
$(V_{\rho,r})$. In panel (b), the solid and dashed lines represent the
variation of $m_{cl}$ and $m_{cr}$, respectively. The rest of the
dimensionless parameters are the same as figure 5.

To examine the effect of vapour density ratio on the evaporation/condensation
rate, in figure 6, we depict the profile of the evaporation flux at different
values of $t/t_{f}$ for $V_{\rho,r}=0.35-0.95$. It can be seen that increasing
the value of $V_{\rho,r}$ considerably decreases the evaporation flux
$(J_{v})$ at the early times $(t/t_{f}<10^{-3})$. It is clear from eq. (21)
that if the vapour concentration at any point along the liquid-gas interface
exceeds saturation, the evaporation flux becomes negative, indicating the
condensation of vapour on the substrate. The opposite happens (evaporation
occurs) when the vapour concentration falls below the saturation value. For
$V_{\rho,r}=0.35$, no condensation takes place, and water continues to
evaporate while the droplet undergoes freezing (figure 6a). The evaporation
flux is highest at the apex of the droplet, with a local maximum arising in
the contact line region of the droplet. This is evident at $t/t_{f}=10^{-5}$
in figure 6a. As expected, the evaporation flux at the apex of the droplet
decreases as time progresses. For $V_{\rho,r}=0.6$, we observe condensation
near the contact line of the droplet at $t/t_{f}=10^{-5}$ (figure 6a).
However, subsequently (for $t/t_{f}\geq 10^{-3}$), the evaporation flux turns
positive along the entire substrate, indicating that the condensation is
eventually suppressed. Inspection of figure 6(a-d) reveals that the
condensation in the contact line region is enhanced with increasing humidity
of the droplet atmosphere. Nevertheless, the increased humidity also leads to
significantly lower evaporation rates. In fact, for very humid environments
(see figure 6d for $V_{\rho,r}=0.95$), evaporation takes place only in a small
part of the droplet (i.e. near the contact line), and vapour mostly condenses
along the liquid-gas interface. We also note that in the early stages of the
freezing process, the relative importance between evaporation and condensation
changes drastically with increasing humidity.

As discussed in §III.1, the vapour condensation in the contact line region
forms the frost halo surrounding the frozen droplet. As shown in figure 6, the
droplet exhibits a frost halo for $V_{\rho,r}\geq 0.6$ for the base
parameters. Thus, to evaluate the effect of the vapour density ratio on the
estimated size of halo $(d_{halo})$, the temporal variation of the thickness
of the region beyond the contact line with net condensate is plotted in figure
7(a) for different values of vapour density ratio ranging from
$V_{\rho,r}=0.60$ to $V_{\rho,r}=0.95$. The corresponding temporal variation
of the left $(m_{cl})$ and right $(m_{cr})$ ends of the condensation halo over
the substrate are shown in figure 7(b). It can be seen that the width of the
frost halo, $d_{halo}$, initially increases with time for all the values of
$V_{\rho,r}$, reaches a maximum before disappearing at later stages of the
freezing process. Inspection of these results reveals that an increase in the
vapour density ratio $V_{\rho,r}$ has a significant impact on the total
lifetime ($t/t_{f}$) of the frost halo. This effect is primarily attributed to
a substantial reduction in evaporation, which prolongs the presence of the net
condensate on the substrate. Figure 7(b) further illustrates that as the value
of $V_{\rho,r}$ increases, the extent of the halo increases significantly.

At this point, it would be useful to compare the evolution of the freezing
front obtained from our model with the experimental observations of Sebilleau
_et al._ [44]. In the latter study, the freezing of a water drop on a copper
substrate was investigated under the conditions of a substrate temperature of
$-15^{\circ}$C and an ambient temperature ranging from $15^{\circ}$C to
$20^{\circ}$C. Their findings revealed a substantial impact of relative
humidity on the freezing front propagation rates, while wettability did not
significantly influence the process. On the other hand, our study explores a
considerably colder ambient temperature range ($-15^{\circ}$C to
$0^{\circ}$C). It is important to note that the freezing front dynamics depend
on the heat transfer from the water-ice interface to the ice, substrate, and
the remaining water and atmosphere, as suggested by Jung  _et al._ [6],
Sebilleau  _et al._ [44]. Thus, in contrast to our investigation, Sebilleau
_et al._ [44] emphasizes the substantial impact of evaporative cooling,
leading to a significant variation in the freezing front propagation rate
under different relative humidity conditions. By performing a simulation for a
scenario where the ambient temperature aligns with the melting temperature
($T_{v}=1$) and a smaller droplet with a contact angle of $14.3^{\circ}$,
while keeping the remaining parameters the same as those of Sebilleau  _et
al._ [44], we observe a close match agreement between our simulated freezing
front propagation and experimental result of Sebilleau  _et al._ [44] at
$40\%$ relative humidity (see figure 8a). However, given that the heat flux
primarily traverses into the ice and subsequently to the substrate, we do not
observe significant variations in the freezing front propagation rates for
different relative humidity values. Furthermore, we explore the impact of the
contact angle of the droplet in figure 8(b) for a typical set of parameters.
Our findings indicate that changes in wettability do not significantly affect
the freezing front propagation. A similar behaviour was observed by Sebilleau
_et al._ [44].

(a) (b)

Figure 8: (a) Comparison of the evolution of the freeing front obtained from
our simulation with the experimental result of Sebilleau  _et al._ [44]. Here,
the results shown by lines are from our simulations with different values of
the relative humidity for a droplet with a contact angle of $14.3^{\circ}$.
The symbols are associated with the experimental results of Sebilleau  _et
al._ [44]. (b) Variation of the freezing front with $(\alpha t)^{1/2}$ for
different values of the contact angle for $RH=0.7$. The values of the
dimensionless parameters are $Ste=1.22\times 10^{-3}$, $T_{v}=1$,
$A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.89$, $\Lambda_{W}=698$,
$\chi=0.01$, $K=8\times 10^{-4}$, $Bi=0.16$, $D_{w}=15$, $V_{\rho,r}=0.70$,
$\epsilon=0.2$, $D_{v}=1.65\times 10^{-6}$, $\Delta=10^{-4}$, $\Psi=0.94$ and
$Pe_{v}=1$.

###  III.3 Effect of liquid volatility

Figure 9: Evolution of the freezing front, $s$ (dot-dashed lines) and shape of
the droplet, $h$ (solid line) for (a) $D_{v}=10^{-4}$ and (b) $D_{v}=10^{-2}$.
The values of the rest of the dimensionless parameters are $Ste=1.22\times
10^{-3}$, $T_{v}=0$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.89$,
$\Lambda_{W}=0.33$, $\chi=0.01$, $V_{\rho,r}=0.85$, $K=8\times 10^{-4}$,
$Bi=0.16$, $D_{w}=15$, $\epsilon=0.2$, $\Delta=10^{-4}$, $\Psi=0.94$ and
$Pe_{v}=1$. The values of $t_{f}$ for $D_{v}=10^{-4}$ and $10^{-2}$ are
$t_{f}=4210$ and 634, respectively.

(a) (b)

(c) (d)

(e)

Figure 10: Variation of the evaporative flux $(J_{v})$ along the substrate for
(a) $D_{v}=10^{-4}$ and (b) $D_{v}=10^{-2}$. Temporal variation of (c) the
width of the frost halo $(d_{halo})$, (d) the left $(m_{cl})$ and right
$(m_{cr})$ ends of the condensation halo over the substrate, and (e) the
location of maximum condensation flux $(x_{halo})$ for different values of
$D_{v}$. In panel (d), the solid and dashed lines represent the variation of
$m_{cl}$ and $m_{cr}$, respectively. The rest of the dimensionless parameters
are the same as figure 9. The values of $t_{f}$ for $D_{v}=10^{-4}$, $10^{-3}$
and $10^{-2}$ are $t_{f}=4210$, 2318 and 634, respectively.

To further examine the effect of liquid volatility, we vary the parameter,
$D_{v}$, which is the ratio between the equilibrium vapour concentration and
liquid density, in figure 9. We consider a wide range for this parameter (not
applicable to a typical water droplet) in order to examine cases covering the
entire range between non-volatile and highly volatile liquids. A few examples
of highly volatile liquids are butane, toluene and octane. The total freezing
time for highly volatile liquids decreases drastically ($t_{f}=634$ for
$D_{v}=10^{-2}$ and $t_{f}=4210$ for $D_{v}=10^{-4}$) due to both the effect
of evaporative cooling as well as the faster depletion of the liquid layer
which vaporizes in the presence of an unsaturated atmosphere. In addition, the
difference in liquid volatility impacts the final shape of the frozen droplet.
Highly volatile droplets assume a much thinner profile with a lower total mass
of ice than less volatile ones.

To examine the effect of $D_{v}$ on the frost halo formation, we plot, in
figure 10(a) and (b), the profiles of the evaporation flux ($J_{v}$) for
$D_{v}=10^{-4}$ and $10^{-2}$ at various time instants. As expected, the
evaporation flux is higher for highly volatile liquids, but we also notice
more substantial condensation in the periphery of the droplet, driven by a
higher amount of vapour in that area. In order to quantify this effect, in
figure 10(c) and (d), we plot the temporal variation of the width of the frost
halo ($d_{halo}$), the leftmost point $(m_{cl})$ and the rightmost point of
the condensation halo $(m_{cr})$ for different values of $D_{v}$. It can be
seen that both the normalised lifespan ($t/t_{f}$) and the width of the halo
increase considerably with the increasing volatility of the liquid. Moreover,
we notice that the condensation occurs closer to the droplet with increasing
volatility (see figure 10(a) and (b)). This is shown more clearly in figure
10(e), where we plot the position of the maximum condensation flux
$(x_{halo})$ for different values of $D_{v}$. Evaporation is accompanied by
cooling due to the effect of latent heat. We observe (not shown) that the
concavity of the freezing front increases with increasing $\chi$ due to the
increased influence of evaporative cooling near the tip of the droplet. This
leads to a local decrease in the temperature, promoting freezing in that
region. The total freezing times for $\chi$ to $0.01$, $0.1$, and $1.0$ are
$t_{f}=4210$, 4110, and 3503, respectively.

###  III.4 Effect of substrate thermal resistance

Figure 11: Evolution of the freezing front, $s$ (dot-dashed lines) and shape
of the droplet, $h$ (solid line) for (a) $\Lambda_{W}=1$ and (b)
$\Lambda_{W}=500$. The values of the rest of the dimensionless parameters are
$Ste=1.22\times 10^{-3}$, $T_{v}=0$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$,
$\chi=0.01$, $\Lambda_{S}=3.89$, $V_{\rho,r}=0.85$, $K=8\times 10^{-4}$,
$Bi=0.16$, $D_{w}=15$, $\epsilon=0.2$, $D_{v}=1.65\times 10^{-6}$,
$\Delta=10^{-4}$, $\Psi=0.94$ and $Pe_{v}=1$.

(a) (b)

(c) (d)

(e)

Figure 12: Variation of the evaporative flux along the substrate for (a)
$\Lambda_{W}=1$ and (b) $\Lambda_{W}=500$. Temporal variation of (c) the width
of the frost halo $(d_{halo})$, (d) the left $(m_{cl})$ and right $(m_{cr})$
ends of the condensation halo over the substrate, and (e) the location of
maximum condensation flux $(x_{halo})$ for different values of $\Lambda_{W}$.
In panel (d), the solid and dashed lines represent the variation of $m_{cl}$
and $m_{cr}$, respectively. The rest of the dimensionless parameters are the
same as figure 11. The values of $t_{f}$ for $\Lambda_{W}=1$, $10$, $100$ and
$500$ are $t_{f}=4240$, 1400, 352 and 240, respectively. Figure 13: Variation
of the mass of ice layer, $m_{ice}$ with ${\rm{log}}(t/t_{f})$ for different
values of $\Lambda_{W}$. The rest of the dimensionless parameters are the same
as figure 11.

Next, we investigate the effect of substrate thermal resistance on the
freezing behavior of drop by changing the thermal conductivity of the
substrate, for a constant thickness of the substrate and keeping
$\lambda_{l}=0.57$ Wm-1K-1 fixed. Here, $V_{\rho,r}=0.85$ and the rest of the
parameters are the same as those in figure 5. The temporal evolutions of the
freezing front and the droplet shape (figures 11) reveal that increasing the
thermal conductivity of the substrate promotes freezing and influences the
formation of the freezing front by making the final droplet shape more cusp-
like morphology. A close inspection of figure 11(a) and (b) reveals that the
droplet exhibits a concave freezing front for a low thermal conductivity
substrate. This is because, for low conductive substrates, the liquid-gas
interface dissipates more heat through convection and evaporative cooling than
the substrate itself, resulting in a lower temperature near the triple contact
line than at the centre of the droplet. Quantitatively, we observe that the
total freezing time of the droplet for $\Lambda_{W}=1$ is almost 18 times that
for $\Lambda_{W}=500$. A similar result was also observed in the experiments
of Jung  _et al._ [6].

To examine the effect of the thermal conductivity ratio on the frost halo
formation, the variation of the evaporation flux $(J_{v})$ along the substrate
for different values of $t/t_{f}$ are plotted in figure 12(a) and (b) for
$\Lambda_{W}=1$ and 500, respectively. In order to quantify this effect, in
figure 12(c), we plot the temporal variation of the width of the frost halo
($d_{halo}$), the leftmost ($m_{cl}$) and rightmost point ($m_{cr}$) of the
net condensate for different values of $\Lambda_{W}$. Figure 12(c) shows that,
as the thermal conductivity of the substrate ($\Lambda_{W}$) increases, the
width of the halo decreases. It is also evident from figure 12(d) that the
halo region is closer to the contact line for the substrate with higher
conductivity. As the freezing front evolves, the temperature at the liquid-gas
interface increases, leading to a slight rise in evaporation for the substrate
with higher thermal conductivity. It is essential to note that the total
freezing time ($t_{f}$) increases as the thermal conductivity of the substrate
decreases. This leads to a significantly higher amount of vaporized liquid
during the whole freezing process for the droplet on a substrate with a low
thermal conductivity than the higher thermal conductivity (see figure 13). Due
to this evaporation, the resulting vapour condenses in the droplet’s vicinity,
with its concentration reaching supersaturation. In figure 12(a), it can be
observed that the condensation occurs far away from the droplet (near $x=2$)
at $t/t_{f}=10^{-3}$. In contrast, it occurs near the contact line of the
droplet in figure 12(b). To better quantify this effect, in figure 12(e), we
plot the position of the maximum condensation flux $(x_{halo})$ for different
values of $\Lambda_{w}$. In our configuration, as the condensation happens at
an earlier time of the freezing, the scaled time $(t/t_{f})$ for the substrate
with lower conductivity represents a higher value in terms of the absolute
dimensionless time $(t)$ when compared to the substrate with higher
conductivity leading to condensation far away from the substrate. Regarding
the effect of the substrate thickness, it can be noted that numerically
$D_{w}=0$ resembles a substrate with infinite thermal conductivity. Increasing
the thickness of the substrate increases the thermal resistance of the
substrate, and thereby, it has the same effect as decreasing the thermal
conductivity of the substrate.

###  III.5 Effect of $T_{v}$

(a) (b)

(c) (d)

(e)
The values of $t_{f}$ for $T_{v}=0$, $0.5$ and $1$ are $t_{f}=5121$, 7919 and
19060, respectively.

Figure 14: Variation of the evaporative flux along the substrate for (a)
$T_{v}=0$ and (b) $T_{v}=1$. Temporal variation of (c) the width of the frost
halo $(d_{halo})$, (d) the left $(m_{cl})$ and right $(m_{cr})$ ends of the
condensation halo over the substrate, and (e) the location of maximum
condensation flux $(x_{halo})$ for different values of $T_{v}$. In panel (d),
the solid and dashed lines represent the variation of $m_{cl}$ and $m_{cr}$,
respectively. The values of the rest of the dimensionless parameters are
$Ste=1.22\times 10^{-3}$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$,
$\Lambda_{S}=3.89$, $\Lambda_{W}=0.33$, $\chi=0.01$, $V_{\rho,r}=0.85$,
$K=8\times 10^{-4}$, $Bi=0.16$, $D_{w}=15$, $\epsilon=0.2$, $D_{v}=1.65\times
10^{-6}$, $\Delta=10^{-4}$, $\Psi=0.94$ and $Pe_{v}=1$.

Finally, we investigate the effect of the ambient temperature ($T_{v}$) on the
freezing behaviour of the droplet. Here, $V_{\rho,r}=0.85$ and the rest of the
parameters are the same as those used to generate figure 5 (‘base’
parameters). Here, we consider a thick substrate with low thermal conductivity
(i.e. $D_{w}=15$ and $\Lambda_{W}=0.33$), resulting in less heat conduction
through the substrate. It can be seen in figure 14(a) that the evaporation is
suppressed except for a region near the contact line where the evaporation
flux is maximized. Due to the effect of latent heat, the local temperature
decreases, favouring freezing near the contact line and leading to a concave
shape of the liquid-ice interface. In contrast, when the ambient temperature
is higher (e.g. for $T_{v}=1$, the gas temperature is equal to the melting
point temperature), the convective heat transfer occurs from the ambient to
the droplet, raising the temperature of the remaining liquid in the droplet.
Besides delaying the freezing process, this results in the enhancement of
vaporization (see figure 14(e)) and leads to a more uniform evaporation flux
along the liquid-gas interface, as shown in figure 14(b). This, in turn,
results in a nearly flat profile of the liquid-ice interface. Figure 14(c) and
(d) depict the temporal variation of the width of the frost halo $(d_{halo})$,
and the left $(m_{cl})$ and right $(m_{cr})$ ends of the condensation halo
over the substrate for different values of $T_{v}$. It can be seen that, as
$T_{v}$ increases, the timescale of the halo formation decreases even though
the maximum width of the halo and the extent of the halo do not change
considerably since the evaporation along the contact line does not change
significantly for different values of $T_{v}$.

##  IV Conclusions

We have numerically investigated the freezing behaviour of a sessile liquid
droplet placed on a cold, horizontal solid substrate using the lubrication
approximation in the finite element method framework. We present a
comprehensive model for freezing front dynamics, encompassing solid-liquid
phase change, heat balance across the liquid-gas interface, heat flux between
the cold substrate and ice, and mass conservation involving the liquid, solid,
and gas phases. The current work mainly focuses on the frost halo formation,
which has only been observed experimentally [6]. We found that evaporation has
a distinct effect on the freezing process of a sessile droplet, in contrast to
when evaporation is neglected [37, 9]. In the initial freezing stages, we
observe a negative evaporation flux surrounding the droplet, indicating that
vapor generated during freezing condenses on the substrate near the contact
line, forming a frost halo. The accumulation of condensate triggers re-
evaporation, causing a temporal shift of the frost halo region away from the
contact line and ultimately, the halo disappears due to the diffusive nature
of water vapor far from the droplet. We found that an increase in relative
humidity extends the lifetime of the frost halo, as it significantly reduces
evaporation, thereby prolonging the presence of net condensate on the
substrate.

We also performed a parametric study to examine the effect of the thermal
resistance of the substrate and ambient temperature on the freezing behaviour
of the droplet. It is observed that greater liquid volatility results in an
increased evaporation flux, and condensation occurs closer to the droplet as a
higher concentration of vapor is present in the surrounding region of the
droplet. Additionally, decreasing the thermal conductivity of the substrate
extends the total freezing time and increases the amount of vaporized liquid.
This, in turn, results in condensation, with the concentration reaching
supersaturation. As expected, reducing the thermal resistance of the
substrate, by increasing the thermal conductivity or reducing the substrate
thickness, promotes freezing. We found another interesting behaviour - the
droplet exhibits a concave freezing front for a low thermal conductivity
substrate. This is because, in the case of a less-conductive substrate, the
liquid-gas interface dissipates more heat through convection and evaporative
cooling than the substrate itself, resulting in a lower temperature near the
triple contact line than at the centre of the droplet.

Acknowledgement: K.C.S. thanks the Science & Engineering Research Board and
IIT Hyderabad, India for their financial supports through grants
CRG/2020/000507 and IITH/CHE/F011/SOCH1, respectively. C.S.S. acknowledges IIT
Ropar for financial support through ISIRD (Grant No. 9-388/2018/IITRPR/3335).
We also acknowledge Vineet Kumar (IIT Ropar) for preliminary results on
sessile droplet freezing in figure 1.

Declaration of interests: The authors report no conflict of interest.

##  Appendix A Lubrication model

By employing the scalings shown in Eq. (27) and incorporating the lubrication
approximation ($\epsilon\ll 1$), we obtain the following dimensionless
governing equations and boundary conditions in the liquid, ice and gas phases.

####  A.0.1 Liquid phase

The dimensionless governing equations in the liquid phase are given by

$\displaystyle\partial^{2}_{z}u$ $\displaystyle=$
$\displaystyle\partial_{x}p,$ (A.1) $\displaystyle\partial_{z}p$
$\displaystyle=$ $\displaystyle 0,$ (A.2)
$\displaystyle\partial_{x}u+\partial_{z}w$ $\displaystyle=$ $\displaystyle 0,$
(A.3) $\displaystyle\partial^{2}_{z}T_{l}$ $\displaystyle=$ $\displaystyle 0.$
(A.4)

The boundary conditions at the liquid-gas interface ($z=h(x,t)$) are

$\displaystyle p$ $\displaystyle=$ $\displaystyle-\kappa_{lg}\gamma_{lg}-\Pi,$
(A.5) $\displaystyle\partial_{z}u$ $\displaystyle=$ $\displaystyle 0,$ (A.6)
$\displaystyle\partial_{z}T_{l}$ $\displaystyle=$ $\displaystyle-\chi
J_{v}+Bi\left(T_{v}-T_{l}\right),$ (A.7)
$\displaystyle\partial_{t}h+u\partial_{x}h-w$ $\displaystyle=$ $\displaystyle-
D_{v}J_{v}.$ (A.8)

Here, $D_{v}=\rho_{ve}(T_{m})/\rho_{l}$ represents the density ratio,
$\chi=\epsilon\rho_{ve}(T_{m})UL_{v}H/(\lambda_{l}\Delta T)$ denotes the
scaled latent heat of vaporization and $Bi=h_{tc}H/\lambda_{l}$ is the Biot
number.

The boundary conditions at the liquid-ice interface (at $z=s(x,t)$) are

$u=0\;\;{\rm and}\;\;T_{l}=T_{f},$ (A.9)

where $T_{f}=T_{m}$.

The dimensionless disjoining pressure is given by

$\Pi=\epsilon^{-2}A_{n}\left[\left(\frac{\beta}{h-s}\right)^{n}-\left(\frac{\beta}{h-s}\right)^{m}\right],$
(A.10)

where $\beta$ is of the same order as the equilibrium precursor layer
thickness [50] and $A_{n}=HA/\gamma_{lgo}$ denotes the dimensionless Hamaker
constant. The interaction between the repulsive and attractive components of
eq. (A.10) dictates the value of the equilibrium contact angle $\theta_{eq}$
[38, 9, 51, 50, 37]. This can be approximated [50]

$\theta_{eq}\approx\sqrt{\beta A_{n}}.$ (A.11)

The full mean curvatures at the liquid-gas and liquid-ice interfaces,
retaining higher order contributions, are given by

$\kappa_{lg}=\frac{h_{xx}}{(1+\epsilon^{2}h_{x}^{2})^{\frac{3}{2}}}~{}{\rm
and}~{}\kappa_{sl}=\frac{s_{xx}}{(1+\epsilon^{2}s_{x}^{2})^{\frac{3}{2}}}.$
(A.12)

####  A.0.2 Solid (ice) phase

Under the lubrication approximation ($\epsilon\ll 1$), the dimensionless
energy conservation equations for the ice phase is given by

$\partial^{2}_{z}T_{s}=0.$ (A.13)

The energy balance at the liquid-ice interface ($z=s(x,t)$) gives

$Ste\left(\Lambda_{s}\partial_{z}T_{s}-\partial_{z}T_{l}\right)=J_{s},$ (A.14)

where $\Lambda_{s}=\lambda_{s}/\lambda_{l}$ is the thermal conductivity ratio
and $Ste=\lambda_{l}\Delta T/(\epsilon\rho_{s}UL_{f}H)$ denotes the Stefan
number. At $z=s(x,t)$, the conservation of mass (eq. 10) gives

$\partial_{t}s-w=D_{s}J_{s},~{}{\rm where}~{}J_{s}=\partial_{t}s.$ (A.15)

At the ice-solid interface ($z=0$), we impose

$T_{s}=T_{w}.$ (A.16)

####  A.0.3 Gas phase

The dimensionless conservation equation for the vapor concentration becomes

$Pe_{v}({\partial_{t}\rho_{v}}+u_{g}\partial_{x}\rho_{v}+w_{g}\partial_{z}\rho_{v})=\partial^{2}_{z}\rho_{v}+\epsilon^{2}\partial^{2}_{x}\rho_{v}.$
(A.17)

At the liquid-gas interface ($z=h(x,t)$), the conservation of mass reduces to

$Pe_{v}J_{v}=-\partial_{z}\rho_{v},$ (A.18)

while the constitutive equation for the evaporation flux gives

$KJ_{v}=1-\rho_{v}+\Delta p+\Psi\left(T_{l}-1\right).$ (A.19)

The boundary condition far from the droplet ($z=D_{g}$) is given by

$\rho_{v}=V_{\rho,r}.$ (A.20)

The various dimensionless numbers appearing in eqs. (A.17 \- A.19) are

$\begin{gathered}Pe_{v}=\frac{\epsilon UH}{D_{m}},~{}K=\frac{\epsilon
U}{\alpha}\left(\frac{2\pi
M}{R_{g}T_{m}}\right)^{1/2},~{}\Delta=\frac{\mu_{l}ULM}{H^{2}\rho_{l}R_{g}T_{m}},\\\
\Psi=\frac{L_{v}M\Delta T}{R_{g}T_{m}^{2}},~{}{\rm
and}~{}V_{\rho,r}={\rho_{vi}\over\rho_{ve}(T_{m})}.\end{gathered}$ (A.21)

####  A.0.4 Solid substrate

The scaled energy conservation equation for the solid substrate is given by

$\partial^{2}_{z}T_{w}=0.$ (A.22)

The boundary condition at the solid-ice interface ($z=0$) is

$\Lambda_{s}\partial_{z}T_{s}=\Lambda_{w}\partial_{z}T_{w},$ (A.23)

where $\Lambda_{w}=\lambda_{w}/\lambda_{l}$ denotes the thermal conductivity
ratio. At the bottom of the substrate (at $z=-D_{w}$), we impose

$T_{w}=0.$ (A.24)

###  A.1 Evolution equations

By integrating eqs. (A.1) and (A.2) with respect to $z$ and using eqs. (A.6),
(A.9) and (A.5), we get

$u=\frac{\partial_{x}p}{2}(z^{2}-s^{2})-h\partial_{x}p(z-s),$ (A.25)
$p=-\kappa_{lg}\gamma_{lg}-\Pi.$ (A.26)

By integrating eq. (A.3) and using eq. (A.8), we get the following evolution
equation:

$\partial_{t}h-\partial_{t}s=-\partial_{x}q_{l}-D_{v}J_{v}-D_{s}J_{s},$ (A.27)

where

$q_{l}=\frac{\partial_{x}p}{2}\left(\frac{h^{3}}{3}-s^{2}h+\frac{2s^{3}}{3}\right)-h\partial_{x}p\left(\frac{h^{2}}{2}-sh+\frac{s^{2}}{2}\right).$
(A.28)

Similarly, by integrating eq. (A.4) and using eqs. (A.7) and (A.9), we get

$T_{l}=(-\chi J_{v}+Bi(T_{v}-T_{l}\big{|}_{h}))(z-s)+T_{f}.$ (A.29)

The temperature distribution in the ice phase is governed by

$T_{s}=\frac{T_{f}}{D_{w}+s\Lambda_{w}/\Lambda_{s}}\left(D_{w}+z\frac{\Lambda_{w}}{\Lambda_{s}}\right).$
(A.30)

Using the above expression along with eq. (A.15) and introducing them into eq.
(A.14), we get

$\partial_{t}s=Ste\left(\frac{\Lambda_{w}T_{f}}{D_{w}+s\Lambda_{w}/\Lambda_{s}}+\chi
J_{v}-Bi(T_{v}-T_{l}\big{|}_{h})\right).$ (A.31)

The temperature profile in the solid substrate is given by

$T_{w}=\frac{T_{f}}{D_{w}+s\;\Lambda_{w}/\Lambda_{s}}(z+D_{w}).$ (A.32)

##  Appendix B Grid independence test

Figure 15: Evolution of the freezing front, $s$ (dot-dashed lines) and
magnified shape of the upper part of the droplet, $h$ (solid line) placed on a
cold substrate. The results obtained using 2001, 4801, and 8001 grid points
were indistinguishable. The rest of the dimensionless parameters are
$V_{\rho,r}=0.85$, $\chi=0.23$, $\Delta=10^{-4}$, $\Psi=0.14$ and
$Pe_{v}=1.0$, $D_{w}=7.5$, $\epsilon=0.2$, $Ste=1.7\times 10^{-4}$,
$T_{v}=0.2$, $A_{n}=6.25$, $D_{g}=2$, $D_{s}=0.9$, $\Lambda_{S}=3.82$,
$\Lambda_{W}=191$, $K=8\times 10^{-4}$, $Bi=2.29$ and $D_{v}=4.85\times
10^{-6}$. The values of the dimensionless total freezing time $(t_{f})$ of the
droplet is 1220, respectively.

##  Appendix C Tip angle

Figure 16: Figure showing the estimation of the tip angle for a typical set of
parameters ($Ste=1.22\times 10^{-3}$, $T_{v}=1$, $A_{n}=6.25$, $D_{g}=2$,
$D_{s}=0.9$, $\Lambda_{S}=3.89$, $\Lambda_{W}=698$, $\chi=0.01$, $K=8\times
10^{-4}$, $Bi=0.16$, $D_{w}=15$, $V_{\rho,r}=0.70$, $\epsilon=0.2$,
$D_{v}=1.65\times 10^{-6}$, $\Delta=10^{-4}$, $\Psi=0.94$ and $Pe_{v}=1$).

## References

  * Myers and Thompson [1998] T. G. Myers and C. P. Thompson, Modeling the flow of water on aircraft in icing conditions, AIAA J. 36, 1010 (1998).
  * Cao  _et al._ [2015] Y. Cao, Z. Wu, Y. Su, and Z. Xu, Aircraft flight characteristics in icing conditions, Prog. Aerosp. Sci. 74, 62 (2015).
  * Mensah and Choi [2015] K. Mensah and J. M. Choi, Review of technologies for snow melting systems, J. Mech. Sci. Technol. 29, 5507 (2015).
  * Pruppacher and Klett [1980] H. R. Pruppacher and J. D. Klett, Microphysics of clouds and precipitation, Nature 284, 88 (1980).
  * Angell [1983] C. A. Angell, Supercooled water, Annu. Rev. Phys. Chem. 34, 593 (1983).
  * Jung  _et al._ [2012a] S. Jung, M. K. Tiwari, and D. Poulikakos, Frost halos from supercooled water droplets, Proc. Natl. Acad. Sci. 109, 16073 (2012a).
  * Hu and Jin [2010] H. Hu and Z. Jin, An icing physics study by using lifetime-based molecular tagging thermometry technique, Int. J. Multiph. Flow 36, 672 (2010).
  * Anderson  _et al._ [1996] D. M. Anderson, M. G. Worster, and S. H. Davis, The case for a dynamic contact angle in containerless solidification, J. Cryst. Growth. 163, 329 (1996).
  * Zadravzil  _et al._ [2006] A. Zadravzil, F. Stepanek, and O. K. Matar, Droplet spreading, imbibition and solidification on porous media, J. Fluid Mech. 562, 1 (2006).
  * Sanz  _et al._ [1987] A. Sanz, J. Meseguer, and L. Mayo, The influence of gravity on the solidification of a drop, J. Cryst. Growth. 82, 81 (1987).
  * Enríquez  _et al._ [2012] O. R. Enríquez, Á. G. Marín, K. G. Winkels, and J. H. Snoeijer, Freezing singularities in water drops, Phys. Fluids 24, 091102 (2012).
  * Marin  _et al._ [2014] A. G. Marin, O. R. Enriquez, P. Brunet, P. Colinet, and J. H. Snoeijer, Universality of tip singularity formation in freezing water drops, Phys. Rev. Lett. 113, 054301 (2014).
  * Peng  _et al._ [2020] H. Peng, Q. Wang, T. Wang, L. Li, Z. Xia, J. Du, B. Zheng, H. Zhou, and L. Ye, Study on dynamics and freezing behaviors of water droplet on superhydrophobic aluminum surface, Appl. Phys. A 126, 1 (2020).
  * Pan  _et al._ [2019] Y. Pan, K. Shi, X. Duan, and G. F. Naterer, Experimental investigation of water droplet impact and freezing on micropatterned stainless steel surfaces with varying wettabilities, Int. J. Heat Mass Transf. 129, 953 (2019).
  * Liu  _et al._ [2021] X. Liu, J. Min, X. Zhang, Z. Hu, and X. Wu, Supercooled water droplet impacting-freezing behaviors on cold superhydrophobic spheres, Int. J. Multiph. Flow 141, 103675 (2021).
  * Jin  _et al._ [2017] Z. Jin, H. Zhang, and Z. Yang, The impact and freezing processes of a water droplet on different cold cylindrical surfaces, Int. J. Heat Mass Transf. 113, 318 (2017).
  * Zhang  _et al._ [2018] H. Zhang, Z. Jin, M. Jiao, and Z. Yang, Experimental investigation of the impact and freezing processes of a water droplet on different cold concave surfaces, Int. J. Therm. Sci. 132, 498 (2018).
  * Ju  _et al._ [2018] J. Ju, Z. Jin, H. Zhang, Z. Yang, and J. Zhang, The impact and freezing processes of a water droplet on different cold spherical surfaces, Exp. Therm. Fluid Sci. 96, 430 (2018).
  * Jung  _et al._ [2011] S. Jung, M. Dorrestijn, A. Raps, D.and Das, C. M. Megaridis, and D. Poulikakos, Are superhydrophobic surfaces best for icephobicity?, Langmuir 27, 3059 (2011).
  * Huang  _et al._ [2012] L. Huang, Z. Liu, Y. Liu, Y. Gou, and L. Wang, Effect of contact angle on water droplet freezing process on a cold flat surface, Exp. Therm. Fluid Sci. 40, 74 (2012).
  * Jung  _et al._ [2012b] S. Jung, M. K. Tiwari, N. V. Doan, and D. Poulikakos, Mechanism of supercooled droplet freezing on surfaces, Nat. Commun. 3, 1 (2012b).
  * Boinovich  _et al._ [2014] L. Boinovich, A. M. Emelyanenko, V. V. Korolev, and A. S. Pashinin, Effect of wettability on sessile drop freezing: when superhydrophobicity stimulates an extreme freezing delay, Langmuir 30, 1659 (2014).
  * Hao  _et al._ [2014] P. Hao, C. Lv, and X. Zhang, Freezing of sessile water droplets on surfaces with various roughness and wettability, Appl. Phys. Lett. 104, 161609 (2014).
  * Chaudhary and Li [2014] G. Chaudhary and R. Li, Freezing of water droplets on solid surfaces: An experimental and numerical study, Exp. Therm. Fluid Sci. 57, 86 (2014).
  * Zhang  _et al._ [2016] H. Zhang, Y. Zhao, R. Lv, and C. Yang, Freezing of sessile water droplet for various contact angles, Int J Therm Sci 101, 59 (2016).
  * Nath  _et al._ [2017] S. Nath, S. F. Ahmadi, and J. B. Boreyko, A review of condensation frosting, Nanoscale Microscale Thermophys. Eng. 21, 81 (2017).
  * Yancheshme  _et al._ [2020] A. A. Yancheshme, G. Momen, and R. J. Aminabadi, Mechanisms of ice formation and propagation on superhydrophobic surfaces: A review, Adv. Colloid Interface Sci. 279, 102155 (2020).
  * Castillo  _et al._ [2021] J. E. Castillo, Y. Huang, Z. Pan, and J. A. Weibel, Quantifying the pathways of latent heat dissipation during droplet freezing on cooled substrates, Int. J. Heat Mass Transf. 164, 120608 (2021).
  * Graeber  _et al._ [2018] G. Graeber, V. Dolder, T. M. Schutzius, and D. Poulikakos, Cascade freezing of supercooled water droplet collectives, ACS Nano 12, 11274 (2018).
  * Vu  _et al._ [2015] T. Vu, G. Tryggvason, S. Homma, and J. C. Wells, Numerical investigations of drop solidification on a cold plate in the presence of volume change, Int. J. Multiph. Flow 76, 73 (2015).
  * Vu  _et al._ [2018] T. V. Vu, K. V. Dao, and B. D. Pham, Numerical simulation of the freezing process of a water drop attached to a cold plate, J. Mech. Sci. Technol. 32, 2119 (2018).
  * Blake  _et al._ [2015] J. Blake, D. Thompson, D. Raps, and T. Strobl, Simulating the freezing of supercooled water droplets impacting a cooled substrate, AIAA J 53, 1725 (2015).
  * Pérez  _et al._ [2021] J. G. Pérez, S. Leclaire, S. Ammar, J. Trépanier, M. Reggio, and A. Benmeddour, Investigations of water droplet impact and freezing on a cold substrate with the lattice boltzmann method, Int. J. Thermofluids 12, 100109 (2021).
  * Wang  _et al._ [2022] C. Wang, Z. Xu, H. Zhang, J. Zheng, P. Hao, F. He, and X. Zhang, A new freezing model of sessile droplets considering ice fraction and ice distribution after recalescence, Phys. Fluids 34, 092115 (2022).
  * Tembely  _et al._ [2018] M. Tembely, R. Attarzadeh, and A. Dolatabadi, On the numerical modeling of supercooled micro-droplet impact and freezing on superhydrophobic surfaces, Int. J. Heat Mass Transf. 127, 193 (2018).
  * Ajaev and Davis [2004] V. S. Ajaev and S. H. Davis, The effect of tri-junction conditions in droplet solidification, J. Cryst. Growth 264, 452 (2004).
  * Tembely and Dolatabadi [2019] M. Tembely and A. Dolatabadi, A comprehensive model for predicting droplet freezing features on a cold substrate, J. Fluid Mech. 859, 566 (2019).
  * Schwartz and Eley [1998] L. W. Schwartz and R. R. Eley, Simulation of droplet motion on low-energy and heterogeneous surfaces, J. Colloid Interface Sci. 202, 173 (1998).
  * Moosman and Homsy [1980] S. Moosman and G. M. Homsy, Evaporating menisci of wetting fluids, J. Colloid Interface Sci. 73, 212 (1980).
  * Ajaev and Homsy [2001] V. S. Ajaev and G. M. Homsy, Steady vapor bubbles in rectangular microchannels, J. Colloid Interface Sci. 240, 259 (2001).
  * Craster  _et al._ [2009] R. V. Craster, O. K. Matar, and K. Sefiane, Pinning, retraction, and terracing of evaporating droplets containing nanoparticles, Langmuir 25, 3601 (2009).
  * Karapetsas  _et al._ [2016] G. Karapetsas, K. C. Sahu, and O. K. Matar, Evaporation of sessile droplets laden with particles and insoluble surfactants, Langmuir 32, 6871 (2016).
  * Stan  _et al._ [2023] C. A. Stan, A. Kalita, S. Marte, T. F. Kaldawi, P. R. Willmott, and S. Boutet, Rocket drops: The self-propulsion of supercooled freezing drops, Phys. Rev. Fluids 8, L021601 (2023).
  * Sebilleau  _et al._ [2021] J. Sebilleau, E. Ablonet, P. Tordjeman, and D. Legendre, Air humidity effects on water-drop icing, pre 104, L032802 (2021).
  * Charitatos and Kumar [2020] V. Charitatos and S. Kumar, A thin-film model for droplet spreading on soft solid substrates, Soft Matter 16, 8284 (2020).
  * Prosperetti and Plesset [1984] A. Prosperetti and M. S. Plesset, The stability of an evaporating liquid surface, Phys. Fluids 27, 1590 (1984).
  * Sultan  _et al._ [2005] E. Sultan, A. Boudaoud, and M. B. Amar, Evaporation of a thin film: Diffusion of the vapour and Marangoni instabilities, J. Fluid Mech. 543, 183 (2005).
  * Nagata  _et al._ [2019] Y. Nagata, T. Hama, E. H. G. Backus, M. Mezger, D. Bonn, M. Bonn, and G. Sazaki, The surface of ice under equilibrium and nonequilibrium conditions, Accounts of chemical research 52, 1006 (2019).
  * Slater and Michaelides [2019] B. Slater and A. Michaelides, Surface premelting of water ice, Nature Reviews Chemistry 3, 172 (2019).
  * Pham and Kumar [2019] T. Pham and S. Kumar, Imbibition and evaporation of droplets of colloidal suspensions on permeable substrates, Phys. Rev. Fluid 4, 034004 (2019).
  * Espín and Kumar [2015] L. Espín and S. Kumar, Droplet spreading and absorption on rough, permeable substrates, J. Fluid Mech. 784, 465 (2015).