2.1 Navier–Stokes equations

The liquid, assumed incompressible of density ( $\unicode[STIX]{x1D70C}_{l}$ ), surface tension ( $\unicode[STIX]{x1D6FE}_{l}$ ) and viscosity ( $\unicode[STIX]{x1D707}_{l}$ ), evolves following the mass and momentum conservation equations, which are given by the Navier–Stokes equations as follows:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{l}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{v}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}_{l}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+\unicode[STIX]{x1D70C}_{l}g\sin \unicode[STIX]{x1D6FC}\boldsymbol{i}-\unicode[STIX]{x1D70C}_{l}g\cos \unicode[STIX]{x1D6FC}\boldsymbol{k}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{v}$ is the velocity, $g$ is the gravitational acceleration; for generality the substrate may be tilted by an angle $\unicode[STIX]{x1D6FC}$ .

In addition to the fluid flow equations, the energy conservation equation has to be solved through both the solid substrate and the droplet. The phase change is accounted for by the coupling at the (water) liquid/ice interface which makes it possible to describe the freezing front evolution.

2.2 Energy conservation equation

The heat transfer equations within the droplet for both liquid and solid phases, as depicted in figure 1, can be written as

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{J}C_{PJ}\left(\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})T\right)=k_{J}\unicode[STIX]{x1D6FB}^{2}T\quad J=l,s, & & \displaystyle\end{eqnarray}$$

where the subscripts $l$ , $s$ stand for liquid and solid phases, respectively.

Under the long range approximation, these equations result in the Laplace equation, where thermal conduction is the main driving heat transfer mechanism. Regarding the heat transfer under phase change, this can be further justified by the fact that, based on the characteristic temperature difference $\unicode[STIX]{x0394}T$ , for example between the droplet and the substrate, the number $\unicode[STIX]{x0394}TC_{Pl}/L_{f}\ll 1$ , where $C_{Pl}$ , $L_{f}$ are the specific and latent heat of fusion, respectively. In other words, the transient heat conduction is negligible compared to the freezing front advance, as pointed out in Ajaev & Davis (Reference Ajaev and Davis2003). Finally, for both liquid and ice, and neglecting the transient contribution, the energy equation can be simplified as follows: for the solid (ice) phase,

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}T_{s}=0, & & \displaystyle\end{eqnarray}$$

and for the liquid phase,

(2.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{l}C_{Pl}(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})T=k_{l}\unicode[STIX]{x1D6FB}^{2}T. & & \displaystyle\end{eqnarray}$$

The advection term kept in the liquid heat transfer equation will be detailed later on when accounting for the effect of the density contrast between the liquid and solid (ice) phases.

The coupling between (2.1), (2.2) and (2.4), (2.5) will be performed under the lubrication approximation. Hence, the momentum equations (2.2) can be simplified assuming that the inertial effect is negligible, and can be rewritten by taking $\boldsymbol{v}=(\boldsymbol{u},w)$ , with the velocity horizontal components $\boldsymbol{u}=(u,v)$ on the substrate plane and a perpendicular component $w$ . In the case of a droplet under the lubrication approximation, the contact angle should be relatively low ( ${<}65^{\circ }$ ); therefore, our model is mainly valid for the hydrophilic configuration, even though an extension proposed in the present work makes it possible to handle droplet aspect ratios ( $H_{0}/R_{0}$ ) close to 1 ( ${\approx}90^{\circ }$ ). Under these hypotheses, the lubrication approximation, also known as the long range or thin film approximation, can be applied (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Alleborn & Raszillier Reference Alleborn and Raszillier2004) and the Navier–Stokes equations (2.1), (2.2) can be simplified as follows:

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\boldsymbol{u}+\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\tilde{\unicode[STIX]{x1D735}}p+\unicode[STIX]{x1D707}_{l}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{u}}{\unicode[STIX]{x2202}z^{2}}+\unicode[STIX]{x1D70C}_{l}g\sin \unicode[STIX]{x1D6FC}\boldsymbol{i}, & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}-g\cos \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D735}}\bullet =(\unicode[STIX]{x2202}\bullet /\unicode[STIX]{x2202}x,\unicode[STIX]{x2202}\bullet /\unicode[STIX]{x2202}y)$ stands for the planar gradient.

Integrating the pressure (2.8) with respect to $z$ leads to $p=-g\cos \unicode[STIX]{x1D6FC}z+\text{const.}$ and, on the other hand, applying the Young–Laplace boundary condition at the liquid–air interface and accounting for the disjoining pressure $\unicode[STIX]{x1D6F1}$ , one can express the pressure as follows:

(2.9) $$\begin{eqnarray}\displaystyle p(h+s)=-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}+P_{0}-\unicode[STIX]{x1D6F1}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ corresponds to the curvature at the droplet surface, $P_{0}$ is the atmospheric pressure, and $\unicode[STIX]{x1D6F1}$ corresponds to the disjoining pressure. The disjoining pressure is present owing to the precursor film, therefore (i) avoiding the singularity at the triple line and (ii) maintaining the droplet at equilibrium on the precursor film with a well-defined contact angle. At the molecular level, this disjoining pressure arises from the molecular interaction between the two interfaces, namely solid–liquid and liquid–gas. Finally, one can find the expression for the pressure within the droplet as

(2.10) $$\begin{eqnarray}p=-g\cos \unicode[STIX]{x1D6FC}(h+s-z)-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}-\unicode[STIX]{x1D72B}+P_{0}.\end{eqnarray}$$

By inserting the pressure equation in (2.7) and integrating twice with respect to $z$ between $[s,h+s]$ , we retrieve the quadratic velocity profile:

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D707}\boldsymbol{u}=-\tilde{\unicode[STIX]{x1D735}}P+\unicode[STIX]{x1D70C}_{l}g\sin \unicode[STIX]{x1D6FC}\boldsymbol{i}(z-s)\left[\frac{z-s}{2}-h\right],\end{eqnarray}$$

where $P=-g\cos \unicode[STIX]{x1D6FD}(h+s)-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}-\unicode[STIX]{x1D6F1}$ . This equation is derived under the no-slip condition at the liquid–ice interface $z=s$ , $\boldsymbol{u}=0$ , and the vanishing shear stress at $z=h+s$ , $\unicode[STIX]{x1D707}_{l}\unicode[STIX]{x2202}\boldsymbol{u}/\unicode[STIX]{x2202}z=0$ . Finally, one can deduce the average velocity over the liquid thickness defined as

(2.12) $$\begin{eqnarray}\boldsymbol{U}=\frac{1}{h}\int _{s}^{h+s}\boldsymbol{u}\,\text{d}z=-[-\tilde{\unicode[STIX]{x1D735}}P+\unicode[STIX]{x1D70C}_{l}g\sin \unicode[STIX]{x1D6FC}\boldsymbol{i}]\frac{h^{2}}{3\unicode[STIX]{x1D707}}.\end{eqnarray}$$

Knowing the velocity field in the liquid phase, we can now turn to the continuity equation of the freezing droplet. Integrating equation (2.6) with respect to $z$ between $[s,s+h]$ , we find

(2.13) $$\begin{eqnarray}\int _{s}^{s+h}\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}z+w|_{z=h+s}-w|_{z=s}=0.\end{eqnarray}$$

Using the Leibniz integral rule,

(2.14) $$\begin{eqnarray}\int _{s}^{s+h}\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}z=\tilde{\unicode[STIX]{x1D735}}.\int _{s}^{s+h}\boldsymbol{u}\,\text{d}z-[\boldsymbol{u}\tilde{\unicode[STIX]{x1D735}}(h+s)|_{h+s}-\boldsymbol{u}\tilde{\unicode[STIX]{x1D735}}s|_{s}],\end{eqnarray}$$

and using the definition of $\boldsymbol{U}$ from (2.12) and the no-slip condition at the liquid–ice interface $\boldsymbol{u}|_{s}\equiv 0$ , equation (2.14) can be simplified as

(2.15) $$\begin{eqnarray}\int _{s}^{s+h}\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}z=\tilde{\unicode[STIX]{x1D735}}.(h\boldsymbol{U})-\boldsymbol{u}\tilde{\unicode[STIX]{x1D735}}(h+s)|_{h+s},\end{eqnarray}$$

and finally noting that the kinematic boundary condition at the top of the droplet can be expressed as

(2.16) $$\begin{eqnarray}w|_{z=h+s}=\frac{\text{d}}{\text{d}t}(h+s)=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h+s)+\boldsymbol{u}\tilde{\unicode[STIX]{x1D735}}(h+s),\end{eqnarray}$$

then by combining (2.13) and (2.15) into (2.16), the continuity equation can be rewritten as follows:

(2.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(h+s)+\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }(h\boldsymbol{U})=w|_{z=s}.\end{eqnarray}$$

The right-hand side in (2.17) to be determined is the freezing front vertical speed $w|_{z=s}$ , which will be expressed by exploiting the mass and energy conservation at the interface between ice and liquid, as detailed below.

2.3 Mass conservation at the freezing front

The conservation of mass between the liquid and solid phases leads to

(2.18) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{l}(\boldsymbol{v}_{l}-\boldsymbol{v}_{i})\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D70C}_{s}(\boldsymbol{v}_{s}-\boldsymbol{v}_{i})\boldsymbol{\cdot }\boldsymbol{n},\end{eqnarray}$$

where the velocity of the solid (ice) phase is assumed to be zero, $\boldsymbol{v}_{s}\equiv 0$ . The normal unit vector at the interface $\boldsymbol{n}=(-\unicode[STIX]{x2202}_{x}s,-\unicode[STIX]{x2202}_{y}s,1)$ , and the interface velocity is given by $\boldsymbol{v}_{i}=(0,0,\unicode[STIX]{x2202}_{t}s)$ , hence (2.18) can be rewritten as

(2.19) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{l}(-u\unicode[STIX]{x2202}_{x}s-v\unicode[STIX]{x2202}_{x}s+w-\unicode[STIX]{x2202}_{t}s)=-\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x2202}_{t}s,\end{eqnarray}$$

which can be further simplified by using the no-slip condition at the interface $u\equiv v\equiv 0$ :

(2.20) $$\begin{eqnarray}w|_{z=s}=\left(1-\frac{\unicode[STIX]{x1D70C}_{s}}{\unicode[STIX]{x1D70C}_{l}}\right)\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

We can see that when $\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D70C}_{l}<1$ , there would be an expansion, and the induced velocity to the liquid phase due to expansion is positive.

Finally, inserting (2.20) in (2.17), one obtains the lubrication equation accounting for the phase change as follows:

(2.21) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }(h\boldsymbol{U})=-\frac{\unicode[STIX]{x1D70C}_{s}}{\unicode[STIX]{x1D70C}_{l}}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

The source term on the right-hand side describes the conversion between ice and liquid phases during the freezing process: an increase of the ice front leads to a decrease of the liquid, and that conversion rate is controlled by the ratio of the ice and liquid density.

2.4 Energy conservation at the freezing front

In addition to the mass conservation, energy conservation at the interface has to be applied to close the problem. At the interface, the energy balance through conduction and phase change is written as

(2.22) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{l}H_{l}(\boldsymbol{v}_{l}-\boldsymbol{v}_{i})\boldsymbol{\cdot }\boldsymbol{n}-k_{l}\unicode[STIX]{x1D735}T\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D70C}_{s}H_{s}(\boldsymbol{v}_{s}-\boldsymbol{v}_{i})\boldsymbol{\cdot }\boldsymbol{n}-k_{s}\unicode[STIX]{x1D735}T\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

Accounting for the properties of both ice and solid, the enthalpy jump at the interface corresponds to $\unicode[STIX]{x0394}H=H_{l}-H_{s}$ , with the enthalpies $H_{l}=C_{Pl}T_{m}+L_{f}$ and $H_{s}=C_{Ps}T_{m}$ of the liquid and solid, respectively, with phase change occurring at a definite temperature ( $T_{m}$ ). Assuming the heat transfer is dominant in the $z$ -direction with the following condition satisfied, $\unicode[STIX]{x2202}_{x}T\unicode[STIX]{x2202}s_{x}\propto \unicode[STIX]{x2202}_{y}T\unicode[STIX]{x2202}s_{y}\ll \unicode[STIX]{x2202}_{z}T$ , equation (2.22) can be expressed as

(2.23) $$\begin{eqnarray}[\unicode[STIX]{x1D70C}_{s}L_{f}+(\unicode[STIX]{x1D70C}_{l}C_{Pl}-\unicode[STIX]{x1D70C}_{s}C_{Ps})T_{m}]\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}=k_{s}\frac{\unicode[STIX]{x2202}T_{s}}{\unicode[STIX]{x2202}z}-k_{l}\frac{\unicode[STIX]{x2202}T_{l}}{\unicode[STIX]{x2202}z}.\end{eqnarray}$$

Equation (2.23) requires solving the heat transfer equation for both the liquid and ice phase. As previously indicated in § 2.2, these equations reduce to the Laplace equation, quasi-steady approximation, which will be solved in the $z$ -direction. The quasi-static temperature distribution provides a reasonable approximation as long as $C_{Pl}\unicode[STIX]{x0394}T<L_{f}$ .

2.5 Equilibrium temperature at freezing front: Gibbs–Thomson relation

Interfacial surface tension and the curvature affect the equilibrium temperature at the interface between the ice and liquid following the Gibbs–Thomson relation with the correction below for the melting temperature:

(2.24) $$\begin{eqnarray}T_{eq}=T_{m}\left(1-\frac{\unicode[STIX]{x1D6FE}_{SL}\unicode[STIX]{x1D705}_{s}}{\unicode[STIX]{x1D70C}_{s}L_{f}}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{SL}$ is the interfacial energy between ice–liquid, $\unicode[STIX]{x1D705}_{s}$ the curvature at the freezing front $s$ . Therefore the boundary condition for the temperature is expressed as

(2.25) $$\begin{eqnarray}T_{l}=T_{s}=T_{eq}\quad \text{at }z=s.\end{eqnarray}$$

It is worth noting that in (2.23), the prefactor on the left-hand side term has been approximated as $(\unicode[STIX]{x1D70C}_{l}C_{pl}-\unicode[STIX]{x1D70C}_{s}C_{ps})T_{m}(1-\unicode[STIX]{x1D6E4}\unicode[STIX]{x1D705})\approx (\unicode[STIX]{x1D70C}_{l}C_{pl}-\unicode[STIX]{x1D70C}_{s}C_{ps})T_{m}$ , since $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6FE}_{SL}/\unicode[STIX]{x1D70C}_{S}L$ and $\unicode[STIX]{x1D70C}_{l}C_{pl}-\unicode[STIX]{x1D70C}_{s}C_{ps}$ are usually small for water.

Knowing the melting temperature at the liquid/solid interface, the heat flux can be evaluated considering the thermal resistance of the substrate, $d_{W}/k_{W}$ , and the solid (ice), $s/k_{s}$ , along with the thermal contact resistance, $r_{c}$ , assumed between the liquid and the substrate figure 2. This contact resistance imposes a flux of $q=(T_{W}-T)/r_{c}$ at $z=0$ , on the substrate. From these considerations, one can deduce the flux at the liquid–ice interface, $z=s$ , as follows:

(2.26) $$\begin{eqnarray}q_{s}=-k_{s}\frac{\unicode[STIX]{x2202}T_{s}}{\unicode[STIX]{x2202}z}=-k_{s}\frac{T_{eq}-T_{W}}{s+d_{W}k_{s}/k_{W}+k_{s}r_{c}},\end{eqnarray}$$

where $T_{W}$ is the temperature imposed at the bottom of the substrate located at $z=-d_{W}$ .

Figure 2.

Schematic of the heat transfer model.

In order to find the heat flux from the liquid $q_{l}=-k_{l}\unicode[STIX]{x2202}_{z}T$ , we have to solve the heat transfer equation in the liquid phase:

(2.27) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{l}C_{Pl}w_{s}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}=k_{l}\frac{\unicode[STIX]{x2202}^{2}T}{\unicode[STIX]{x2202}z^{2}}.\end{eqnarray}$$

Under the quasi-static approximation, neglecting $\unicode[STIX]{x1D70C}_{l}C_{Pl}\unicode[STIX]{x2202}_{t}T$ as previously explained, and taking the value for the moving interface given by (2.20), $w_{s}=(1-\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D70C}_{l}){\dot{s}}$ , we can deduce the equation to be verified by the temperature field within the liquid phase as

(2.28) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}T}{\unicode[STIX]{x2202}z^{2}}-\frac{\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{s}}{\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FC}_{l}}{\dot{s}}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}=0,\end{eqnarray}$$

where the liquid thermal diffusivity is $\unicode[STIX]{x1D6FC}_{l}=k_{l}/\unicode[STIX]{x1D70C}C_{Pl}$ . This equation is slightly different from Laplace’s equation due to the contribution of the density difference. A similar contribution accounting for the density contrast between the solid and liquid phases can be found in Carslaw & Jaeger (Reference Carslaw and Jaeger1959) regarding the change of volume of the Stefan problem.

In addition to (2.25), the second boundary condition, neglecting the evaporation, is the heat transfer at the liquid interface between the droplet and the surrounding air, which is assumed to follow Newton’s law of cooling as follows:

(2.29) $$\begin{eqnarray}-k_{l}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}=h_{Tc}(T-\tilde{T}_{a}).\end{eqnarray}$$

$\tilde{T}_{a}$ is the surrounding temperature at the top of the static droplet, and $h_{Tc}$ is the heat transfer coefficient, which may account for the cooling of the initial (relatively hot) deposited droplet on the substrate before freezing. Since before depositing the droplet there is a thermal equilibrium between the substrate at $T_{W}$ and the atmosphere at $T_{a}$ , the temperature of the surroundings relevant for heat transfer between the droplet and the air that we consider is the one situated at a distance from the substrate corresponding to the height ( $H_{0}$ ) of the droplet, which we denote by $\tilde{T}_{a}$ . It is worth noting that this configuration is different from the convection resulting from plunging a solid sphere into a fluid at a given temperature; the deposited droplet case is different due to the presence of the substrate on which the temperature field depends. The difficulty in estimating the correct value of $h_{Tc}$ is well known and discussed in the literature. The correlations provided in the literature, for instance in the case of a sphere, do not apply to our configuration of a deposited droplet on a cold substrate. In order to compare our results with the experiment, one numerically determines $h_{Tc}$ along with $\tilde{T}_{a}$ (see appendix A for details). It is worth noting that the determination of this coefficient may affect the usability of the present model.

After solving (2.28), via the boundary conditions (2.25) and (2.29), one finds that the liquid heat flux in (2.23) is expressed as

(2.30) $$\begin{eqnarray}q_{l}=-k_{l}\frac{\unicode[STIX]{x2202}T_{l}}{\unicode[STIX]{x2202}z}=\frac{-\unicode[STIX]{x1D713}k_{l}h_{T}(T_{a}-T_{eq})}{h_{T}(\exp (\unicode[STIX]{x1D713}h)-1)+\unicode[STIX]{x1D713}k_{l}\exp (\unicode[STIX]{x1D713}h)},\end{eqnarray}$$

where we set for simplicity $\unicode[STIX]{x1D713}=((\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{s})/\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FC}_{l}){\dot{s}}$ , with ${\dot{s}}=\unicode[STIX]{x2202}s/\unicode[STIX]{x2202}t$ .

Combining (2.26) and (2.30) into (2.23) leads to

(2.31) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\tilde{L}_{f}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}=k_{s}\frac{T_{eq}-T_{W}}{s+d_{W}k_{s}/k_{W}+k_{s}r_{c}}-\frac{\unicode[STIX]{x1D713}k_{l}h_{T}(T_{a}-T_{eq})}{h_{T}(\exp (\unicode[STIX]{x1D713}h)-1)+\unicode[STIX]{x1D713}k_{l}\exp (\unicode[STIX]{x1D713}h)},\end{eqnarray}$$

where we denote the effective latent heat $\tilde{L}_{f}=L_{f}+T_{m}(\unicode[STIX]{x1D70C}_{l}C_{pl}-\unicode[STIX]{x1D70C}_{s}C_{ps})/\unicode[STIX]{x1D70C}_{s}$ . The solidification models in the literature do not account for this effect, which should not be neglected since it corrects the latent heat by approximately a factor of 3. Furthermore, that correction could become important when investigating the droplet freezing time, as sought in the present work. Besides, we can simplify (2.31), noting that $\unicode[STIX]{x1D713}h\ll 1$ , i.e. the contribution to the advection due the density contrast is low.

2.6 Dimensionless general equations for droplet freezing

In order to apply the model to the investigation of water droplet freezing, the governing equations will be expressed in cylindrical coordinates $(r,z)$ and in dimensionless form for convenience. Therefore from (2.21), we deduce

(2.32) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(rh\boldsymbol{U})=-\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}, & \displaystyle\end{eqnarray}$$

(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}=-h^{2}B_{p}\frac{\unicode[STIX]{x2202}\tilde{P}}{\unicode[STIX]{x2202}r}, & \displaystyle\end{eqnarray}$$

where $B_{p}=\unicode[STIX]{x1D6FE}{H_{0}}^{3}t_{c}/3\unicode[STIX]{x1D707}{R_{0}}^{4}$ and the density ratio $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D70C}_{l}$ . The coordinate $r$ is scaled by the droplet initial radius $R_{0}$ , and the heights $h$ and $s$ are scaled by the droplet initial height $H_{0}$ , while the temperature is scaled by $\unicode[STIX]{x0394}T=T_{m}-T_{W}$ , with $T_{m}$ the melting temperature of water. It is worth noting that, even if the basal radius of the droplet is constant, the droplet–air interface is subject to change during the freezing process; and the description of the droplet evolution equation, equation (2.32), is essential for the numerical model to track the droplet–air interface from the onset of freezing to the tip formation. Regarding the modelling of droplet freezing, two time scales are involved: the first one ( $t_{S}$ ) is related to droplet deformation while the freezing is concerned with the second one ( $t_{F}$ ). These two time scales are given by

(2.34a,b ) $$\begin{eqnarray}t_{S}=\frac{3\unicode[STIX]{x1D707}R_{0}^{4}}{\unicode[STIX]{x1D6FE}H_{0}^{3}}\quad \text{and}\quad t_{F}=\frac{\unicode[STIX]{x1D70C}_{s}\tilde{L}H_{0}^{2}}{k_{s}\unicode[STIX]{x0394}T}.\end{eqnarray}$$

The two characteristic times are relevant since the modelling is concerned with the freezing of a (deformable) droplet, as opposed to a static object. Only $t_{F}$ would have been relevant if the freezing had to be performed on a static configuration such as in the case of the Stefan problem. Therefore, we introduce $t_{c}$ , which combines the two time scales for the modelling of the freezing droplet and defines it as the product between the two characteristic times as

(2.35) $$\begin{eqnarray}t_{c}^{2}=\frac{3\unicode[STIX]{x1D707}R_{0}^{4}}{\unicode[STIX]{x1D6FE}H_{0}^{3}}\frac{\unicode[STIX]{x1D70C}_{s}\tilde{L}H_{0}^{2}}{k_{s}\unicode[STIX]{x0394}T}=\frac{3\unicode[STIX]{x1D707}\unicode[STIX]{x1D70C}_{s}\tilde{L}R_{0}^{4}}{\unicode[STIX]{x1D6FE}k_{s}H_{0}\unicode[STIX]{x0394}T}.\end{eqnarray}$$

The time scale $t_{c}=(3\unicode[STIX]{x1D707}\unicode[STIX]{x1D70C}_{s}\tilde{L}_{f}{R_{0}}^{4}/(\unicode[STIX]{x1D6FE}k_{s}H_{0}\unicode[STIX]{x0394}T))^{1/2}$ combines both the visco-capillary time scale $t_{S}=3\unicode[STIX]{x1D707}{R_{0}}^{4}/\unicode[STIX]{x1D6FE}{H_{0}}^{3}$ and the time scale of the freezing process $t_{F}=\unicode[STIX]{x1D70C}_{s}\tilde{L}{H_{0}}^{2}/k_{s}\unicode[STIX]{x0394}T$ . This is the time scale that accounts for the physics of our problem and is the one adopted throughout our modelling for droplet freezing. It is worth noting that the constant $B_{p}$ , which overall represents the visco-capillary contribution during the freezing process, is quite dependent on our choice of the characteristic time scale. For example, by choosing $t_{c}=t_{S}$ , i.e. the time is measured in the unit of the spreading time scale, then $B_{p}=1$ , and one retrieves the classical form for modelling droplet spreading under the lubrication approximation (Schwartz & Eley Reference Schwartz and Eley1998).

In order to have a similar droplet shape to the experiment, the expression of the full mean curvature in the capillary and disjointing pressure contribution is adopted, which enables our model to extend the one-dimensional approximation beyond the leading-order asymptotics. By doing so, the equilibrium shape of the droplet can be retrieved as a spherical cap. The use of the full mean curvature in the one-dimensional analysis is a powerful tool to capture the nonlinear dynamics of interfacial flows, as demonstrated in Eggers (Reference Eggers1993), Eggers & Dupont (Reference Eggers and Dupont1994) and Tembely et al. (Reference Tembely, Vadillo, Mackley and Soucemarianadin2012). It is worth noting that, traditionally, the thin film model for modelling droplet spreading does not make use of the full mean curvature expression, as applied in the present work (2.36). For a partially wetted surface, a disjoining pressure accounting for the intermolecular force is added to the static pressure inside the droplet (O’Brien & Schwartz Reference O’Brien, Schwartz and Hubbard2002; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). Here the two-term disjoining pressure, corresponding to a Lennard–Jones interaction potential ( $n=9$ , $m=3$ ), in agreement with the Frumkin–Derjaguin model (Schwartz & Eley Reference Schwartz and Eley1998) relating the static contact angle to interfacial energies, is used for the substrate wettability. The disjoining pressure prevents the thin film from thinning below the prescribed thickness $h^{\ast }$ . Taking the full mean curvature yields the pressure contribution term as follows:

(2.36) $$\begin{eqnarray}\displaystyle \tilde{P}=-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\frac{r(\unicode[STIX]{x2202}_{r}h+\unicode[STIX]{x2202}_{r}s)}{[1+\unicode[STIX]{x1D700}^{2}(\unicode[STIX]{x2202}_{r}h+\unicode[STIX]{x2202}_{r}s)^{2}]^{1/2}}-\unicode[STIX]{x1D6F1}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D700}=H_{0}/R_{0}$ , and $\unicode[STIX]{x1D6F1}=A[(h^{\ast }/h)^{n}-(h^{\ast }/h)^{m}]$ with $A=(1-\cos \unicode[STIX]{x1D703}_{E})(n-1)(m-1){R_{0}}^{2}/h^{\ast }(n-m){H_{0}}^{2}$ .

The disjoining pressure serves to maintain the droplet at the equilibrium position while the freezing front advances. It has little effect on the tip singularity as long as $h^{\ast }$ is small. At isothermal conditions, a higher value of $h^{\ast }$ causes the droplet to spread more (Schwartz & Eley Reference Schwartz and Eley1998).

Finally, accounting for the properties of the liquid water, ice and substrate, the enthalpy jump at the interface, equation (2.31), leads to the following evolution equation for the freezing front in dimensionless form:

(2.37) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t}=Ste\left[\frac{1-GT}{s+W+k_{s}R_{c}/H_{0}}-\unicode[STIX]{x1D706}\frac{B_{i}(\unicode[STIX]{x1D717}_{a}+GT)}{B_{i}h+1+\unicode[STIX]{x1D713}h}\right], & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D717}_{a}=T_{a}-T_{m}/(T_{m}-T_{W})$ and $\unicode[STIX]{x1D706}=k_{l}/k_{s}$ , while the Stefan number is defined as $Ste=t_{c}k_{s}\unicode[STIX]{x0394}T/\unicode[STIX]{x1D70C}_{s}\tilde{L_{f}}{H_{0}}^{2}$ , and the Biot number $Bi=h_{Tc}H_{0}/k_{l}$ . A summary of the main dimensionless parameters of the model is given in table 1. It is worth noting that a modified latent heat $\tilde{L}_{f}=L_{f}+T_{m}(\unicode[STIX]{x1D70C}_{l}C_{Pl}-\unicode[STIX]{x1D70C}_{s}C_{Ps})/\unicode[STIX]{x1D70C}_{s}$ is used in $Ste$ , which is often not considered in the literature of solidification models. Neglecting the density and heat capacity contrast in $\tilde{L}_{f}$ results in underestimating its value by almost a factor of three. Consequently, this effect is essential in order to accurately predict the water droplet features and freezing time, as considered below. In addition to the ratio of thermal conduction resistances of the substrate and the droplet $W=d_{W}k_{s}/H_{0}k_{W}$ , equation (2.37) includes the effect of the velocity induced during phase change, $\unicode[STIX]{x1D713}=(\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{s}){\dot{s}}/\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FC}_{l}$ , as a result of the density contrast between the liquid and ice.

Table 1.

Main dimensionless numbers used in the model.

Finally, the Gibbs–Thomson effect, termed GT in (2.37), yields

(2.38) $$\begin{eqnarray}\displaystyle GT=\frac{(\unicode[STIX]{x1D6FE}_{LS}/\unicode[STIX]{x1D70C}_{s}\tilde{L})T_{m}H_{0}}{{R_{0}}^{2}\unicode[STIX]{x0394}T}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{r\unicode[STIX]{x2202}_{r}s}{[1+\unicode[STIX]{x1D700}^{2}(\unicode[STIX]{x2202}_{r}s)^{2}]^{1/2}}\right). & & \displaystyle\end{eqnarray}$$

The GT contribution mainly serves to ensure the convergence and stability of the numerical solution, while its actual effect on the freezing seems to be marginal.

The previous nonlinear PDEs, equations (31)–(34), are implemented and solved using the method of lines suitable for stiff initial value and nonlinear problems. The UMFPACK linear solver and a maximum backward differentiation formula (BDF) order of 5 are used since these methods have been shown to be efficient to solve such initial value stiff nonlinear partial differential equations (Sellier & Panda Reference Sellier and Panda2010). On the boundary of the domain, the following conditions are imposed: $\unicode[STIX]{x1D735}P\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{n}=0$ . In one dimension, the mesh line between 0 and 2.5 is meshed with 1920 elements, and the precursor film thickness is taken as $h^{\ast }=0.006$ , which is found to be a good compromise between numerical accuracy and computational cost, while a Heaviside function is used to prevent precursor film freezing. When the droplet is assumed to be a spherical cap, the equation for the initial profile in dimensionless form $h_{0}(r)$ is expressed as follows:

(2.39) $$\begin{eqnarray}h_{0}(r)=\sqrt{\frac{\tilde{R}}{H_{0}^{2}}-\left(\frac{R_{0}}{H_{0}}\right)^{2}r^{2}}-\frac{\tilde{R}-H_{0}}{H_{0}},\end{eqnarray}$$

where $\tilde{R}=(H_{0}^{2}+R_{0}^{2})/2H_{0}$ , and the initial contact angle is given by $\tan \unicode[STIX]{x1D703}_{0}=2H_{0}R_{0}/(R_{0}^{2}-H_{0}^{2})$ . Note that the spherical cap approximation is fully valid as long as the Bond number $Bo=\unicode[STIX]{x1D70C}gH_{0}^{2}/\unicode[STIX]{x1D6FE}<1$ , or the droplet initial height ( $H_{0}$ ) is less than the capillary length.

3.1 Droplet freezing modelling

We investigated water droplet freezing at equilibrium on a cold solid substrate by solving the coupled PDEs (2.32)–(2.38). The model predictions are compared with one of the most detailed experimental results (to the best of our knowledge) of water droplet freezing, that reported in Hu & Jin (Reference Hu. and Jin2010). A $500~\unicode[STIX]{x03BC}\text{m}$ diameter droplet is deposited on a brass substrate of 1.5 mm thickness maintained at $-2\,^{\circ }\text{C}$ by circulating nitrogen gas. The use of the molecular tagging thermometry techniques (TAG) enabled the authors to quantify the water droplet freezing phenomena, including the freezing front shape along with the volume expansion. The parameters used in the model are given in table 2.

Table 2.

Parameters of the simulations.

In figure 3, the initial and final shapes of the frozen droplet are depicted both numerically and experimentally. The prediction of our model captures the main features of droplet freezing, such as the volume expansion and the pointy tip formation, with the tip angle at approximately $128^{\circ }$ , close to the $130^{\circ }$ found by Marín et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014). However, this angle is subject to variability, which is rationalized below. It is worth noting that, initially, under our one-dimensional approximation the droplet is taken as a spherical cap, which is a good approximation due to the millimetre size of the droplet (less than the capillary length); by this choice, the governing equations (i.e. (2.32)–(2.38)) are also satisfied as a result of taking the full mean curvature (see (2.36)).

Figure 3.

Experimental results (symbols) and theoretical prediction (solid line) of the initial and final, cusp-like, shape of a freezing water droplet on a solid substrate. The experimental results are from Hu & Jin (Reference Hu. and Jin2010) who used the molecular tagging thermometry techniques (TAG) for water droplet freezing. In addition, the tip angle, of approximately $128^{\circ }$ , is found to be in agreement with the recent finding in Marín et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014).

The evolution of the concave freezing front with time during the freezing is also evidenced, as highlighted in figure 4. This concave ice front is one of the features of water droplet freezing on a hydrophilic surface and is well captured by our model. Importantly, the pointy ice tip is also well retrieved.

Figure 4.

Evolution of predicted transient profile of a freezing droplet on a cold substrate at different times $0.2T_{s}$ , $0.4T_{s}$ , $0.6T_{s}$ and $T_{s}$ , where $T_{s}$ is the solidification time of the droplet.

The tip singularity, which our model retrieves, has remained a challenge in modelling solidification of water droplets. Towards the end of the solidification time ( $T_{s}$ ), the interface shape changes and the cusp singularity appears once the entire droplet has turned to ice. This singularity is the result of the upward liquid expansion during the freezing process, which ultimately focuses on the tip of the droplet, leading to the pointy ice shape. In addition, it is worth noting that the unfrozen liquid part within the droplet keeps its spherical cap shape thanks to the surface tension along with the adopted full mean curvature expression.

One of the interesting parameters in droplet freezing is the solidification time, which could determine the type of ice formation, glaze or rime ice. We tested the predictive capability of our approach regarding the effect of substrate temperature on the freezing time following the same range as in Hu & Jin (Reference Hu. and Jin2010). Using the physical and geometrical parameters in table 2, the comparisons between the proposed model and the experiments are very satisfactory over the range of temperatures from $-2\,^{\circ }\text{C}$ to $-5\,^{\circ }\text{C}$ as shown in figure 5. The analysis indicates that the substrate temperature is the main parameter which controls the droplet freezing time. A correlation shows that the evolution of the freezing time versus the substrate temperature follows a trend in line with the finding in Hu & Jin (Reference Hu. and Jin2010). Interestingly, from (2.35), the choice of the freezing time $T_{s}\propto 1/\unicode[STIX]{x0394}T$ yields the freezing front evolution to be precisely an inverse function of the substrate temperature, in line with the experimental data ( $T_{s}=60.42\unicode[STIX]{x0394}T^{-0.94}$ ). It is worth noting that the accuracy in predicting the freezing time is highly related, as previously explained in § 2, to the fact that a modified latent heat ( $\tilde{L_{f}}$ ) is used in the model. In addition to the operating parameter $T_{W}$ , the Stefan number ( $Ste$ ), based on the physical properties of the droplet, is another influencing parameter which can control water droplet freezing. In fact, the rate of phase change during the solidification phenomenon can be represented by the Stefan number, i.e. the ratio of sensible to latent heats; and we obtain in figure 6 a linear correlation between the freezing time and $1/Ste$ .

Figure 5.

Effect of the substrate temperature on the freezing time.

Figure 6.

Numerical simulation of the freezing time with respect to the inverse Stefan number $1/Ste$ . The inset corresponds to the freezing front evolution at different Stefan numbers.

Figure 7.

Comparative simulations of the effect of (a) temperature on a freezing droplet of aspect ratio $H_{0}/R_{0}=0.5$ and (b) the impact of the interaction between the droplet and atmosphere, as observed through the Biot number ( $Bi$ ), on the freezing front shape at different stages of droplet freezing ( $H_{0}/R_{0}=0.71$ ). The reference case corresponds to $Bi_{0}=0.4$ . The temporal sequence, $T_{f}/20$ , is similar in each comparative case.

Finally, in contrast to the substrate temperature, which has little effect on the freezing droplet features except the freezing time (see figure 7 a), the interaction between the solidifying droplet and the surrounding medium is influenced by both the Biot number ( $Bi$ ) and the temperature of the surroundings. The typical effect of $Bi$ on the freezing front shape is illustrated in figure 7(b). The right-hand side corresponds to a Biot number 20 times higher than the one on the left-hand side. The higher the Biot number or the heat transfer rate between droplet and surroundings, the more concave upwards is the freezing front. At lower Biot number, the outer droplet surface can be considered as adiabatic and there is no heat flux entering through it, therefore the temperature gradients inside the droplet are small. However at high $Bi$ number, the heat is evacuated at a higher rate through the droplet–air interface, which becomes cooler than the centre of the droplet; the freezing front then follows a curved isotherm freezing line. In addition, the simulation shows that the global shape of the droplet can be very much affected by the droplet heat transfer to the surrounding atmosphere; the frozen droplet shape at lower $Bi$ shows an inflection while this feature seems absent from the case of higher $Bi$ number. Furthermore, due to heat exchange between the droplet and the surroundings, the freezing front is no longer perpendicular to the interface. This finding suggests that the tip formation universality (Marín et al. Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014) may not hold when accounting for the interaction with the surroundings, and this may explain the variability observed experimentally in the value of the frozen droplet tip angle in Marín et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014).

3.2 Rationalizing the tip angle variability

Finally, based on our model, we propose to address and rationalize the variability observed in the tip angle of the frozen liquid water droplet (Marín et al. Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014). We provide below the derivation of the tip angle evolution, from (2.32):

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}\approx -\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

which leads after integration, to $h=-\unicode[STIX]{x1D70C}s+\text{const.}$ Initially the droplet shape is $h_{0}(r)$ while $s=0$ , so there is no ice formation yet, therefore $h(r,t)=-\unicode[STIX]{x1D70C}s(r,t)+h_{0}(r)$ . On the other hand, since the interface angle is given by $\unicode[STIX]{x2202}_{r}[h+s]=-\!\tan \unicode[STIX]{x1D703}$ , one can deduce

(3.2) $$\begin{eqnarray}h_{0}^{\prime }+(1-\unicode[STIX]{x1D70C})s^{\prime }=-\!\tan \unicode[STIX]{x1D703}.\end{eqnarray}$$

In order to describe the evolution of the tip angle at the final stage, one makes use of the freezing front evolution equation (2.37), neglecting however the GT effect and the induced velocity due the freezing, deriving with respect to $r$ :

(3.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}=Ste\left[-\frac{s^{\prime }}{(s+\unicode[STIX]{x1D71B})^{2}}+\unicode[STIX]{x1D706}\frac{B_{i}^{2}\unicode[STIX]{x1D717}_{a}h^{\prime }}{(B_{i}h+1)^{2}}\right].\end{eqnarray}$$

For simplicity, we set $\unicode[STIX]{x1D71B}=d_{w}k_{s}/h_{0}k_{w}+k_{s}r_{c}/h_{0}$ .

Close to the final stage of freezing, $r\rightarrow 0$ , $h\rightarrow 0$ , $s\rightarrow 1/\unicode[STIX]{x1D70C}$ , since $h_{0}(r\rightarrow 0)=1$ and $h_{0}^{\prime }(r\rightarrow 0)=0$ :

(3.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}=Ste\left[-\frac{s^{\prime }}{(1/\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D71B})^{2}}+\unicode[STIX]{x1D706}B_{i}^{2}\unicode[STIX]{x1D717}_{a}h^{\prime }\right].\end{eqnarray}$$

Finally the equation satisfied by $s^{\prime }$ is expressed as follows:

(3.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}=-Ste\left[\frac{1}{(1/\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D71B})^{2}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}B_{i}^{2}\unicode[STIX]{x1D717}_{a}\right]s^{\prime }.\end{eqnarray}$$

Noticing that close to the tip, $(1-\unicode[STIX]{x1D70C})s^{\prime }=-\!\tan \unicode[STIX]{x1D703}$ , and by setting $\unicode[STIX]{x1D709}=\tan \unicode[STIX]{x1D703}$ , the following equation can be written:

(3.6) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}t}=-Ste\left[\frac{1}{(1/\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D71B})^{2}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}B_{i}^{2}\unicode[STIX]{x1D717}_{a}\right]\unicode[STIX]{x1D709},\end{eqnarray}$$

which can be readily solved as

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D709}=C\exp \left[-Ste\left(\frac{1}{(1/\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D71B})^{2}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}B_{i}^{2}\unicode[STIX]{x1D717}_{a}\right)t\right].\end{eqnarray}$$

Furthermore, close to the final stage $r\rightarrow 0$ , $t-ts\rightarrow 0$ , the tip angle can be expressed by $\unicode[STIX]{x1D6FC}_{tip}=\unicode[STIX]{x03C0}-2\unicode[STIX]{x1D709}$ , while $\unicode[STIX]{x1D709}=\tan \unicode[STIX]{x1D703}$ verifies the following equation:

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}_{0}\left[1+Ste\left(\frac{1}{(1/\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D71B})^{2}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}B_{i}^{2}\unicode[STIX]{x1D717}_{a}\right)(t_{s}-t)\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D709}_{0}$ corresponds to the tip angle at the final freezing time $T_{s}$ , while for simplicity we set $\unicode[STIX]{x1D71B}=d_{W}k_{s}/h_{0}k_{W}+k_{s}R_{c}/h_{0}$ .

Interestingly, equation (3.8) may rationalize as well the seemingly uncorrelated result of figure 2(b) in Marín et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014), where an increase of the contact angle or height-to-radius aspect ratio seems to induce a higher tip angle. In fact, by combining both the aspect ratio and the substrate temperature effects, the prefactor in (3.8) seems to predict qualitatively the correct change in the tip angle. In fact, based on the operating conditions, aspect ratio and substrate temperature in figure 2(b) from Marín et al. (Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014), the prefactor calculation leads to $\unicode[STIX]{x1D709}_{H/R=0.46,\unicode[STIX]{x0394}T=20.8C}=2.6>\unicode[STIX]{x1D709}_{H/R=1.16,\unicode[STIX]{x0394}T=44.1C}=1.1>\unicode[STIX]{x1D709}_{H/R=4.96,\unicode[STIX]{x0394}T=34.4C}=1.0$ , equivalently the tip angle $\unicode[STIX]{x1D6FC}_{H/R=0.46,\unicode[STIX]{x0394}T=20.8C}<\unicode[STIX]{x1D6FC}_{H/R=1.16,\unicode[STIX]{x0394}T=44.1C}<\unicode[STIX]{x1D6FC}_{H/R=1.16,\unicode[STIX]{x0394}T=44.1C}$ . As can be inferred by (3.8), the interplay between all the physical parameters may rationalize the reason for the variability observed experimentally. It is worth noting that our approach, based on a scaling argument, can only provide qualitative results regarding the tip angle variability.