1. Introduction
Freezing of liquid droplets on a cold surface is a fundamental multiphysics phenomenon relevant to a wide range of natural and engineering processes, including cloud microphysics (Pruppacher & Klett Reference Pruppacher and Klett1980; Hanna, Schultz & Irving Reference Hanna, Schultz and Irving2008), freezing rain formation (Zerr Reference Zerr1997; Mensah & Choi Reference Mensah and Choi2015), cryopreservation (Mazur Reference Mazur1970), food processing and storage (James, Purnell & James Reference James, Purnell and James2015; Celli, Ghanem & Su-Ling Reference Celli, Ghanem and Su-Ling2016; Zhao & Takhar Reference Zhao and Takhar2017; Zhu, Zhang & Sun Reference Zhu, Zhang and Sun2021), freeze-drying technologies (Franks Reference Franks1998; Assegehegn et al. Reference Assegehegn, Brito-de la Fuente, Franco and Gallegos2019), thermal management (Sharma et al. Reference Sharma, Parth, Shobhana, Bobin and Hardik2022), controlled fabrication of porous materials and microstructures (Qian & Zhang Reference Qian and Zhang2011; Li, Lu & Walz Reference Li, Lu and Walz2012), and icing phenomena on aircraft and civil infrastructure (Myers & Thompson Reference Myers and Thompson1998; Cao et al. Reference Cao, Wu, Su and Xu2015), to name a few.
The solidification of a droplet is driven by a complex interplay of latent heat release, thermally induced phase change, interfacial wetting and capillary forces, gravitational and viscous effects, and substrate interactions. The nonlinear interaction of these processes makes droplet freezing a rich and challenging problem in interfacial fluid mechanics. Over the past decades, extensive research has yielded a detailed understanding of the freezing behaviour of sessile, stationary droplets on flat surfaces (Jung et al. Reference Jung, Tiwari, Doan and Poulikakos2012a ,Reference Jung, Tiwari and Poulikakos b ; Wang et al. Reference Wang, Xu, Zhang, Zheng, Hao, He and Zhang2022; Kavuri et al. Reference Kavuri, Karapetsas, Sharma and Sahu2023, Reference Kavuri, Karapetsas, Sharma and Sahu2025). These studies have established that freezing proceeds in two distinct stages. The first is the recalescence stage, during which rapid ice nucleation and crystal growth occur as supercooling drives solidification, releasing latent heat and triggering a transient temperature rise (Feuillebois et al. Reference Feuillebois, Lasek, Creismeas, Pigeonneau and Szaniawski1995; Hindmarsh, Russell & Chen Reference Hindmarsh, Russell and Chen2003; Hu & Jin Reference Hu and Jin2010; Jung, Tiwari & Poulikakos Reference Jung, Tiwari and Poulikakos2012b ). This temperature rise has been attributed to the rapid restructuring and percolation of the hydrogen-bond network in water (Tavakoli, Davis & Kavehpour Reference Tavakoli, Davis and Kavehpour2015). Once the droplet reaches its melting temperature, the process transitions to a slower, typically conduction-controlled freezing stage in which the solid–liquid interface gradually advances from the substrate towards the apex of the droplet (Zadražil, Stepanek & Matar Reference Zadražil, Stepanek and Matar2006; Jung et al. Reference Jung, Tiwari and Poulikakos2012b ; Kavuri et al. Reference Kavuri, Karapetsas, Sharma and Sahu2023). Concurrently, vapour condensation on the cold substrate surrounding the droplet leads to the formation of frost halos (Jung et al. Reference Jung, Tiwari and Poulikakos2012b ; Kavuri et al. Reference Kavuri, Karapetsas, Sharma and Sahu2023).
A striking characteristic of droplet freezing is the development of a pointed ice tip during the final stage of solidification, where volume expansion arising from the density difference between liquid water and ice produces an inflection at the apex and a distinct cusp-like morphology (Sanz, Meseguer & Mayo Reference Sanz, Meseguer and Mayo1987; Anderson, Worster & Davis Reference Anderson, Worster and Davis1996; Zadražil et al. Reference Zadražil, Stepanek and Matar2006). This phenomenon, governed by interfacial thermodynamics and curvature constraints and strongly influenced by environmental and material properties such as humidity, salinity and substrate thermal properties, has led to extensive experimental research and sustained interdisciplinary interest in the freezing of sessile droplets (Anderson et al. Reference Anderson, Worster and Davis1996; Hu & Jin Reference Hu and Jin2010; Jung et al. Reference Jung, Tiwari, Doan and Poulikakos2012a , Reference Jung, Tiwari and Poulikakos2012b ; Marin et al. Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014; Zhang et al. Reference Zhang, Zhao, Lv and Yang2016; Zhang et al. Reference Zhang, Wu and Min2017, Reference Zhang, Liu, Min and Wu2019). Further studies have revealed a strong sensitivity of freezing behaviour to ambient conditions, with Sebilleau et al. (Reference Sebilleau, Ablonet, Tordjeman and Legendre2021) and Kavuri et al. (Reference Kavuri, Karapetsas, Sharma and Sahu2023, Reference Kavuri, Karapetsas, Sharma and Sahu2025) demonstrating that relative humidity and temperature play a key role in controlling frost halo formation and freezing duration. Using a two-dimensional coupled level-set and volume-of-fluid model together with an enthalpy-porosity approach, Wan et al. (Reference Wan, Liu, Chu, Wang, Wang, Wang and Sun2023) reported that increasing salinity slows the freezing process, resulting in longer solidification times for seawater droplets. Recently, Lyu et al. (Reference Lyu, Zhu, Legendre and Sun2023) demonstrated that an environmental medium with high thermal conductivity enhances external cooling, alters the freezing-front dynamics, promotes earlier crust formation and ultimately modifies the final frozen morphology of the droplet. The influence of surface roughness has also been explored (Suzuki et al. Reference Suzuki, Nakajima, Yoshida, Sakai, Hashimoto, Kameshima and Okada2007; Hao, Lv & Zhang Reference Hao, Lv and Zhang2014; Fuller, Kant & Pitchumani Reference Fuller, Kant and Pitchumani2024). Their studies showed that textured or hydrophobic surfaces introduce greater thermal resistance at the droplet–substrate interface, increase the freezing temperature and prolong the overall freezing time.
The impact and solidification of water droplets on cold substrates have also been extensively investigated. Thiévenaz, Séon & Josserand (Reference Thiévenaz, Séon and Josserand2019, Reference Thiévenaz, Josserand and Séon2020a
,
Reference Thiévenaz, Séon and Josserandb
) examined the coupled impact and freezing dynamics, elucidating contact-line pinning, dewetting and retraction behaviour, as well as the influence of interfacial ice growth on droplet spreading. Furthermore, Knight (Reference Knight1967), Demmenie et al. (Reference Demmenie, Reus, Kolpakov, Woutersen, Bonn and Shahidzadeh2023, Reference Demmenie, Gorin, Kolpakov, Smith, Kellay and Bonn2025), Thiévenaz, Josserand & Séon (Reference Thiévenaz, Josserand and Séon2020a
) and Sarlin et al. (Reference Sarlin, Papa, Grivet, Rosenbaum, Huerre, Séon and Josserand2025) demonstrated the partial wetting of water on ice, reporting a finite equilibrium contact angle of approximately
$12^\circ$
near the melting point, along with its strong temperature dependence due to contact-line pinning induced by crystallisation at lower temperatures. This is consistent with the thermodynamic interpretation of MacDowell (Reference MacDowell2026), which shows that finite contact angles can coexist with a nanometric premelting film. In addition, Huerre, Josserand & Séon (Reference Huerre, Josserand and Séon2025) provide a comprehensive review of the coupling between capillary flow and solidification, emphasising the key mechanisms and open questions in freezing dynamics.
Despite significant advances in understanding the freezing of stationary droplets on horizontal substrates, a major gap remains between laboratory studies and real-world scenarios, where droplets rarely remain stationary either before or during freezing. In many practical situations, droplets slide, deform and accelerate under gravity on inclined natural or engineered surfaces, such as aircraft wings, wind turbine blades, marine vessels, power transmission lines, heat exchangers and fluid transport pipes (Geer Reference Geer1939; Makkonen Reference Makkonen1998; Gent, Dart & Cansdale Reference Gent, Dart and Cansdale2000; Chuvilin et al. Reference Chuvilin, Naletova, Miklyaeva, Kozlova and Instanes2001; Kraj & Bibeau Reference Kraj and Bibeau2010; Cao et al. Reference Cao, Wu, Su and Xu2015). Under such conditions, droplet shapes are asymmetric, and gravity-driven motion induces persistent differences between advancing and receding sides, modifies internal pressure and curvature distributions, and alters heat transfer. While the hydrodynamics of sliding droplets have been studied extensively (Podgorski, Flesselles & Limat Reference Podgorski, Flesselles and Limat2001; Kim, Lee & Kang Reference Kim, Lee and Kang2002; Le Grand, Daerr & Limat Reference Le Grand, Daerr and Limat2005; Karapetsas, Sahu & Matar Reference Karapetsas, Sahu and Matar2013; Puthenveettil, Senthilkumar & Hopfinger Reference Puthenveettil, Senthilkumar and Hopfinger2013; Savva & Kalliadasis Reference Savva and Kalliadasis2013; Thampi, Adhikari & Govindarajan Reference Thampi, Adhikari and Govindarajan2013; Benilov & Benilov Reference Benilov and Benilov2015; Park & Kumar Reference Park and Kumar2017), virtually nothing is known about how transient motion, geometric asymmetry and dynamic contact-line behaviour influence the subsequent freezing process and the final frozen morphology. Consequently, existing literature, which overwhelmingly assumes droplets are fixed or immobilised, leaves open fundamental questions about whether the well-established morphologies and universal freezing trajectories observed in sessile droplets persist when motion precedes or accompanies freezing, highlighting the critical role of surface inclination in influencing freezing dynamics.
The shape of a frozen droplet, including the asymmetry of the cusp-like morphology at its tip, is strongly influenced by its initial equilibrium shape (Jin et al. Reference Jin, Wang, Sui and Yang2016; Zeng et al. Reference Zeng, Lyu, Legendre and Sun2022), highlighting the importance of understanding the initial droplet morphology, particularly on inclined or more complex non-horizontal substrates. Numerous studies have investigated the dynamics of droplets sliding on inclined substrates, focusing on contact line behaviour, velocity and shape evolution, albeit in the absence of freezing and evaporation (Podgorski et al. Reference Podgorski, Flesselles and Limat2001; Kim et al. Reference Kim, Lee and Kang2002; Thiele et al. Reference Thiele, Neuffer, Bestehorn, Pomeau and Velarde2002; Le Grand et al. Reference Le Grand, Daerr and Limat2005; Karapetsas et al. Reference Karapetsas, Sahu and Matar2013; Puthenveettil et al. Reference Puthenveettil, Senthilkumar and Hopfinger2013; Savva & Kalliadasis Reference Savva and Kalliadasis2013; Thampi et al. Reference Thampi, Adhikari and Govindarajan2013; Benilov & Benilov Reference Benilov and Benilov2015; Park & Kumar Reference Park and Kumar2017). Podgorski et al. (Reference Podgorski, Flesselles and Limat2001) demonstrated that droplets of silicone oil and water on fluoro-polymer-coated surfaces exhibit rounded, cornered or pearling shapes depending on the balance between viscous and surface tension forces, which also governs the dynamic contact angle at the trailing edge. On rough surfaces, droplets begin to slide when the substrate is tilted beyond a critical angle. Kim et al. (Reference Kim, Lee and Kang2002) showed that ethylene glycol, glycerine and glycerine–water mixture droplets on polycarbonate surfaces exhibit sliding velocities consistent with a linear balance between gravitational and viscous–capillary forces. Le Grand et al. (Reference Le Grand, Daerr and Limat2005) reported that, with increasing Bond number (the ratio between the gravitational and surface tension forces), a millimetre-sized drop on a partially wetting inclined plane elongates and undergoes morphological transitions from rounded to cornered, cusped and ultimately pearling shapes. In addition, Puthenveettil et al. (Reference Puthenveettil, Senthilkumar and Hopfinger2013) found that at high Reynolds numbers, droplet velocity scales with the Bond number through a balance of inertial, viscous and capillary forces, leading mercury droplets to form cornered shapes and water droplets to develop rivulets.
In parallel with these experimental investigations, several theoretical and numerical studies (Thiele et al. Reference Thiele, Neuffer, Bestehorn, Pomeau and Velarde2002; Karapetsas et al. Reference Karapetsas, Sahu and Matar2013; Savva & Kalliadasis Reference Savva and Kalliadasis2013; Thampi et al. Reference Thampi, Adhikari and Govindarajan2013; Benilov & Benilov Reference Benilov and Benilov2015; Park & Kumar Reference Park and Kumar2017) have examined droplet dynamics on inclined substrates using long-wave, lubrication and diffuse-interface models. Thiele et al. (Reference Thiele, Neuffer, Bestehorn, Pomeau and Velarde2002) identified universal and non-universal regimes of sliding drops governed by substrate inclination, while Karapetsas et al. (Reference Karapetsas, Sahu and Matar2013) showed that thermo-capillarity and temperature-dependent wettability significantly influence contact line motion and spreading. Thampi et al. (Reference Thampi, Adhikari and Govindarajan2013) distinguished rolling and sliding modes based on droplet geometry, viscosity contrast and slip length. Benilov & Benilov (Reference Benilov and Benilov2015) derived asymptotic expressions for thin drop velocities using a Navier-slip lubrication framework and validated their model with experiments. Savva & Kalliadasis (Reference Savva and Kalliadasis2013) investigated two-dimensional droplets on chemically heterogeneous substrates using a long-wave evolution equation with slip, revealing stick–slip dynamics, substrate-induced hysteresis and uphill motion under strong chemical gradients. Park & Kumar (Reference Park and Kumar2017) studied droplets on substrates with topographical defects using a thin-film evolution equation with a precursor layer, illustrating that pinning–depinning and residual droplet formation are controlled by interactions of the receding contact line with defect geometry.
The freezing of droplets on inclined surfaces has also attracted increasing attention in recent years due to its relevance to several of the applications discussed earlier. Starostin et al. (Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2022, Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2023) experimentally examined freezing on horizontal, inclined, inverted and silicone-oil-lubricated surfaces, showing that the tip angles of frozen droplets remain nearly constant at
$143^\circ \pm 11^\circ$
, and analysed how asymmetric cooling governs tip orientation. Dang et al. (Reference Dang, Song, Zhang, Pekař and Hosseini2022) modelled the freezing of sessile droplets on inclined substrates and demonstrated that smaller droplets are more sensitive to inclination, whereas larger droplets exhibit only modest changes, with freezing occurring more readily on steeper surfaces, an effect that is particularly pronounced for small droplets. More recently, Kumar, Subudhi & Das (Reference Kumar, Subudhi and Das2025) conducted an experimental study on the sequential freezing and melting of impact-deposited and sessile droplets, investigating contact-line arrest, tip singularity, melting modes and interfacial dynamics over a broad range of temperatures, impact velocities and substrate inclination angles. Furthermore, in the studies of Dang et al. (Reference Dang, Song, Zhang, Pekař and Hosseini2022), Starostin et al. (Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2022, Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2023) and Kumar et al. (Reference Kumar, Subudhi and Das2025), the contact line becomes arrested after the initial spreading or deformation stage, before the onset of solidification. Consequently, the freezing process proceeds from an essentially stationary droplet configuration.
In the present work, we theoretically investigate the freezing dynamics of a droplet sliding on a smooth substrate, where the droplet continues to move until contact-line arrest occurs due to solidification. Although droplets typically do not slide on smooth surfaces unless the inclination exceeds a critical angle or the droplet is sufficiently large for gravity to overcome capillary adhesion (Quéré, Azzopardi & Delattre Reference Quéré, Azzopardi and Delattre1998), we assume that contact-angle hysteresis is sufficiently small to permit sliding over the range of inclinations considered. Under these conditions, asymmetric freezing can substantially alter post-solidification contact angles, shift the cusp singularity away from the surface normal, and is strongly influenced by surface wettability, substrate inclination, freezing rate and the interplay between gravitational and capillary forces. To address these unresolved issues, we develop a unified theoretical–numerical framework that couples lubrication-based hydrodynamics, dynamic contact-line motion and conduction-controlled solidification. In contrast to earlier studies that assume a frozen geometry or prescribe the motion of the freezing front, our model captures the fully self-consistent evolution of the liquid shape, freezing front, contact angles and wetting state from the onset of solidification to complete freezing. By systematically varying substrate inclination and surface wettability, we uncover a rich hierarchy of freezing behaviours that are not observed in stationary droplets. Our results reveal that substrate inclination induces strong and persistent asymmetry in both droplet morphology and freezing dynamics. During the early stages, the advancing and receding contact angles diverge, resulting in an asymmetric liquid layer above the freezing front and a shift of the final cusp relative to the contact line normal, with the degree of asymmetry increasing monotonically with inclination angle of the substrate. Remarkably, even brief droplet motion prior to contact-line arrest leaves a lasting imprint on the frozen structure. As freezing progresses, a dynamical cross-over emerges in which gravity dominates the early behaviour, while the diminishing liquid volume enhances capillary effects that gradually restore symmetry and drive the system towards a universal tip angle. These findings demonstrate that the freezing of moving droplets is not a trivial extension of freezing of sessile droplets, but a qualitatively distinct regime in which transient hydrodynamics, geometry and phase change remain strongly coupled, providing the first mechanistic insight into asymmetric freezing governed by substrate tilt, surface wettability and motion-induced geometric memory.
The remainder of this paper is structured as follows. In § 2, we present the problem formulation, governing equations, scaling considerations and the numerical approach used within the framework of the lubrication approximation. In § 3, we discuss the results, focusing on the mechanism governing the freezing of a sliding drop, supported by an extensive parametric study. Finally, concluding remarks are provided in § 4.
2. Problem formulation
We examine the freezing dynamics of a sessile droplet on an inclined substrate using numerical simulations based on the lubrication approximation. A schematic of the two-dimensional (2-D) droplet at a given time,
$t$
, during the freezing process, illustrating the physical parameters used in the model, is shown in figure 1. A Cartesian coordinate system
$(x,z)$
is employed in the theoretical modelling, whose origin is located at the centre of the droplet on the substrate. The substrate, inclined at an angle
$\alpha$
, has a thickness
$d_w$
, thermal conductivity
$\lambda _w$
and specific heat capacity
$C_{pw}$
, with its bottom surface maintained at a constant temperature
$T_c$
. Initially (at
$t = 0$
), when freezing begins, the droplet on the inclined substrate is considered to have reached its equilibrium shape and steady sliding velocity, with the initial recalescence phase completed and the liquid at the melting temperature (
$T_m$
). Since our focus is on a droplet that is already sliding, this assumption ensures that the droplet is in a steady state before freezing starts. Prior to the propagation of the freezing front, a very rapid recalescence phase occurs (Feuillebois et al. Reference Feuillebois, Lasek, Creismeas, Pigeonneau and Szaniawski1995; Hindmarsh et al. Reference Hindmarsh, Russell and Chen2003; Hu & Jin Reference Hu and Jin2010; Jung et al. Reference Jung, Tiwari and Poulikakos2012b
; Wang et al. Reference Wang, Xu, Zhang, Zheng, Hao, He and Zhang2022), lasting only a few milliseconds. During this phase, the release of latent heat forms an ice crystal scaffold (or mushy phase) and induces a transient rise in temperature. As the recalescence time scale is much shorter than the total freezing time, it is reasonable to assume that this phase has already concluded and that the liquid is at
$T_m$
when the freezing front begins to propagate. This initial recalescence explains why the unfrozen liquid does not exhibit a contact angle of
$12^\circ$
at the liquid–ice interface immediately after freezing onset, as noted in previous studies (Knight Reference Knight1967; Thiévenaz et al. Reference Thiévenaz, Josserand and Séon2020a
; Demmenie et al. Reference Demmenie, Reus, Kolpakov, Woutersen, Bonn and Shahidzadeh2023, Reference Demmenie, Gorin, Kolpakov, Smith, Kellay and Bonn2025; Huerre et al. Reference Huerre, Josserand and Séon2025; Sarlin et al. Reference Sarlin, Papa, Grivet, Rosenbaum, Huerre, Séon and Josserand2025).
Schematic of a sessile droplet freezing on a solid substrate inclined at an angle
$\alpha$
. The droplet length and maximum height are denoted by
$L$
and
$H$
, respectively. The quantities
$S_\infty$
and
$H_\infty$
represent the initial thickness of the microscopic ice layer and the precursor film. The advancing and receding contact angles of the droplet are
$\theta _a$
and
$\theta _r$
. The advancing and receding contact angles at the liquid–ice interface are
$\theta _{\textit{ice},a}$
and
$\theta _{\textit{ice},r}$
. The left and right contact lines of the unfrozen droplet are denoted by
$x_{cl}$
and
$x_{cr}$
. The ambient and substrate temperatures are
$T_g$
and
$T_c$
, respectively. The cusp angle of the droplet, measured with respect to the
$z$
-axis, is denoted by
$\theta _{\textit{Cusp}}$
.
At
$t = 0$
, a uniformly thin ice layer of thickness
$S_\infty$
has already formed along the substrate, indicating the onset of the slower solidification stage governed by heat conduction. The droplet length, maximum height, initial ice layer thickness, precursor film thickness, and the advancing and receding contact angles at the droplet–substrate and liquid–ice interfaces are denoted by
$L$
,
$H$
,
$S_\infty$
,
$H_\infty$
,
$\theta _a$
,
$\theta _r$
,
$\theta _{\textit{ice},a}$
and
$\theta _{\textit{ice},r}$
, respectively. The angle formed by the singularity peak of the droplet with the vertical (
$z$
-axis) at the droplet centre is denoted by
$\theta _{{Cusp}}$
. The liquid is treated as an incompressible, Newtonian fluid with constant density
$\rho _l$
, dynamic viscosity
$\mu _l$
, specific heat capacity
$C_{pl}$
and thermal conductivity
$\lambda _l$
. The liquid–ice and liquid–gas interfaces are represented by
$z = s(x,t)$
and
$z = h(x,t)$
, respectively. The surface tension at the liquid–gas interface is denoted by
$\gamma _{lg}$
. The density
$\rho _s$
, specific heat capacity
$C_{ps}$
and thermal conductivity
$\lambda _s$
of the solid (ice) phase are also assumed to be constant.
Within the lubrication approximation, the droplet is very thin, with its maximum height
$H$
much smaller than its length
$L$
(
$H \ll L$
). Although this approximation strictly applies only to thin droplets, previous studies show that lubrication-based models provide reasonably accurate predictions for contact angles up to
$60^\circ$
(Tembely & Dolatabadi Reference Tembely and Dolatabadi2019; Charitatos & Kumar Reference Charitatos and Kumar2020). Consequently, the small aspect ratio
$\epsilon = 2H/L \ll 1$
justifies using lubrication theory to derive the evolution equations describing the freezing dynamics of the sessile droplet, including the liquid, ice and gas phases. In the present study, we choose the equilibrium contact angle range of
$12^\circ$
–
$30^\circ$
to ensure that this thin-droplet assumption remains valid. The classical no-slip condition at a moving contact line leads to a stress singularity, which can be alleviated by introducing a thin precursor film. This approach has been widely used for droplet motion (Schwartz & Eley Reference Schwartz and Eley1998; Espín & Kumar Reference Espín and Kumar2015; Park & Kumar Reference Park and Kumar2017), droplet evaporation and imbibition (Karapetsas, Sahu & Matar Reference Karapetsas, Sahu and Matar2016; Pham & Kumar Reference Pham and Kumar2019; Wang et al. Reference Wang, Karapetsas, Valluri, Sefiane, Williams and Takata2021; Williams et al. Reference Williams, Karapetsas, Mamalis, Sefiane, Matar and Valluri2021), and droplet freezing (Zadražil et al. Reference Zadražil, Stepanek and Matar2006; Tembely & Dolatabadi Reference Tembely and Dolatabadi2019; Kavuri et al. Reference Kavuri, Karapetsas, Sharma and Sahu2023, Reference Kavuri, Karapetsas, Sharma and Sahu2025). This approach successfully captures key features such as cusp formation, freezing front propagation and volume expansion, while avoiding numerical difficulties at the contact line. The dimensional governing equations describing the freezing process of a sessile droplet on an inclined substrate are presented in the following.
2.1. Dimensional governing equations
2.1.1. Liquid phase
The liquid-phase dynamics is governed by the conservation equations for mass, momentum and energy, which are expressed as
\begin{align} \rho _l \left ( \frac {\partial {\boldsymbol{v}}}{\partial t} + {\boldsymbol{v}} \boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol{v}} \right ) &= - \boldsymbol{\nabla }p + \mu _l {\nabla} ^2 {\boldsymbol{v}} + \rho _l g \,{\mathrm{sin}} \alpha \, e_x-\rho _l g \, {\mathrm{cos}}\alpha \, {e_z}, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}} &= 0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho _l C_{pl} \left ( \frac {\partial T_l}{\partial t} + {\boldsymbol{v}} \boldsymbol{\cdot }\boldsymbol{\nabla }T_l \right ) &= \lambda _l {\nabla} ^2 T_l, \\[-12pt] \nonumber \end{align}
where
$\boldsymbol{v}$
,
$p$
and
$T_l$
denote the velocity, pressure and temperature in the liquid phase, respectively;
$\boldsymbol{\nabla}$
is the gradient operator;
$g$
represents the acceleration due to gravity; and
$e_x$
and
$e_z$
denote the components of the unit vectors along the
$x$
- and
$z$
-directions, respectively.
At the free surface (
$z=h(x,t)$
), the outward unit normal to the liquid–gas interface (
${\boldsymbol{n}}_l$
) is defined as
\begin{equation} {\boldsymbol{n}}_l = (-h_{x}e_x + e_z)/\sqrt {\left(h_{x}^2+1\right)}. \end{equation}
The tangential components of the velocity are identical in both phases and are given by
${\boldsymbol{v}}_{\tau } = {\boldsymbol{v}} - ({\boldsymbol{v}} \boldsymbol{\cdot }{\boldsymbol{n}}_l) {\boldsymbol{n}}_l = {\boldsymbol{v}}_{lg} - ({\boldsymbol{v}}_{lg} \boldsymbol{\cdot }{\boldsymbol{n}}_l) {\boldsymbol{n}}_l$
. At
$z = h(x,t)$
, the velocity field satisfies the local mass, force, and energy balance between the liquid and gas phases (Karapetsas et al. Reference Karapetsas, Sahu and Matar2016), which are given by
\begin{align} \rho _l ({\boldsymbol{v}} -{\boldsymbol{v}}_{lg}) \boldsymbol{\cdot }{\boldsymbol{n}}_l &= \rho _g ({\boldsymbol{v}}_g -{\boldsymbol{v}}_{lg}) \boldsymbol{\cdot }{\boldsymbol{n}}_l, \\[-28pt] \nonumber \end{align}
\begin{align} {\boldsymbol{n}}_l \boldsymbol{\cdot }\left [-p {\boldsymbol{I}} + \mu _l \left ( \boldsymbol{\nabla }{\boldsymbol{v}} + \boldsymbol{\nabla }{\boldsymbol{v}} ^T \right ) \right ] &= -p_g {\boldsymbol{n}}_l + \gamma _{lg} \kappa _{lg} {\boldsymbol{n}}_l + \varPi {\boldsymbol{n}}_l , \\[-28pt] \nonumber \end{align}
\begin{align} \lambda _l \boldsymbol{\nabla }T_l \boldsymbol{\cdot }{\boldsymbol{n}}_l - \lambda _g \boldsymbol{\nabla }T_g \boldsymbol{\cdot }{\boldsymbol{n}}_l -h_{tc}(T_{g}-T_{lg})&= 0, \\[-12pt] \nonumber \end{align}
where
$\rho _g$
,
$\lambda _g$
,
${\boldsymbol{v}}_g$
,
$T_g$
, denote the density, thermal conductivity, velocity field and temperature of the gas phase, respectively;
$v_{lg}$
denotes the velocity of the liquid–gas interface, and
$h_{tc}$
denotes the convective heat transfer coefficient. Here,
$\boldsymbol{I}$
represents the identity tensor,
$T_{lg}$
denotes the temperature at the liquid–gas interface,
$\kappa _{lg} = - \boldsymbol{\nabla} _{s,l} \boldsymbol{\cdot }{\boldsymbol{n}}_l$
is the mean curvature of the free surface and
$\boldsymbol{\nabla} _{s,l} = (I - {\boldsymbol{n}}_l {\boldsymbol{n}}_l) \boldsymbol{\cdot }\boldsymbol{\nabla}$
is the surface gradient operator. The disjoining pressure (
$\varPi$
), which accounts for the van der Waals interactions, is given by
\begin{equation} \varPi = A \left [ \left ( \frac {B}{h - s} \right )^{n} - \left ( \frac {B}{h - s} \right )^{m} \right ], \end{equation}
where
$A = A_{{Ham}} / B^{n} \geqslant 0$
quantifies the magnitude of the intermolecular interaction energy between the liquid–gas and liquid–ice interfaces, and
$B$
is the precursor layer thickness. Here,
$n \gt m \gt 1$
and
$A_{{Ham}}$
represents the dimensional Hamaker constant.
At the liquid–ice interface
$(z = s(x,t))$
, the velocity field satisfies the kinematic boundary condition:
\begin{equation} {\boldsymbol{v}} = {\boldsymbol{v}}_{ls} - \left ( \frac {J_s}{\rho _l} \right ) {\boldsymbol{n}}_s, \end{equation}
where
$J_s$
denotes the freezing mass flux and
${\boldsymbol{v}}_{ls}$
is the velocity of the liquid–ice interface. The outward unit normal to the liquid–ice interface (
${\boldsymbol{n}}_s$
) is given by
\begin{equation} {\boldsymbol{n}}_s = {\left ( {-s_x e_x + e_z} \right ) / {\sqrt {s_x^2 + 1}}}. \end{equation}
A no-slip condition is imposed at the liquid–ice interface, requiring the tangential component of the velocity to satisfy
\begin{equation} {\boldsymbol{v}} \boldsymbol{\cdot }{\boldsymbol{t}}_s = 0, \end{equation}
where
${\boldsymbol{t}}_s$
is the outward unit tangential vector along the liquid–ice interface, defined as
\begin{equation} {\boldsymbol{t}}_s = {\left (e_x + s_x e_z \right ) / \sqrt {s_x^2 + 1}}. \end{equation}
2.1.2. Solid (ice) phase
The temperature evolution in the solid (ice) phase is governed by
\begin{equation} \rho _s C_{ps} \frac {\partial T_s}{\partial t} = \lambda _s {\nabla} ^2 T_s, \end{equation}
where
$T_s$
represents the temperature in the solid phase.
At the solid substrate (
$z=0$
), the continuity of temperature is imposed, which is given by
\begin{equation} T_s = T_w, \end{equation}
where
$T_w$
denotes the substrate temperature at
$z=0$
.
At the freezing front (
$z=s(x,t)$
), the temperature satisfies
\begin{equation} T_s = T_l = T_f, \end{equation}
where the equilibrium temperature at the freezing front,
$T_f$
, is assumed to be equal to the melting temperature,
$T_m$
.
The conservation of mass and energy at the liquid–solid interface (at
$z = s(x,t)$
) yields
\begin{align} J_s = \rho _l ({\boldsymbol{v}}_{ls} - {\boldsymbol{v}})\boldsymbol{\cdot }{\boldsymbol{n}}_s &= \rho _s ({\boldsymbol{v}}_{ls} - {\boldsymbol{v}}_s)\boldsymbol{\cdot }{\boldsymbol{n}}_s, \end{align}
\begin{align} \rho _s ({\boldsymbol{v}}_{ls} - {\boldsymbol{v}}_s)\boldsymbol{\cdot }{\boldsymbol{n}}_s H_s - \lambda _s \boldsymbol{\nabla }T_s \boldsymbol{\cdot }{\boldsymbol{n}}_s &= \rho _l ({\boldsymbol{v}}_{ls} - {\boldsymbol{v}})\boldsymbol{\cdot }{\boldsymbol{n}}_s H_l - \lambda _l \boldsymbol{\nabla }T_l \boldsymbol{\cdot }{\boldsymbol{n}}_s, \end{align}
where
${\boldsymbol{v}}_s$
is the velocity of the ice phase, which is stationary in our study (
${\boldsymbol{v}}_s = 0$
);
$H_s$
and
$H_l$
denote the enthalpy of the solid and liquid phases, respectively.
Combining (2.16) and (2.17), and employing
${\boldsymbol{v}}_s = 0$
, we obtain
\begin{equation} J_s \Delta H_{sl} - \lambda _s \boldsymbol{\nabla }T_s \boldsymbol{\cdot }{\boldsymbol{n}}_s + \lambda _l \boldsymbol{\nabla }T_l \boldsymbol{\cdot }{\boldsymbol{n}}_s = 0, \end{equation}
where
$\Delta H_{sl} = H_s - H_l$
represents the enthalpy jump across the liquid–ice interface. Then, substituting
$H_s = C_{ps}(T_f - T_m) + L_f(T_m)$
and
$H_l = C_{pl}(T_f - T_m)$
, we get
\begin{equation} \Delta H_{sl} = (C_{ps} - C_{pl})(T_f - T_m) + L_f, \end{equation}
where
$L_f$
is the latent heat of fusion defined with
$T_m$
as the reference temperature. In the following, (2.18) is used to determine the evolution of the freezing front,
$s(x,t)$
.
2.1.3. Solid substrate
The energy equation governing heat transfer in the solid substrate is
\begin{equation} \rho _w C_{pw} \frac {\partial T_w}{\partial t} = \lambda _w {\nabla} ^2 T_w, \end{equation}
where
$\rho _w$
denotes the density of the substrate. This equation is supplemented by the continuity of heat flux at the ice–substrate interface (
$z = 0$
):
\begin{equation} \lambda _s \frac {\partial T_s}{\partial z} = \lambda _w \frac {\partial T_w}{\partial z}, \end{equation}
and the temperature condition at the bottom boundary of the substrate (
$z = -d_w$
):
\begin{equation} T_w = T_c. \end{equation}
The fluid properties and ranges of physical parameters considered in this study are summarised in table 1. All simulations are performed for water droplets on a highly conductive copper substrate, with the bottom surface maintained at approximately
$-30\,^\circ\mathrm{C}$
and the ambient temperature close to the melting point of water. Under these conditions, evaporation is minimal and is therefore neglected in the present model. In the following, we non-dimensionalise the governing equations and boundary conditions.
Representative values of the physical parameters used in our simulations. The listed properties correspond to the water–air system and copper substrate.
2.2. Scaling
The above-mentioned governing equations and boundary conditions are non-dimensionalised using the following scalings, where a tilde denotes the dimensionless variable:
\begin{align}&\qquad\quad ( x, z, h, D_{w} ) = l ( \tilde {x}, \epsilon \tilde {z}, \epsilon \tilde {h} , \epsilon \tilde {D_{w}}), \, t = \frac {l}{U} \tilde {t}, \, ( u, w ) = U ( \tilde {u}, \epsilon \tilde {w} ),\nonumber\\& (u_g, w_g ) = U ( \tilde {u_{g}}, \epsilon \tilde {w_{g}} ), \, (p,\varPi ) = \frac {\mu _l U l}{H^2} (\tilde {p},\tilde {\varPi }), \, T_i = \Delta T \; \tilde {T}_i + T_c \; (i=l,s,w), \nonumber\\& \qquad\qquad\qquad J_s = \epsilon \rho _s U \; \tilde {J}_s, \boldsymbol{\nabla }= \frac {1}{l} \tilde {\boldsymbol{\nabla }}, \boldsymbol{\nabla} _{s,i} = \frac {1}{l} \tilde {\boldsymbol{\nabla} }_{s,i} \; (i=l,s), \end{align}
where
$\Delta T = T_m - T_c$
,
$l = L/2$
(half-length of the droplet),
$\epsilon = H/l$
,
$\widetilde {\boldsymbol{\nabla }} = {\boldsymbol{e}}_x \widetilde {\partial }_x + {\boldsymbol{e}}_z \epsilon ^{-1} \widetilde {\partial }z$
and
$\widetilde {\boldsymbol{\nabla }}{s,i} = (\boldsymbol{I} - {\boldsymbol{n}} {\boldsymbol{n}}) \boldsymbol{\cdot }\widetilde {\boldsymbol{\nabla }}$
$(i = l, s)$
. The velocity scale is chosen as
$U = \epsilon ^3 \gamma _{lg} / \mu _l$
, such that
$Ca/\epsilon ^2 = 1$
, wherein
$Ca = \mu _l U / (\epsilon \gamma _{lg})$
is the capillary number. Hereafter, the tilde notation is suppressed, and the partial derivatives
$\partial / \partial x$
,
$\partial / \partial z$
and
$\partial / \partial t$
are represented by the subscripts
$x$
,
$z$
and
$t$
, respectively. By applying these scalings and invoking the lubrication approximation (
$\epsilon \ll 1$
), we obtain the dimensionless governing equations and boundary conditions for the liquid, ice and gas phases.
2.2.1. Liquid phase
The governing equations for the liquid phase, written in dimensionless form, are as follows:
\begin{align} \partial ^2_z u + \frac {Bo\,{\mathrm{sin}}\alpha }{\epsilon } &= \partial _x p, \\[-32pt] \nonumber \end{align}
\begin{align} \partial _z p + Bo \, {\mathrm{cos}}\alpha &= 0, \\[-32pt] \nonumber \end{align}
\begin{align} \partial _x u + \partial _z w &= 0, \\[-32pt] \nonumber \end{align}
\begin{align} \partial ^2_z T_l &= 0. \\[0pt] \nonumber \end{align}
The boundary conditions at the liquid–gas interface (
$z = h(x,t)$
) are
\begin{align} p &= - \kappa _{lg} - \varPi , \\[-32pt] \nonumber \end{align}
\begin{align} \partial _z u &= 0, \\[-32pt] \nonumber \end{align}
\begin{align} \partial _z T_l &= Bi(T_v - T_l), \\[-32pt] \nonumber \end{align}
\begin{align} \partial _t h + u\partial _x h - w &= 0. \\[0pt] \nonumber \end{align}
Here,
$ \textit{Bo} = \rho _l g l^{2}/\gamma _{lg}$
is the Bond number, representing the ratio of gravitational to surface tension forces, and
$Bi = h_{tc} H/\lambda _l$
is the Biot number, which compares the internal thermal resistance due to conduction to the external resistance associated with convective heat transfer.
The boundary conditions at the liquid–ice interface (at
$z = s(x,t)$
) are
\begin{equation} u = 0, \quad T_l = T_f, \end{equation}
where
$T_f = T_m = 1$
.
The dimensionless disjoining pressure can be expressed as
\begin{equation} \varPi = A_{n} \epsilon ^{-2} \left [ \left ( \frac {\beta }{h-s} \right )^{n} - \left ( \frac {\beta }{h-s} \right )^{m} \right ], \end{equation}
where
$\beta$
is of the same order as the equilibrium precursor film thickness (Pham & Kumar Reference Pham and Kumar2019) and
$A_{n} = H A / \gamma _{lg}$
is the dimensionless Hamaker constant. The competition between the attractive and repulsive contributions in (2.33) determines the equilibrium contact angle,
$\theta _{\textit{eq}}$
(Schwartz & Eley Reference Schwartz and Eley1998; Zadražil et al. Reference Zadražil, Stepanek and Matar2006; Espín & Kumar Reference Espín and Kumar2015; Pham & Kumar Reference Pham and Kumar2019; Tembely & Dolatabadi Reference Tembely and Dolatabadi2019), which may be approximated by (Pham & Kumar Reference Pham and Kumar2019)
\begin{equation} \theta _{\textit{eq}} \approx \sqrt {\beta A_{n}}. \end{equation}
The mean curvatures of the liquid–gas and liquid–ice interfaces, retaining higher-order terms (
$\epsilon ^{2} h_x^{2}$
and
$\epsilon ^{2} s_x^{2}$
), are given by
\begin{equation} \kappa _{lg} = \frac {h_{\textit{xx}}}{\left(1+\epsilon ^{2} h_x^{2}\right)^{3/2}} \quad \textrm {and} \quad \kappa _{sl} = \frac {s_{\textit{xx}}}{\left(1+\epsilon ^{2} s_x^{2}\right)^{3/2}}, \end{equation}
respectively.
2.2.2. Solid (ice) phase
The dimensionless energy conservation equation in the ice phase is given by
\begin{equation} \partial ^2_z T_s = 0. \end{equation}
The corresponding energy balance at the liquid–ice interface (
$z = s(x,t)$
) is
\begin{equation} Ste \left ( \varLambda _s \partial _z T_s - \partial _z T_l \right ) = J_s, \end{equation}
where
$\varLambda _s = \lambda _s / \lambda _l$
is the thermal conductivity ratio and
${\textit{Ste}} = \lambda _l \Delta T / (\epsilon \rho _s U L_f H)$
is the Stefan number, which represents the ratio of sensible (conductive) heat transport in the liquid to the latent heat released at the solidification front.
The non-dimesionalisation of the mass conservation at
$z = s(x,t)$
(2.9) yields
\begin{equation} \partial _t s - w = D_s J_s, \end{equation}
where
$J_s = \partial _t s$
,
$D_s = \rho _s / \rho _l$
denotes the density ratio of the solid to the liquid phase.
At the ice–substrate interface (
$z = 0$
), the temperature satisfies
\begin{equation} T_s = T_w. \end{equation}
2.2.3. Solid substrate
The dimensionless energy equation for the solid substrate is given by
\begin{equation} \partial ^2_z T_w = 0. \end{equation}
At the solid–ice interface (
$z = 0$
), the temperature field satisfies
\begin{equation} \varLambda _s \partial _z T_s = \varLambda _w \partial _z T_w, \end{equation}
where
$\varLambda _w = \lambda _w / \lambda _l$
is the thermal conductivity ratio for the substrate.
At the bottom surface of the substrate
$z = -D_w$
, where
$D_w = d_w/H$
is the scaled wall thickness, the temperature is prescribed as
\begin{equation} T_w = 0. \end{equation}
2.3. Evolution equations
Integrating (2.24) and (2.25) with respect to
$z$
, and applying the boundary conditions in (2.29), (2.32) and (2.28), yields
\begin{align} u &= \left ( \partial _x p - \frac {Bo \, {\mathrm{sin}}\alpha }{\epsilon } \right ) \left ( \frac {z^2 - s^2}{2} - h(z - s) \right ), \end{align}
\begin{align} p &= Bo (h - z){\mathrm{cos}}\alpha - \kappa _{lg} - \varPi . \end{align}
Integrating the continuity equation (2.26) and applying the kinematic condition (2.31) along with the mass conservation equation (2.38) yields the evolution equation
\begin{equation} \partial _t h - \partial _t s = -\partial _x q_l - D_s J_s, \end{equation}
where the volumetric flux is given by
\begin{equation} q_l = \left ( \partial _x p - \frac {Bo \, {\mathrm{sin}}\alpha }{\epsilon } \right ) \left ( -\frac {h^3}{3} + s h^2 - s^2 h + \frac {s^3}{3} \right ). \end{equation}
Similarly, integrating (2.27) and applying the thermal boundary conditions ((2.30) and (2.32)) gives
\begin{equation} T_l = \left [ Bi \left ( T_v - T_l\big |_{h} \right ) \right ] (z - s) + T_f. \end{equation}
The temperature field in the ice phase is given by
\begin{equation} T_s = \frac {T_f}{D_w + s \varLambda _w/\varLambda _s} \left ( D_w + z \frac {\varLambda _w}{\varLambda _s} \right ). \end{equation}
Substituting this expression, along with (2.38), in (2.37) yields
\begin{equation} \partial _t s = Ste \left ( \frac {\varLambda _w T_f}{D_w+s \varLambda _w/\varLambda _s} - Bi\left (T_{v}-T_{l}\big |_{h}\right )\right ). \end{equation}
Finally, the temperature profile in the solid substrate is given by
\begin{equation} T_w = \frac {T_f}{D_w + s \varLambda _w/\varLambda _s} (z + D_w). \end{equation}
2.4. Initial and boundary conditions
In the theoretical model incorporated in the present study based on lubrication approximation, the droplet is placed on a thin precursor film that lies above a thin ice layer. By suitably selecting the dimensionless Hamaker constant
$A_{n}$
((2.33) and (2.34)), we ensure that the droplet attains the prescribed equilibrium contact angle with the substrate. This choice effectively balances the attractive and repulsive contributions of the disjoining pressure, thereby setting the equilibrium contact angle (Schwartz & Eley Reference Schwartz and Eley1998; Zadražil et al. Reference Zadražil, Stepanek and Matar2006; Espín & Kumar Reference Espín and Kumar2015; Pham & Kumar Reference Pham and Kumar2019; Tembely & Dolatabadi Reference Tembely and Dolatabadi2019). The initial conditions imposed on the domain are
\begin{eqnarray} h(x,t=0) = \max \big (h_{\infty }+s_{\infty }-x^{2}+2xx_{0}+1-x_{0}^{2}, \, h_{\infty }+s_{\infty }\big ). \end{eqnarray}
The precursor layer thickness far from the droplet is defined as
$h_{\infty } = (h - s_{\infty }) = H_{\infty }/H$
, and is assumed to be uniform, flat and of zero mean curvature. The parameter
$x_{0}$
denotes the initial position of the centre of the droplet. In all simulations,
$\beta = 0.01$
and the initial thickness of the ice layer is taken as
$s_{\infty } = 10^{-3}$
. We confirm that reducing the value of the latter parameter by one or two orders of magnitude does not affect the results.
To prevent unphysical solidification within the precursor film, we introduce a thickness-dependent term that suppresses freezing in this region. This is essential, as solidification of the precursor layer would reintroduce the stress singularity, defeating its purpose. Following Zadražil et al. (Reference Zadražil, Stepanek and Matar2006), Tembely & Dolatabadi (Reference Tembely and Dolatabadi2019) and Kavuri et al. (Reference Kavuri, Karapetsas, Sharma and Sahu2023, Reference Kavuri, Karapetsas, Sharma and Sahu2025), we incorporate a thickness-dependent Stefan number,
${\textit{Ste}}(x)$
, which inhibits freezing within the precursor film and is defined as
\begin{equation} Ste(x) = \frac {1}{2}\big (1 + \tanh \big [4 \times 10^{3}\big ((h - s) - 1.4\beta \big )\big ]\big ) Ste. \end{equation}
The boundary conditions applied at the edges of the computational domain are
\begin{align} h_{x}(0,t) &= h_{xxx}(0,t) = h_{x}(x_{\pm \infty },t) = h_{xxx}(x_{\pm \infty },t) = 0, \end{align}
\begin{align} s_{x}(0,t) &= s_{xxx}(0,t) = s_{x}(x_{\pm \infty },t) = s_{xxx}(x_{\pm \infty },t) = 0, \end{align}
where
$\pm x_{\infty }$
denotes the far-field boundary of the domain.
2.5. Numerical method and validation
The dimensionless governing equations are discretised using the Galerkin finite element method, with weak forms formulated for each equation. The nonlinear system is solved iteratively using the Newton–Raphson method and time integration is performed with an implicit Euler scheme featuring adaptive time stepping. The adaptive time step is controlled by the maximum residual from the preceding iteration, consistent with standard adaptive implicit Euler procedures. Linear algebra operations are carried out using the LAPACK library and the complete solver is implemented in FORTRAN. Our model is validated for both droplet sliding on a smooth inclined substrate and droplet freezing on a horizontal substrate. We first simulated a non-freezing droplet sliding down an inclined surface, a configuration previously investigated by Park & Kumar (Reference Park and Kumar2017). As shown in Appendix A.1, figure 16, the final equilibrium droplet profile closely matches the corresponding result reported in figure 2(c) of Park & Kumar (Reference Park and Kumar2017). We further validated the freezing component of the model by reproducing the benchmark case of Zadražil et al. (Reference Zadražil, Stepanek and Matar2006). Appendix A.2, figure 17, illustrates the temporal evolution of the droplet interface,
$h$
(solid line), and the freezing front,
$s$
(dash-dotted line), showing good agreement with figure 15 of Zadražil et al. (Reference Zadražil, Stepanek and Matar2006). All simulations are performed on a one-dimensional mesh along the
$x$
-direction, spanning a domain of thirty-two dimensionless units, discretised using 25 601 grid points. Extensive domain-independence and grid-convergence tests were conducted to optimise the domain size and grid resolution, taking into account computational time. Appendix A.3 reports the grid-convergence results (figure 18), indicating that grids with 12 001, 25 601 and 30 001 points produce indistinguishable solutions. Additional details can be found from Williams et al. (Reference Williams, Karapetsas, Mamalis, Sefiane, Matar and Valluri2021), Kavuri et al. (Reference Kavuri, Karapetsas, Sharma and Sahu2023) and Wang et al. (Reference Wang, Karapetsas, Valluri and Inoue2024). Now, we describe the procedure used to determine the contact angles at both the substrate and the liquid–ice interface.
2.6. Calculation of contact angles
The apparent advancing (
$\theta _a$
) and receding (
$\theta _r$
) contact angles are obtained by locating the points of maximum slope on the droplet profile and drawing tangents from these points to the substrate before freezing begins (
$t/t_{\kern-1pt f} = 0$
), where
$t_{\kern-1pt f}$
denotes the total freezing time. As shown in figure 2(a), the point of maximum slope on the receding side is marked as A and that on the advancing side as B. The enlarged insets next to panel (a) highlight these points and illustrate how the corresponding contact angles are measured. After freezing is complete (
$t/t_{\kern-1pt f} = 1$
), the horizontal positions of points A and B are retained and relabelled as C and D, representing the receding and advancing sides, respectively. Since the locations of maximum slope after freezing are no longer close to the contact line, the final contact angles are evaluated at the same horizontal positions identified prior to freezing. The insets in figure 2(b) illustrate the procedure used to determine these angles. As noted earlier, we employ a penalty function (2.52) to preserve the precursor layer thickness far from the droplet and to maintain accuracy in predicting the freezing process. This method, however, limits the freezing rate when the liquid layer becomes comparable to the precursor thickness. As a result, the model becomes less accurate near the end of solidification and is unable to resolve the formation of the sharp tip. Recognising this limitation and noting that the tip forms only at the very final freezing stage, we estimate the tip angle using the local geometry of the ice–gas interface near the tip following Kavuri et al. (Reference Kavuri, Karapetsas, Sharma and Sahu2023). As shown in figure 2(b), tangents are drawn near the tip (labelled as E) and the corresponding inset illustrates this construction. For the set of parameters considered in figure 2, the resulting tip angle is approximately
$137^{\circ }$
, which closely matches the experimental value reported by Marin et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014).
Illustration of the procedure used to determine contact angles during (a) the initial and (b) final stages of droplet freezing. Insets at points A–D highlight the measurement locations, while the tip angle after complete freezing (point E) is shown in the right inset. The advancing (
$\theta _a$
) and receding (
$\theta _r$
) contact angles at
$t = 0$
(points A and B) are obtained by drawing tangents at the locations of maximum slope in the profile of
$h_x = \partial h/\partial x$
, following Schwartz & Eley (Reference Schwartz and Eley1998) and Espín & Kumar (Reference Espín and Kumar2015). The blue solid line in panel (a) represents
$0.1 h_x$
. At
$t = t_{\kern-1pt f}$
, tangents at the same locations (points C and D shown in panel b) yield the final contact angles. The tip angle near the cusp (point E) is approximately
$137^\circ$
. Here,
$\theta _{\textit{eq}} = 26.5^\circ$
,
$\alpha = 75^\circ$
and the remaining parameters are listed in table 2.
Illustration of the procedure used to determine the liquid–ice contact angle (
$\theta _{\textit{ice}}$
) at an early stage of freezing (
$t/t_{\kern-1pt f} = 0.1$
). The receding (
$\theta _{\textit{ice},r}$
) and advancing (
$\theta _{\textit{ice},a}$
) contact angles (points A and B) are identified by drawing tangents at the locations of maximum slope in the profile of
$h_x = \partial h/\partial x$
, following Schwartz & Eley (Reference Schwartz and Eley1998) and Espín & Kumar (Reference Espín and Kumar2015). The blue solid line represents
$0.1 h_x$
. Here,
$\theta _{\textit{eq}} = 26.5^\circ$
,
$\alpha = 75^\circ$
and the remaining parameters are provided in table 2.
Dimensionless groups and their typical values considered in the present study.
The contact angle at the liquid–ice interface (
$\theta _{\textit{ice}}$
) is calculated using a procedure analogous to that described earlier. The apparent advancing (
$\theta _{\textit{ice},a}$
) and receding (
$\theta _{\textit{ice},r}$
) angles are obtained by identifying the points of maximum slope within the liquid layer above the ice and drawing tangents from these points to the ice surface. As shown in the middle panel of figure 3, the points of maximum slope on the receding and advancing sides are labelled A and B, respectively. The corresponding magnified insets highlight these locations and illustrate the procedure used to determine the receding and advancing contact angles. Additionally, to clearly illustrate the advancing and receding contact angles, both the contact length (diameter) and the drop height are presented using consistent scaling in figures 2 and 3. Accordingly, we plot
$\epsilon h$
versus
$x$
rather than
$h$
versus
$x$
, since the contact length and drop height are scaled differently,
$(x, z, h) = l (\tilde {x}, \epsilon \tilde {z}, \epsilon \tilde {h})$
(see, (2.23)).
3. Results and discussion
We begin by examining the asymmetric shape of the droplet and analysing the contact angles at the two interfaces, where
$\theta$
denotes the droplet–substrate contact angle and
$\theta _{{\textit{ice}}}$
represents the liquid–ice interface contact angle for a droplet placed on an inclined surface undergoing freezing. For comparison, the freezing dynamics of a droplet on a horizontal substrate is also presented. In all simulations, the droplet is first allowed to attain its equilibrium shape prior to the onset of freezing. This is achieved by suppressing freezing by setting the Stefan number,
${\textit{Ste}} = 0$
. Once the droplet reaches a steady shape while moving on the substrate, indicating equilibrium, its geometry and the pressure distribution along the liquid–gas interface are extracted and used as initial conditions for the subsequent simulations in which freezing is activated. The dimensionless parameters used in the freezing simulations are
${\textit{Ste}} = 1.49 \times 10^{-3}$
,
$ \textit{Bo} = 0.4$
,
$A_{n} = 10$
,
$D_{s} = 0.9$
,
$\varLambda _{S} = 3.89$
,
$\varLambda _{W} = 675$
,
$Bi = 0.16$
,
$D_{w} = 2.94$
and
$\epsilon = 0.2$
, with their definitions summarised in table 2. This set of parameters corresponds to a water droplet exhibiting a contact angle of
$17.5^\circ$
(2.34) on a horizontal cold substrate maintained at
$T_c=-30^\circ$
. It is important to note that, although Knight (Reference Knight1967), Thiévenaz et al. (Reference Thiévenaz, Josserand and Séon2020a
), Huerre et al. (Reference Huerre, Josserand and Séon2025), Sarlin et al. (Reference Sarlin, Papa, Grivet, Rosenbaum, Huerre, Séon and Josserand2025) and Demmenie et al. (Reference Demmenie, Reus, Kolpakov, Woutersen, Bonn and Shahidzadeh2023, Reference Demmenie, Gorin, Kolpakov, Smith, Kellay and Bonn2025) report that the contact angle of water on ice near the melting temperature is approximately
$12^\circ$
, the liquid layer above the ice does not necessarily exhibit this equilibrium value once the recalescence phase has occurred. While the freezing front itself remains at the melting temperature, the unfrozen liquid layer can display larger apparent contact angles than
$12^\circ$
, as observed experimentally by Marin et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014), Zhang et al. (Reference Zhang, Zhao, Lv and Yang2016, Reference Zhang, Liu, Min and Wu2019), Jung et al. (Reference Jung, Tiwari and Poulikakos2012b
), Wang et al. (Reference Wang, Xu, Zhang, Zheng, Hao, He and Zhang2022). In the discussion that follows, we will investigate a wide range of contact angles and examine its effect on the freezing process. For completeness, we present droplet shapes for
$\theta _{\textit{eq}} = 12^\circ$
and compare them with
$\theta _{\textit{eq}} = 17.5^\circ$
in figure 19 of Appendix B.
Evolution of the droplet shape,
$h$
(solid lines), and the freezing front,
$s$
(dash-dotted lines), on a substrate with (a)
$\alpha = 0^\circ$
and (b)
$\alpha = 75^\circ$
. The normalised time,
$t/t_{\kern-1pt f}$
, corresponding to each row is indicated in panel (a). Here, the equilibrium contact angle is fixed at
$\theta _{\textit{eq}} = 17.5^\circ$
, with all other parameters listed in table 2. The corresponding total freezing time (
$t_{\kern-1pt f}$
) is reported in table 3.
Temporal evolution of (a, c) the droplet–substrate contact angle (
$\theta$
) and (b, d) the liquid–ice interface contact angle (
$\theta _{\textit{ice}}$
). Panels (a, b) correspond to a horizontal substrate (
$\alpha = 0^\circ$
), while panels (c, d) correspond to an inclined substrate (
$\alpha = 75^\circ$
). In each panel, the contact angles on the advancing (
$\theta _a$
) and receding (
$\theta _r$
) sides are depicted. The equilibrium contact angle is fixed at
$\theta _{\textit{eq}} = 17.5^\circ$
, with all other dimensionless parameters listed in table 2. The corresponding total freezing time (
$t_{\kern-1pt f}$
) is reported in table 3.
The temporal evolution of the droplet shape (
$h$
) and the freezing front (
$s$
) for substrates with
$\alpha = 0^\circ$
and
$\alpha = 75^\circ$
is shown in figures 4(a) and 4(b), respectively. It can be seen that on a horizontal substrate (figure 4
a), the droplet solidifies symmetrically with respect to the normal (dotted line) passing through its centre. In contrast, on an inclined substrate with
$\alpha = 75^\circ$
(figure 4
b), the initial droplet shape at
$t/t_{\kern-1pt f} = 0$
is asymmetric. This asymmetry persists during freezing, resulting in an asymmetrically solidified droplet, with the cusp formed at complete solidification (
$t/t_{\kern-1pt f} = 1$
) oriented towards the advancing side. Moreover, the final frozen morphology also deviates significantly from the initial droplet shape.
Figures 5(a, b) and 5(c, d) show the temporal evolution of the droplet–substrate contact angle (
$\theta$
) and the liquid–ice interface contact angle (
$\theta _{{\textit{ice}}}$
) for the same set of parameters used in figure 4. For a droplet freezing on a horizontal substrate, the substrate contact angle (
$\theta$
) remains nearly constant until the contact line freezes at
$t/t_{\kern-1pt f} = 10^{-3}$
. Beyond this point,
$\theta$
decreases and then remains constant until solidification is complete, as shown in figure 5(a). The liquid–ice interface angle (
$\theta _{\textit{ice}}$
), however, continues to increase until
$t/t_{f} = 0.65$
, by which time, the freezing front has already approached the droplet apex (figure 5
b). It can be seen that on a horizontal substrate, both
$\theta$
and
$\theta _{{\textit{ice}}}$
are symmetric on either side of the droplet and overlap, as evident in figures 5(a) and 5(b). In contrast, figure 5(c, d) illustrates that for the inclined substrate, both the advancing (
$\theta _a$
) and receding (
$\theta _r$
) contact angles initially increase slightly before dropping sharply once the contact line freezes. Throughout the freezing process, the advancing side consistently maintains a higher contact angle than the receding side (figure 5
c). The earlier decrease in
$\theta _r$
compared with
$\theta _a$
arises from the difference in height between the points of maximum slope on the advancing and receding sides of the droplet, where the contact angles are determined by drawing a tangent at the point of maximum slope, as described in § 2.6. Initially, the contact angles at the liquid–ice interface on both sides (
$\theta _{\textit{ice},a}$
and
$\theta _{\textit{ice},r}$
) match their respective substrate contact angles and increase as freezing progresses, converging near the end of the process (figure 5
d). This convergence indicates a loss of asymmetry between the advancing and receding sides, which explains why, as reported by Starostin et al. (Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2022, Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2023), the tip angle of frozen droplets remains nearly the same, around
$143^\circ \pm 11^\circ$
, on both horizontal and inclined substrates. This behaviour is expected since the effective Bond number (
$ \textit{Bo}_m = Bo \,{\mathrm{sin}}\alpha$
) for the remaining liquid towards the end of freezing becomes small, thus diminishing the effect of gravity on the liquid profile near the tip just before completion of freezing. With gravity playing a negligible role, the contact angles at the liquid–ice interface converge, leaving the tip angle largely independent of the inclination angle (
$\alpha$
).
Total dimensionless freezing time (
$t_{\kern-1pt f}$
) for different values of the inclination angle
$(\alpha )$
and equilibrium contact angle
$(\theta _{\textit{eq}})$
. Here,
$ \textit{Bo} = 0.4$
and the remaining dimensionless parameters are listed in table 2.
The equilibrium contact angle (
$\theta _{\textit{eq}}$
) estimated from (2.34) differs slightly from the apparent contact angles obtained by drawing a tangent at the point of maximum slope, as described in § 2.6 of the revised manuscript. This is illustrated in figures 5(a) and 5(c), where (2.34) gives
$\theta _{\textit{eq}} = 17.5^\circ$
, while the contact angle on the horizontal substrate, determined using the procedure in § 2.6, is approximately
$20^\circ$
. For the inclined substrate (
$\alpha = 75^\circ$
), the average of the advancing and receding angles,
$(\theta _a + \theta _r)/2$
, computed using the same method, is approximately
$13.5^\circ$
. The discrepancy arises from the absence of a sharply defined contact line. The apparent contact angle is therefore inferred rather than directly measured, introducing inherent ambiguity. Throughout this study, we consistently use the method in § 2.6 to estimate contact angles on both the substrate and the liquid–ice interface. In the following, we present a mechanistic framework to explain the observed asymmetry during droplet freezing by performing a large set of simulations in which the dimensionless parameters listed in table 2 are systematically varied. For convenience, tables 3–6 summarise the total dimensionless freezing time (
$t_{\kern-1pt f}$
) for different values of the inclination angles
$(\alpha )$
, equilibrium contact angles
$(\theta _{\textit{eq}})$
, Bond number
$(Bo)$
and Stefan number
$(Ste)$
.
Total dimensionless freezing time (
$t_{\kern-1pt f}$
) for different values of the inclination angle
$(\alpha )$
and Bond number
$(Bo)$
at
$\theta _{\textit{eq}} = 26.5^\circ$
. The remaining dimensionless parameters are listed in table 2.
Total dimensionless freezing time (
$t_{\kern-1pt f}$
) for different values of the equilibrium contact angle
$(\theta _{\textit{eq}})$
at
$\alpha =60^\circ$
. The remaining dimensionless parameters are listed in table 2.
Total dimensionless freezing time (
$t_{\kern-1pt f}$
) for different values of the the inclination angle
$(\alpha )$
and Stefan number
$(Ste)$
at
$\theta _{\textit{eq}}=26.5^\circ$
. The remaining dimensionless parameters are listed in table 2.
3.1. Effect of the sliding force
Effect of the effective Bond number (
$ \textit{Bo}_m = Bo \, {\mathrm{sin}}\alpha$
) on (a) the cusp angle
$\theta _{{Cusp}}$
and (b) the droplet–substrate contact angle hysteresis
$(\theta _a-\theta _r)$
at
$t/t_{\kern-1pt f} = 0$
(solid lines) and
$t/t_{\kern-1pt f} = 1$
(dash-dotted lines) for different equilibrium contact angles. The remaining dimensionless parameters are listed in table 2. The total freezing times (
$t_{\kern-1pt f}$
) for different values of
$\alpha$
and
$ \textit{Bo}$
are provided in tables 3 and 4.
The freezing behaviour of a sliding droplet is strongly influenced by both the Bond number (
$ \textit{Bo}$
) and the inclination angle (
$\alpha$
), as they together determine the component of gravity acting tangentially along the inclined surface. To capture this combined effect, we employ the effective Bond number (
$ \textit{Bo}_m = Bo \, {\mathrm{sin}}\alpha$
). This parameter represents the tangential gravitational contribution that drives droplet motion and provides a resultant measure of the coupled influence of gravity and inclination of the substrate on the freezing dynamics. The effect of
$ \textit{Bo}_m$
on the cusp angle,
$\theta _{ {Cusp}}$
(defined as the angle between the frozen tip formed at the end of solidification and the
$z$
-axis, as shown in figure 1) and the contact angle hysteresis (
$\theta _a-\theta _r$
) measured on the droplet–substrate interface at
$t/t_{\kern-1pt f} = 0$
and
$t/t_{\kern-1pt f} = 1$
for different values of the equilibrium contact angle
$(\theta _{\textit{eq}})$
are illustrated in figure 6(a,b).
The substrate wettability is characterised by
$\theta _{\textit{eq}}$
, with higher
$\theta _{\textit{eq}}$
corresponding to lower wettability. It can be seen from figure 6(a) that the cusp angle of the droplet increases with
$ \textit{Bo}_m$
for both values of
$\theta _{\textit{eq}}$
considered. This behaviour can be attributed to the continuous sliding of the liquid droplet and the increasing asymmetry at higher inclination angles, which leads to a larger deviation of the frozen tip from the normal (i.e. a higher
$\theta _{ {Cusp}}$
). This trend is observed for both
$\theta _{\textit{eq}}=17.5^\circ$
and
$\theta _{\textit{eq}}=26.5^\circ$
. For
$\theta _{\textit{eq}}=26.5^\circ$
, the droplet–substrate contact angle hysteresis (
$\theta _a-\theta _r$
), measured before freezing (
$t/t_{\kern-1pt f}=0$
) and after complete solidification (
$t/t_{\kern-1pt f}=1$
), exhibits a non-monotonic dependence on
$ \textit{Bo}_m$
, as shown in figure 6(b). This non-monotonic dependence is attributed to the long tail of the droplet that forms and eventually freezes when the sliding force is high. Figure 7(a–c) shows the evolution of droplet shapes for effective Bond numbers
$ \textit{Bo}_m = 0.19$
, 0.92 and 1.2, respectively. For each
$ \textit{Bo}_m$
, the droplet profiles are presented at the onset of freezing (
$t/t_{\kern-1pt f} = 0$
), at an intermediate stage (
$t/t_{\kern-1pt f} = 0.1$
) and at the completion of freezing (
$t/t_{\kern-1pt f} = 1$
). It can be seen that as
$ \textit{Bo}_m$
increases from 0.92 to 1.2, the initial droplet shape becomes increasingly elongated, with pronounced swelling on the advancing side and tapering on the receding side. The stronger sliding force at higher
$ \textit{Bo}_m$
enhances droplet elongation, leading to a reduction in the droplet–substrate contact angle hysteresis (
$\theta _a-\theta _r$
) at both the initial and final stages of freezing.
Droplet shapes at different normalised times for various effective Bond numbers
$ \textit{Bo}_m \, (= Bo \, {\mathrm{sin}}\alpha )$
: (a)
$ \textit{Bo}_m = 0.19$
, (b)
$ \textit{Bo}_m = 0.92$
and (c)
$ \textit{Bo}_m = 1.2$
. At
$t/t_{\kern-1pt f} = 1.0$
, the freezing front coincides with the droplet surface. The equilibrium contact angle is
$\theta _{\textit{eq}} = 26.5^\circ$
. The remaining dimensionless parameters and the corresponding total freezing time (
$t_{\kern-1pt f}$
) are provided in tables 2 and 4, respectively.
Temporal evolution of (a) the advancing (
$\theta _{a}$
, red dashed lines) and receding (
$\theta _{r}$
, black solid lines) contact angles on the substrate; (b) the droplet shape,
$h$
(solid lines), and the freezing front,
$s$
(dash-dotted lines), for a droplet freezing on an inclined substrate with
$\alpha = 75^\circ$
; and (c) variations of the capillary velocity component (
$u_{\textit{ca}}$
), gravitational velocity component (
$u_{g}$
) and average velocity (
$u_{a}$
) along the
$x$
-axis at a representative early freezing stage (
$t/t_{f} = 0.01$
) for a droplet on a substrate inclined at
$\alpha = 75^\circ$
, where the dotted line indicates
$u = 0$
. The equilibrium contact angle is
$\theta _{\textit{eq}} = 17.5^\circ$
. The remaining dimensionless parameters and the corresponding total freezing time (
$t_{\kern-1pt f}$
) are provided in tables 2 and 3, respectively.
Interestingly, for the substrate with the highest wettability (
$\theta _{\textit{eq}} = 17.5^\circ$
), the pre-freezing droplet–substrate contact angle hysteresis (
$(\theta _a - \theta _r)$
at
$t/t_{\kern-1pt f} = 0$
) exceeds the post-freezing value (
$(\theta _a - \theta _r)$
at
$t/t_{\kern-1pt f} = 1$
) for
$ \textit{Bo}_m \gt 0.28$
, corresponding to an inclination angle
$\alpha \gt 45^\circ$
for
$ \textit{Bo} = 0.4$
, as shown in figure 6(b). To understand this behaviour, in figure 8(a), we plot the temporal evolution of the advancing (
$\theta _a$
) and receding (
$\theta _r$
) contact angles for substrate inclinations of
$\alpha = 15^\circ$
,
$45^\circ$
and
$75^\circ$
. It is observed that
$\theta _a$
increases while
$\theta _r$
decreases as the substrate inclination increases. However, the increase in
$\theta _a$
from
$\alpha = 45^\circ$
to
$75^\circ$
is noticeably smaller than the increase observed when the inclination increases from
$\alpha = 15^\circ$
to
$45^\circ$
. For the case of
$\alpha = 75^\circ$
, the evolution of the droplet shape (
$h$
) and the freezing front (
$s$
), shown in figure 8(b), indicates that the advancing contact line remains pinned throughout freezing, whereas the receding contact line progressively moves downstream. Consequently, the increase in
$\theta _a$
from
$\alpha = 45^\circ$
to
$75^\circ$
in figure 8(a) is relatively small. To further elucidate this behaviour, we evaluate the vertically averaged velocity
$(u_a(x))$
and the bulk average velocity
$(\bar {u}_a)$
of the unfrozen liquid droplet using the following expressions:
\begin{align} u_{a} (x) = \frac {\int _{s}^{h}u\,{\rm d}z}{\int _{s}^{h}\,{\rm d}z} &= -\frac {1}{3}\left (p_{x}-\frac {Bo \, {\mathrm{sin}}\alpha }{\epsilon }\right )(h-s)^{2}, \end{align}
\begin{align} \bar {u}_{a} = \frac {\int _{x_{cl}}^{x_{cr}}\int _{s}^{h}u\,{\rm d}z\,{\rm d}x}{\int _{x_{cl}}^{x_{cr}}\int _{s}^{h}\,{\rm d}z\,{\rm d}x} &= -\frac {1}{3}\frac {\int _{x_{cl}}^{x_{cr}}\left (p_{x}-\frac {Bo \, {\mathrm{sin}}\alpha }{\epsilon }\right )(h-s)^{3}\,{\rm d}x}{\int _{x_{cl}}^{x_{cr}}(h-s)\,{\rm d}x}. \end{align}
Here,
$x_{cl}$
and
$x_{cr}$
denote the left and right contact lines of the droplet (see figure 1). To gain further insight, the velocities are analysed by separating the contributions from capillary forces and gravitational forces. Within this framework, both the vertically averaged velocity
$(u_a(x))$
and the centre-of-mass velocity
$(\bar {u}_a)$
of the unfrozen liquid droplet are decomposed into these two distinct components.
\begin{align} u_a(x) &= u_{\textit{ca}}(x) + u_g(x) = -\frac {1}{3} p_x (h - s)^2 + \frac {1}{3}\left (\frac {Bo {\mathrm{sin}}\alpha }{\epsilon }\right )(h - s)^2, \nonumber\\ u_{\textit{ca}}(x) &= -\frac {1}{3} p_x (h - s)^2, \quad u_g(x) = \frac {1}{3}\left (\frac {Bo {\mathrm{sin}}\alpha }{\epsilon }\right )(h - s)^2,\\[-28pt] \nonumber \end{align}
\begin{align} \bar {u}_a &= \bar {u}_{\textit{ca}} + \bar {u}_g, \\[0pt] \nonumber \end{align}
where
\begin{equation} \bar {u}_{\textit{ca}} = -\frac {1}{3} \frac {\int _{x_{cl}}^{x_{cr}} p_x (h-s)^3 \, {\rm d}x} {\int _{x_{cl}}^{x_{cr}} (h-s) \, {\rm d}x}, \quad \bar {u}_g = \frac {1}{3} \frac {\int _{x_{cl}}^{x_{cr}} \frac {Bo {\mathrm{sin}}\alpha }{\epsilon } (h-s)^3 \, {\rm d}x} {\int _{x_{cl}}^{x_{cr}} (h-s) \, {\rm d}x}. \end{equation}
Temporal evolution of (a) the advancing (
$\theta _{\textit{ice},a}$
, red dashed lines) and receding (
$\theta _{\textit{ice},r}$
, black solid lines) contact angles at the liquid–ice interface, (b) the difference between the advancing and receding contact angles at the liquid–ice interface (
$\Delta \theta _{\textit{ice}} = \theta _{\textit{ice},a} - \theta _{\textit{ice},r}$
), and (c) the bulk velocity (
$\bar {u}_a$
) during the early stage of freezing for a droplet on substrates with different inclination angles (
$\alpha$
). The equilibrium contact angle is
$\theta _{\textit{eq}}=17.5^\circ$
. The remaining dimensionless parameters and the values of the total freezing time (
$t_{\kern-1pt f}$
) for different inclination angles are provided in tables 2 and 3, respectively.
The capillary contributions
$u_{\textit{ca}}$
and
$\bar {u}_{\textit{ca}}$
capture the influence of interfacial forces that act to redistribute fluid within the medium, while the gravitational components
$u_{g}$
and
$\bar {u}_g$
account for the fluid motion induced by the weight of the liquid. In figure 8(c), we plot the capillary (
$u_{\textit{ca}}$
) and gravitational (
$u_{g}$
) components of the average velocity (
$u_{a}$
) for an inclination of
$\alpha = 75^\circ$
. As demonstrated (and also as predicted by 3.3), gravitational forces become significant only in regions where the droplet is relatively thick, whereas capillary forces are highly sensitive to the local curvature of the droplet interface. On the advancing side, capillary and gravitational contributions effectively counteract each other, resulting in the pinning of the contact line. Additionally, once the first layer of ice forms, the remaining liquid rests on a mesa region, further restricting the motion of the droplet on the advancing side. In contrast, on the receding side, capillary forces dominate and drive the retraction of the contact line. It is also observed that variations in the advancing and receding contact angles occur only during the early stages of freezing (up to
$t/t_{\kern-1pt f} = 0.01$
). Beyond this point, both angles remain nearly constant (see figure 8
a) as the liquid adjacent to the substrate solidifies.
In addition to the contact line between the ice and the solid substrate, two additional contact lines form at the liquid–ice interface. Figure 9(a) shows the temporal evolution of the corresponding contact angles on the advancing (
$\theta _{\textit{ice},a}$
) and receding sides (
$\theta _{\textit{ice},r}$
), respectively, for different substrate inclinations (
$\alpha$
). For a horizontal substrate (
$\alpha = 0^\circ$
), the liquid–ice contact angles on both sides of the droplet remain equal throughout the freezing process. In contrast, for inclinations of
$\alpha = 15^\circ$
,
$45^\circ$
and
$75^\circ$
, the advancing and receding contact angles differ, with
$\theta _{\textit{ice},a} \gt \theta _{\textit{ice},r}$
, and remain unequal until
$t/t_{\kern-1pt f} \approx 0.65$
. Figure 9(b) illustrates that the difference
$\Delta \theta _{\textit{ice}} = \theta _{\textit{ice},a} - \theta _{\textit{ice},r}$
increases with
$\alpha$
, reflecting the larger tangential gravitational force acting on the droplet at higher inclinations. This trend is consistent with the temporal evolution of the bulk velocity
$\bar {u}_a$
in figure 9(c), which increases with
$\alpha$
, but vanishes after
$t/t_{\kern-1pt f} = 0.1$
. For all non-zero inclinations, the difference between
$\theta _{\textit{ice},a}$
and
$\theta _{\textit{ice},r}$
increases up to
$t/t_{\kern-1pt f} = 0.01$
, primarily due to the initial contact angle hysteresis caused by the sliding motion of the droplet prior to freezing, combined with the pinning behaviour of the contact-line on the advancing side and capillary-driven retraction on the receding side (see figure 8
b, c and related discussion). As freezing proceeds, the loss of liquid volume from solidification reduces the effect of gravity, consequently decreasing the velocity of the centre-of-mass
$\bar {u}_a$
. Capillary forces then become dominant, and once the remaining liquid volume is sufficiently small, gravitational effects become negligible. This leads to a nearly symmetric liquid–ice interface, with the advancing and receding contact angles equalising, as shown in figure 9(a).
3.2. Interplay between wettability and inclination
Variation of (a) the cusp angle (
$\theta _{{Cusp}}$
) at the droplet tip, measured relative to the
$z$
-axis, as a function of substrate inclination (
$\alpha$
) for different equilibrium contact angles. Here, the experimental results from Kumar et al. (Reference Kumar, Subudhi and Das2025) are included for comparison. (b) Variation of the difference between advancing and receding contact angles (
$\theta _a - \theta _r$
) with inclination angle (
$\alpha$
) for different equilibrium contact angles at
$t/t_{\kern-1pt f} = 0$
(black solid lines) and
$t/t_{\kern-1pt f} = 1$
(magenta dashed lines). (c) Variation of
$\theta _{{Cusp}}$
with
$\theta _{\textit{eq}}$
. (d) Variation of (
$\theta _a - \theta _r$
) with
$\theta _{\textit{eq}}$
at
$t/t_{\kern-1pt f} = 0$
(black solid line) and
$t/t_{\kern-1pt f} = 1$
(magenta dashed line) on the left axis, with the change in
$(\theta _a - \theta _r)$
between
$t/t_{\kern-1pt f} = 1$
and
$t/t_{\kern-1pt f} = 0$
shown on the right axis (blue solid line). In panels (c) and (d),
$\alpha = 60^\circ$
. The remaining dimensionless parameters are listed in table 2. The values of the total freezing time (
$t_{\kern-1pt f}$
) for different inclination angles and equilibrium contact angles are provided in tables 3 and 5, respectively.
In the previous subsection, we focused on droplets with
$\theta _{\textit{eq}} = 17.5^\circ$
. Here, we examine the effect of inclination (
$\alpha$
) and wettability (
$\theta _{\textit{eq}}$
) on the asymmetric freezing of a sliding drop, keeping all other parameters fixed. In figure 10(a, b), we examine the dependence of the cusp angle (
$\theta _{{Cusp}}$
), and the difference between the advancing and receding contact angles
$(\theta _a-\theta _r)$
at
$t/t_{\kern-1pt f}=0$
(black solid lines) and
$t/t_{\kern-1pt f}=1$
(red dashed lines) as a function of the inclination angle (
$\alpha$
) for droplets with
$\theta _{\textit{eq}} = 17.5^\circ$
,
$21.2^\circ$
,
$26.5^\circ$
and
$30.6^\circ$
. As observed in figure 10(a) for a given inclination angle,
$\theta _{{Cusp}}$
decreases with decreasing substrate wettability. Interestingly, the rate of decrease progressively diminishes and by
$\theta _{\textit{eq}}=30.6^\circ$
, the curves exhibit a noticeable reduction in slope over the range considered, although the system remains in a transitional regime. This trend is also illustrated in figure 10(c), where we plot the cusp angle as a function of the equilibrium contact angle for an inclination angle
$\alpha =60^\circ$
. The corresponding values of
$\theta _{ {Cusp}}$
are comparable in magnitude to those reported by Kumar et al. (Reference Kumar, Subudhi and Das2025) for a sessile water droplet on an inclined copper substrate at different inclination angles and an equilibrium contact angle of
$70^\circ$
; however, this comparison is qualitative, as the physical and wettability conditions differ. This behaviour can be attributed to the fact that as the wettability decreases, the droplet tends to retain its original shape due to the increased resistance to spreading. As a result, the influence of gravitational forces acting on the sliding droplet becomes less effective in deforming its shape.
Temporal evolution of droplet shape (
$h$
, solid) and freezing front (
$s$
, dash-dotted) on an inclined substrate (
$\alpha = 75^\circ$
) for (a)
$\theta _{\textit{eq}} = 17.5^\circ$
and (b)
$\theta _{\textit{eq}} = 26.5^\circ$
. The remaining dimensionless parameters are listed in table 2 and the values of the total freezing time (
$t_{\kern-1pt f}$
) for different equilibrium contact angles are provided in table 5.
Evolution of the droplet shape,
$h$
(solid lines), and the freezing front position,
$s$
(dash-dotted lines) for (a)
$\theta _{\textit{eq}} = 15.3^\circ$
and (b)
$\theta _{\textit{eq}} = 21.2^\circ$
. The normalised time,
$t/t_{\kern-1pt f}$
, corresponding to each row is indicated in panel (a). Here,
$\alpha = 60^\circ$
, and the remaining dimensionless parameters and the corresponding total freezing times (
$t_{\kern-1pt f}$
) are provided in tables 2 and 5, respectively.
In figure 10(b), we observe that the difference between the advancing (
$\theta _{a}$
) and receding (
$\theta _{r}$
) contact angles increases with the inclination angle (
$\alpha$
). For larger equilibrium contact angles (
$\theta _{\textit{eq}} = 21.2^\circ$
,
$26.5^\circ$
and
$30.6^\circ$
), the hysteresis (
$\theta _{a}-\theta _{r}$
) is larger after freezing (
$t/t_{f} = 1$
) than before freezing (
$t/t_{f} = 0$
). This behaviour is not observed for the most wettable case (
$\theta _{\textit{eq}} = 17.5^\circ$
). To further examine the role of wettability, figure 10(d) shows the variation of
$\theta _{a}-\theta _{r}$
with
$\theta _{\textit{eq}}$
for
$\alpha =60^\circ$
. In the unfrozen state (
$t/t_{\kern-1pt f}=0$
), the contact angle hysteresis (
$\theta _a-\theta _r$
) decreases with increasing
$\theta _{\textit{eq}}$
(i.e. with decreasing wettability). However, the opposite trend is observed at the end of the freezing process (
$t/t_{\kern-1pt f}=1$
). We also plot
$\Delta (\theta _a-\theta _r)=(\theta _a-\theta _r)_{t/t_{\kern-1pt f}=0}-(\theta _a-\theta _r)_{t/t_{\kern-1pt f}=1}$
in figure 10(d). The value of
$\Delta (\theta _a - \theta _r)$
is negative for highly wettable substrates and positive for less hydrophilic ones, and it increases monotonically with
$\theta _{\textit{eq}}$
. To clarify this behaviour, we compare the droplet profile (
$h$
, solid) and the freezing front (
$s$
, dash-dotted) at
$\alpha = 75^\circ$
for
$\theta _{\textit{eq}} = 17.5^\circ$
and
$\theta _{\textit{eq}} = 26.5^\circ$
in figures 11(a) and 11(b). The droplet with
$\theta _{\textit{eq}}=17.5^\circ$
becomes more flattened at high inclination, whereas for
$\theta _{\textit{eq}}=26.5^\circ$
, the liquid shifts more under gravity. This produces a thinner ice layer on the receding side at
$t/t_{\kern-1pt f}=0.01$
and leads to
$(\theta _a-\theta _r)_{t/t_{\kern-1pt f}=1}\gt (\theta _a-\theta _r)_{t/t_{\kern-1pt f}=0}$
for
$\theta _{\textit{eq}}=26.5^\circ$
, while the opposite behaviour occurs for
$\theta _{\textit{eq}}=17.5^\circ$
. Keeping the inclination constant, the next subsection presents a detailed analysis of the influence of wettability on the asymmetric freezing of a sliding drop.
Temporal evolution of (a) the advancing (
$\theta _a$
, red dashed lines) and receding (
$\theta _r$
, black solid lines) contact angles on the substrate, (b) the advancing (
$\theta _{\textit{ice},a}$
, red dashed lines) and receding (
$\theta _{\textit{ice},r}$
, black solid lines) contact angles at the liquid–ice interface, and (c) the difference between the advancing and receding liquid–ice contact angles (
$\Delta \theta _{\textit{ice}}$
), for
$\theta _{\textit{eq}}=15.3^\circ$
and
$21.2^\circ$
. Here,
$\alpha = 60^\circ$
, and the remaining dimensionless parameters and the corresponding total freezing times (
$t_{\kern-1pt f}$
) are given in tables 2 and 5, respectively.
Temporal evolution of (a) the bulk velocity (
$\bar {u}_a$
) during the early stages of freezing for
$\theta _{\textit{eq}}=15.3^\circ$
and for
$\theta _{\textit{eq}}=21.2^\circ$
, (b) the capillary (
$\bar {u}_{\textit{ca}}$
) and gravitational (
$\bar {u}_g$
) components of the bulk velocity (
$\bar {u}_a$
) for
$\theta _{\textit{eq}}=15.3^\circ$
, and (c) the location of the maximum height of the drop (
$x_m$
) for a droplet on substrates with different equilibrium contact angles. Here,
$\alpha = 60^\circ$
, and the remaining dimensionless parameters and the corresponding total freezing times (
$t_{\kern-1pt f}$
) are listed in tables 2 and 5, respectively.
To further understand this behaviour, we examine two cases for
$\alpha = 60^\circ$
: one with a negative
$\Delta (\theta _a-\theta _r)$
(
$\theta _{\textit{eq}}=15.3^\circ$
) and another with a positive value (
$\theta _{\textit{eq}}=21.2^\circ$
). Figure 12 shows the temporal evolution of the droplet shape,
$h$
(solid lines), and the freezing front position,
$s$
(dash-dotted lines), for these cases. For the more hydrophilic substrate (
$\theta _{\textit{eq}}=15.3^\circ$
), the droplet shape before freezing exhibits a local concavity (negative curvature) near the receding side (see figures 12
a and 12
b) at
$t/t_{\kern-1pt f} = 0$
). This arises from the imbalance between gravity and surface tension: the tangential component of gravity stretches the droplet along the incline, while surface tension tends to restore a compact shape, leading to tail formation. This effect is stronger on hydrophilic substrates because enhanced spreading increases the wetted area and flattens the liquid–air interface, especially near the trailing edge. The resulting reduction in interfacial curvature lowers the capillary pressure resisting deformation, allowing gravity to elongate the droplet more easily. As freezing proceeds, the front advances into this concave trailing region, so the receding contact angle of the liquid–ice interface,
$\theta _{\textit{ice},r}$
, remains smaller than in the less hydrophilic case. In contrast, the frozen tail becomes thicker in the more hydrophilic case due to the faster freezing of the thin trailing liquid layer, producing a larger ice–substrate receding contact angle
$\theta _r$
(see figures 12
a and 12
b at
$t/t_{\kern-1pt f} = 1$
). Consequently,
$\Delta (\theta _a-\theta _r)\lt 0$
for sufficiently hydrophilic substrates.
A similar behaviour can also be observed in figure 13(a), which shows the temporal evolution of the contact angles at the droplet–substrate interface. For
$\theta _{\textit{eq}}=15.3^\circ$
, the difference between the advancing (
$\theta _a$
) and receding (
$\theta _r$
) contact angles is larger before the start of freezing (
$t/t_{\kern-1pt f} = 0$
) than after freezing is completed (
$t/t_{\kern-1pt f} = 1$
). In contrast, for
$\theta _{\textit{eq}}=21.2^\circ$
, the difference between
$\theta _a$
and
$\theta _r$
is greater after freezing (
$t/t_{\kern-1pt f} = 1$
) than before freezing begins (
$t/t_{\kern-1pt f} = 0$
). We further examine the temporal evolution of the receding (
$\theta _{\textit{ice},r}$
) and advancing (
$\theta _{\textit{ice},a}$
) contact angles of the liquid–ice interface for
$\theta _{\textit{eq}} = 15.3^\circ$
and
$\theta _{\textit{eq}} = 21.2^\circ$
in figure 13(b). For both cases,
$\theta _{\textit{ice},r}$
and
$\theta _{\textit{ice},a}$
increase slightly at early times. For
$\theta _{\textit{eq}} = 15.3^\circ$
, when the freezing front reaches the region of negative curvature at
$t/t_{\kern-1pt f} \approx 0.03$
(see figure 12
a),
$\theta _{\textit{ice},r}$
decreases, followed shortly by a drop in
$\theta _{\textit{ice},a}$
. Both angles then recover and converge near the end of freezing. In contrast, for
$\theta _{\textit{eq}}=21.2^\circ$
, both contact angles increase monotonically and converge only towards the end of freezing. This behaviour is reflected in figure 13(c), where
$\Delta \theta _{\textit{ice}}$
for
$\theta _{\textit{eq}}=15.3^\circ$
exhibits a pronounced maximum at
$t/t_{\kern-1pt f}\approx 0.03$
before decaying to zero, whereas for
$\theta _{\textit{eq}}=21.2^\circ$
, it increases smoothly to an early peak (
$t/t_{\kern-1pt f}\approx 0.004$
) and then decreases towards zero.
To rationalise this behaviour, we examine the bulk velocity (
$\bar {u}_a$
) of the unfrozen liquid for
$\theta _{\textit{eq}} = 15.3^\circ$
and
$21.2^\circ$
(figure 14
a). While
$\bar {u}_a$
decays to zero after
$t/t_{\kern-1pt f} \approx 0.1$
in both cases, it is larger for
$\theta _{\textit{eq}} = 21.2^\circ$
, explaining the pronounced thinning of the frozen receding edge (figure 12) and the larger post-freezing contact-angle contrast (figure 10
d). At
$t/t_{\kern-1pt f} \approx 0.03$
,
$\bar {u}_a$
briefly becomes negative for
$\theta_{eq} = 15.3$
as the freezing front passes a concave region, reducing capillary pressure at the receding edge and driving liquid towards it. The capillary (
$\bar {u}_{\textit{ca}}$
) and gravitational (
$\bar {u}_g$
) contributions (figure 14
b) show that liquid temporarily accumulates near the receding edge, opposite to gravity. This is reflected in the
$x$
-location of the maximum droplet height
$x_m$
(figure 14
c), which briefly decreases before resuming motion in the gravitational direction after
$t/t_{\kern-1pt f} \approx 0.05$
.
3.3. Effect of the Stefan number (
${\textit{Ste}}$
)
Effect of the Stefan number (
${\textit{Ste}}$
) on (a) the cusp angle (
$\theta _{ {Cusp}}$
) and (b) the difference between the advancing and receding contact angles
$(\theta _a-\theta _r)$
at
$t/t_{\kern-1pt f}=1$
, for a droplet on substrates with different inclination angles (
$\alpha$
). Here, the equilibrium contact angle is
$\theta _{\textit{eq}}=26.5^\circ$
. The remaining dimensionless parameters and the corresponding total freezing times (
$t_{\kern-1pt f}$
) are listed in tables 2 and 6, respectively.
Finally, we examine the effect of the Stefan number (
${\textit{Ste}}$
) on the cusp angle (
$\theta _{{Cusp}}$
) and the droplet–substrate contact-angle hysteresis
$(\theta _a-\theta _r)$
, evaluated at
$t/t_{\kern-1pt f}=0$
and
$t/t_{\kern-1pt f}=1$
, for different substrate inclinations (
$\alpha$
). As shown in figure 15(a), the cusp angle
$\theta _{{Cusp}}$
decreases monotonically with increasing Stefan number. A larger
${\textit{Ste}}$
implies that a greater amount of sensible heat is available relative to the latent heat of solidification, enabling more efficient removal of latent heat and thereby accelerating the freezing process. The faster advance of the freezing front leaves less time for the unfrozen liquid to slide downslope under gravity, resulting in the formation of a smaller cusp angle. As expected, prior to freezing (
$t/t_{\kern-1pt f} = 0$
),
$(\theta _a - \theta _r)$
is independent of the Stefan number. In contrast, after freezing is complete (
$t/t_{\kern-1pt f}=1$
),
$(\theta _a-\theta _r)$
decreases with increasing
${\textit{Ste}}$
, as the enhanced freezing rate at higher Stefan numbers more effectively suppresses the sliding motion of the unfrozen liquid. Consistent with this interpretation, an increase in the Stefan number leads to a substantial reduction in the dimensionless total freezing time, which decreases from
$t_{\kern-1pt f}=344$
for
${\textit{Ste}}=1.49\times 10^{-3}$
to
$t_{\kern-1pt f}=5$
for
${\textit{Ste}}=0.1$
, for
$\alpha =45^\circ$
and
$ \textit{Bo}=0.4$
.
4. Conclusions
We investigate the asymmetric freezing dynamics of a liquid droplet on an inclined cold substrate through numerical simulations based on the lubrication approximation. The model has been validated against previously published numerical and experimental studies for both non-freezing and freezing droplets on cold substrates, demonstrating good quantitative agreement. In particular, we observe that the predicted tip angle is consistent with experimental observations for a horizontal substrate reported by Marin et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014), as well as for inclined substrates in the studies of Starostin et al. (Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2022, Reference Starostin, Strelnikov, Dombrovsky, Shoval, Gendelman and Bormashenko2023). For a fixed inclination angle, the model predicts that
$\theta _{ {Cusp}}$
, the angle between the cusp singularity and the vertical, decreases with decreasing wettability, with the reduction gradually approaching saturation. At
$\theta _{\textit{eq}} = 30.6^\circ$
, the predictions converge to a lower bound that closely matches the experimental measurements of Kumar et al. (Reference Kumar, Subudhi and Das2025) for a sessile water droplet on an inclined copper substrate (
$\theta _{\textit{eq}} \approx 70^\circ$
). The computed frozen droplet shapes qualitatively capture key experimental features reported in the literature (Anderson et al. Reference Anderson, Worster and Davis1996; Hu & Jin Reference Hu and Jin2010; Marin et al. Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014; Zhang et al. Reference Zhang, Liu, Min and Wu2019), including volume expansion, the evolution of the freezing front and cusp formation. Furthermore, the equilibrium droplet shapes prior to freezing are in good agreement with the numerical results of Park & Kumar (Reference Park and Kumar2017) and are also consistent with experimental observations by Podgorski et al. (Reference Podgorski, Flesselles and Limat2001).
By systematically varying the substrate inclination, wettability, effective Bond number and Stefan number, we elucidate the coupled roles of gravity, capillarity and solidification in governing droplet motion, interfacial deformation and the resulting frozen morphology. Our results demonstrate that droplet motion prior to and during freezing fundamentally alters the final structure, thereby challenging conclusions drawn from studies of sessile droplets. While inclination and wettability primarily control the degree of asymmetry, the extent of deformation can also be influenced by substrate and ambient temperatures, an aspect that warrants further investigation. The interplay between gravity and capillarity emerges as the dominant factor shaping the frozen structure. Inclination and wettability jointly dictate the degree of asymmetry, while freezing kinetics modulate the extent of deformation. Highly wetting substrates exhibit unexpected behaviour, including transient motion opposite to gravity, underscoring the complexity of early-stage dynamics.
Notably, the model predicts a transient motion opposite to gravity on highly wetting substrates during the early stages of freezing, a phenomenon that, to our knowledge, has not yet been reported experimentally and may motivate future investigations. By decomposing the liquid motion into capillary- and gravity-driven components, we provide mechanistic insight into contact-line pinning, receding-edge thinning and the emergence of asymmetric ice–liquid contact angles. An increase in the Stefan number accelerates freezing, suppresses sliding-induced deformation, and reduces both the cusp angle and the contrast in post-freezing contact angles. These findings have practical implications for anti-icing strategies, thermal management, and surface design in aerospace and energy systems. Future work should extend this framework to three-dimensional droplets, textured or chemically heterogeneous substrates and coupled evaporation–freezing scenarios to better capture real-world conditions.
Acknowledgements
We sincerely thank the anonymous reviewers for their insightful comments and constructive suggestions, which significantly improved the quality and depth of the manuscript.
Funding
This work was supported by the Suzuki Next Bharat Fellowship Program (Grant No. NBV/CHE/F011/2025-26/S401). G.K. would also like to acknowledge the ThermEnTrans project funded by the European Union’s Horizon Europe research and innovation programme under grant agreement ID: 101236597.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation of the numerical model
A.1. Comparison with Park & Kumar (2017)
We validate our model by first simulating a non-freezing droplet sliding on a smooth inclined substrate, a configuration previously examined by Park & Kumar (Reference Park and Kumar2017). Figure 16 compares the equilibrium droplet shape obtained from our simulations at an inclination angle of
$\alpha = 25^\circ$
with the corresponding case shown in figure 2(c) of Park & Kumar (Reference Park and Kumar2017). To ensure a consistent comparison, freezing effects are neglected by setting
${\textit{Ste}} = 0$
. The remaining dimensionless parameters are
$\theta _{\textit{eq}} \approx 20^\circ$
,
$ \textit{Bo} = 0.9$
,
$A_n = 121.8$
and
$\epsilon = 0.44$
. Since the exact values of the scaled Hamaker constant (
$A_n$
) and the aspect ratio (
$\epsilon$
) are not reported by Park & Kumar (Reference Park and Kumar2017), this missing information may account for the minor discrepancies observed between the two results.
Comparison of the equilibrium shape of a droplet sliding on a substrate inclined at
$\alpha = 25^\circ$
, as obtained from the present simulations, with the corresponding results reported by Park & Kumar (Reference Park and Kumar2017). For comparison with Park & Kumar (Reference Park and Kumar2017), the effect of freezing is neglected in the present formulation by setting
${\textit{Ste}} = 0$
. The other dimensionless parameters are
$\theta _{\textit{eq}} \approx 20^\circ$
,
$ \textit{Bo} = 0.9$
,
$A_{n} = 121.8$
and
$\epsilon = 0.44$
.
Evolution of the freezing front,
$s$
(dash-dotted lines), and the shape,
$h$
(solid line), of the droplet undergoing freezing on a horizontal (
$\alpha = 0^\circ$
) cold substrate. The dimensionless parameters considered in the simulations are
$\epsilon = 0.2$
,
${\textit{Ste}} = 0.04$
,
$T_{v} = 0.5$
,
$A_{n} = 8$
,
$D_{s} = \varLambda _{S} = Bi = Bo = 1$
and
$D_{w} = 0$
. This set of parameters corresponds to those used by Zadražil et al. (Reference Zadražil, Stepanek and Matar2006). Here, the dimensionless total freezing time of the droplet is found to be
$t_{f} = 14$
, which is in good agreement with the value reported by Zadražil et al. (Reference Zadražil, Stepanek and Matar2006).
Evolution of the freezing front,
$s$
(dash-dotted lines), and the droplet shape,
$h$
(solid lines), on a cold substrate inclined at
$\alpha = 75^\circ$
. The remaining dimensionless parameters are
$ \textit{Bo} = 0.4$
,
${\textit{Ste}} = 1.49 \times 10^{-3}$
,
$A_n = 10$
,
$D_s = 0.9$
,
$\varLambda _S = 3.89$
,
$\varLambda _W = 675$
,
$Bi = 0.16$
,
$D_w = 2.94$
and
$\epsilon = 0.2$
. The total dimensionless freezing time is
$t_{\kern-1pt f} = 201$
. The results obtained using 12 001, 25 601 and 30 001 grid points are indistinguishable, demonstrating grid independence.
A.2. Validation using the configuration of Zadražil et al. (2006)
We further validate our freezing model by simulating a configuration previously investigated by Zadražil et al. (Reference Zadražil, Stepanek and Matar2006). Figure 17 shows the temporal evolution of the droplet shape,
$h$
(solid line), and the freezing front,
$s$
(dash-dotted line), for a droplet on a cold substrate, corresponding to the configuration studied by Zadražil et al. (Reference Zadražil, Stepanek and Matar2006). The droplet evolution observed here is consistent with the behaviour reported in figure 15 of their manuscript. The remaining dimensionless parameters are
$\epsilon = 0.2$
,
${\textit{Ste}} = 0.04$
,
$T_v = 0.5$
,
$A_n = 8$
,
$D_s = \varLambda _S = Bi = Bo = 1$
and
$D_w = 0$
. The dimensionless total freezing time obtained in figure 17 is
$t_{\kern-1pt f} = 14$
. To reproduce a droplet shape consistent with the simulations of Zadražil et al. (Reference Zadražil, Stepanek and Matar2006), we adopt the same initial droplet profile at the onset of freezing. In their study, the droplet undergoes spreading and imbibition prior to freezing; accordingly, we use the droplet shape obtained after these processes have completed as the initial condition for our simulations. This initial profile is taken from figure 14 of Zadražil et al. (Reference Zadražil, Stepanek and Matar2006). The total dimensionless time reported by Zadražil et al. (Reference Zadražil, Stepanek and Matar2006) is
$t_{\kern-1pt f} = 16$
, of which
$t = 2.08$
corresponds to spreading and imbibition, with the remaining duration associated with freezing (
$t_{\kern-1pt f} = 13.92$
). Consequently, the freezing time alone reported by Zadražil et al. (Reference Zadražil, Stepanek and Matar2006) agrees closely with the total freezing time
$t_{\kern-1pt f} = 14$
obtained in our simulations. The only differences between our simulations and those of Zadražil et al. (Reference Zadražil, Stepanek and Matar2006) are the thickness of the precursor layer, which is
$\beta = 0.01$
in the present study compared with
$\beta = 0.005$
of Zadražil et al. (Reference Zadražil, Stepanek and Matar2006), and the form of the penalty function employed in our model (2.52), which differs from that used in their work.
Evolution of the droplet shape,
$h$
(solid lines), and the freezing front,
$s$
(dash-dotted lines), on a substrate inclined at
$\alpha = 75^\circ$
for (a)
$\theta _{\textit{eq}} = 12^\circ$
and (b)
$\theta _{\textit{eq}} = 17.5^\circ$
. The normalised time,
$t/t_{\kern-1pt f}$
, corresponding to each row is indicated in panel (a). All other parameters are provided in table 2. The corresponding total freezing times are
$t_{\kern-1pt f} = 44$
and
$t_{\kern-1pt f} = 201$
for panels (a) and (b), respectively.
A.3. Grid convergence test
We perform a grid independence study for an equilibrium contact angle
$\theta _{\textit{eq}} = 17.5^\circ$
and an inclination
$\alpha = 75^\circ$
(near the extreme) using 12 001, 25 601 and 30 001 grid points. The remaining parameters (‘base’ parameters) used in the simulations are
$ \textit{Bo} = 0.4$
,
${\textit{Ste}} = 1.49 \times 10^{-3}$
,
$A_n = 10$
,
$D_s = 0.9$
,
$\varLambda _S = 3.89$
,
$\varLambda _W = 675$
,
$Bi = 0.16$
,
$D_w = 2.94$
and
$\epsilon = 0.2$
, as listed in table 2. The total freezing time (
$t_{\kern-1pt f}$
) remains unchanged across all cases, with
$t_{\kern-1pt f} = 201$
. We also note here that this solver has been extensively validated and employed in our previous studies (Kavuri et al. Reference Kavuri, Karapetsas, Sharma and Sahu2023, Reference Kavuri, Karapetsas, Sharma and Sahu2025).
Appendix B. Droplet shape evolution at
$\alpha = 75^\circ$
for different
$\theta _{\textit{eq}}$
Here, in figure 19, we compare the evolution of the droplet shape
$h$
(solid line) and the freezing front
$s$
(dash-dotted line) for droplets with
$\theta _{\textit{eq}} = 12^\circ$
and
$\theta _{\textit{eq}} = 17.5^\circ$
on a substrate inclined at
$\alpha = 75^\circ$
. The equilibrium contact angle
$\theta _{\textit{eq}} = 12^\circ$
corresponds to experimentally reported values for water on ice near the melting temperature. It can be seen that, for
$\theta _{\textit{eq}} = 12^\circ$
, the droplet becomes significantly more elongated along the inclined substrate, with pronounced swelling on the advancing side and tapering on the receding side. These shapes differ markedly from those obtained for
$\theta _{\textit{eq}} = 17.5^\circ$
, highlighting the strong influence of small equilibrium contact angles on the droplet morphology during freezing.
