2.1 Presentation of the problem

We consider a material existing in two phase states: liquid and solid. Then, an infinite flat solid substrate, filling the half-space $z<0$ , with a temperature $T_{s}$ below the melting temperature $T_{m}$ , is in contact with the solid layer of the material (for $0<z<h(t)$ ), that is growing in its melt ( $z>h(t)$ ). Here $h(t)$ is the position of the solidification front (see schematic figure 1). This configuration belongs to the large class of Stefan problems, a particular kind of boundary value problem where a phase boundary can move with time. It is named after Josef Stefan who first, in 1889, solved the original configuration in which a solidification front propagates between two phases (liquid and solid) of the same material, without substrate (Stefan Reference Stefan1891; Brillouin Reference Brillouin1930). In figure 1, this original Stefan problem amounts to replacing the substrate by the solid phase.

Figure 1.

Summary of the model hypotheses: a finite layer of solid lies between the semi-infinite melt ( $z>h(t)$ ) and the semi-infinite substrate ( $z<0$ ). The temperature of the whole melt is set constant at the melting point ( $T=T_{m}$ ), while the temperature of the substrate tends to $T_{S}$ when $z$ goes to  $-\infty$ . The temperature in the solid phases is given by a set of two heat equations, with a specific diffusion coefficient $D_{k}$ for each phase, coupled by the temperature and heat flux continuity at $z=0$ . At the solidification front ( $z=h(t)$ ), the Stefan condition imposes the downward thermal flux be equal to the latent heat liberated by the freezing.

Configurations belonging to Stefan problems, even restricted to the solidification–melting phase change, are abundant, with an impressively wide range of applications. We can cite as examples, the solidification of the earth, supposed to be molten at the origin, which triggered the first work of interest in this area, by Lamé & Clapeyron (Reference Lamé and Clapeyron1831). The formation of ice crystals, like snowflakes, that grow out from seeds in an environment of supersaturated water vapour. In this problem the ingredients of the original Stefan problem are not sufficient, Gibbs–Thomson equation and three-dimensional (3-D) effects need to be added (Langer Reference Langer1980). The modelling of the dynamics of sea ice on the surface of the polar oceans (Worster Reference Worster2000) or the shape of an icicle (Neufeld, Goldstein & Worster Reference Neufeld, Goldstein and Worster2010), where the liquid motions and the associated convective heat transfer have to be considered. Likewise for various lava flows that cool and gradually solidify until they come to rest (Griffiths Reference Griffiths2000). Cooling can occur from the surrounding atmosphere (or water) or from the underlying solid (Huppert Reference Huppert1989), and in ancient times, some of the lavas were hotter than today and were even capable of melting the underlying rock and shaping their own thermal erosion bed (Huppert Reference Huppert1986). In industry, Stefan problems are frequently studied in the context of melting or solidification of metals or metal alloys (Viskanta Reference Viskanta1988), where the properties of the solid (in particular, its mechanical and thermal properties) are functions of the kinetics of solidification, such as Czochralski crystal growth, used for example in the fabrication of semi-conductor wafers (Nishinaga Reference Nishinaga2014), laser welding (Cline & Anthony Reference Cline and Anthony1977; von Allmen & Blatter Reference von Allmen and Blatter2013), or the synthesis of nanoparticles from melting metal films (Font, Afkhami & Kondic Reference Font, Afkhami and Kondic2017). Finally, this class of problems is also a fantastic playground for mathematicians (Rubinstein Reference Rubinstein1971; Gupta Reference Gupta2003); but to the best of our knowledge, we are not aware of a system of equations modelling solid growth by a sudden contact between a liquid and a cold solid substrate.

In the following, we assume the melt stays at rest, at constant temperature everywhere ( $T_{m}$ ). We neglect the variation of heat capacity and of thermal conductivity with temperature. We also neglect the thermal expansion, and more generally the variation of density. In other words, the thermal parameters of the media are the following: the latent heat of solidification  $L$ , the heat capacity $C_{pk}$ , the thermal conductivity  $\unicode[STIX]{x1D706}_{k}$ , the heat diffusion coefficient $D_{k}=\unicode[STIX]{x1D706}_{k}/(\unicode[STIX]{x1D70C}_{k}C_{pk})$ and therefore the density  $\unicode[STIX]{x1D70C}_{k}$ , are taken as constant. The subscript $k$ stands for $l$ in the liquid phase, $i$ in the solid phase and $s$ for the substrate (see figure 1).

The choice of a one-dimensional model will be discussed in § 2.3 and is supported by the geometry of the impacted drop that will be considered later. We assume that the liquid is at the melting temperature  $T_{m}$ ; this approximation will also be discussed in § 2.3. It can be justified a priori by both, the small thickness of the liquid layer and the little energy needed to cool the water down to its melting temperature,  $C_{pl}(T_{d}-T_{m})\sim 4000\times 20\sim 8\times 10^{4}~\text{J}~\text{kg}^{-1}$ , where $T_{d}$ is the initial temperature of the liquid drop, compared to the latent heat for solidification ( $L\sim 3\times 10^{5}~\text{J}~\text{kg}^{-1}$ ).

Under these assumptions, the mathematical problem that needs to be solved is the following (see figure 1): a constant temperature $T_{m}$ in the liquid, two heat equations for the temperature field  $T(z,t)$

(2.1a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}=D_{s}\frac{\unicode[STIX]{x2202}^{2}T}{\unicode[STIX]{x2202}z^{2}}\quad \text{for}~z\leqslant 0;\quad \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}=D_{i}\frac{\unicode[STIX]{x2202}^{2}T}{\unicode[STIX]{x2202}z^{2}}\quad \text{for}~0\leqslant z\leqslant h(t),\end{eqnarray}$$

and four boundary conditions at the two interfaces

(2.2a,b ) $$\begin{eqnarray}T(0^{-},t)=T(0^{+},t);\quad \unicode[STIX]{x1D706}_{s}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}(0^{-},t)=\unicode[STIX]{x1D706}_{i}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}(0^{+},t),\end{eqnarray}$$

which impose both the continuity of the temperature and of the heat flux at the substrate/solid interface ( $z=0$ ), and

(2.3a,b ) $$\begin{eqnarray}T(h(t)^{-},t)=T_{m};\quad \unicode[STIX]{x1D706}_{i}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}(h(t)^{-},t)=\unicode[STIX]{x1D70C}_{i}L\frac{\text{d}h}{\text{d}t},\end{eqnarray}$$

which impose both the continuity of the temperature, and the law of motion of the solidification front ( $z=h(t)$ ). This energy conservation law under liquid–solid phase change is called the Stefan condition. Here, it dictates that the solidification front velocity is proportional to the rate at which latent heat can be transported in the solid phase. The solidification front is thus controlled by the diffusion in the solid and the substrate through the Stefan condition. Finally, we impose a constant temperature $T_{s}$ deep in the substrate

(2.4) $$\begin{eqnarray}\lim _{z\rightarrow -\infty }T(z,t)=T_{s},\end{eqnarray}$$

and we complement this set of equations by the initial conditions taken at $t=0$

(2.5a,b ) $$\begin{eqnarray}T(z,0)=T_{s}\quad \text{for}~z\leqslant 0\quad \text{and}\quad h(0)=0,\end{eqnarray}$$

indicating that at $t=0$ the liquid is suddenly put in contact with the substrate.

2.2 Solution of our unidimensional solidification model

Similarity analysis shows that this diffusive problem exhibits a self-similar structure, with the usual self-similar variables involved in diffusion problems, namely $z/\sqrt{D_{i}t}$ in the ice and $z/\sqrt{D_{s}t}$ in the substrate, even in the presence of the moving solidification front. In this case, the solidification front location also follows a square-root in time law, showing the diffusive property of the dynamics

(2.6) $$\begin{eqnarray}h(t)=\sqrt{D_{eff}t},\end{eqnarray}$$

where $D_{eff}$ is the effective diffusion coefficient that determines the growth of the solid layer. It is really the quantity of interest that we need to compute and thus relates to the different thermal properties of our problem.

Introducing the self-similar variable in the set of equations (2.1)–(2.5), we obtain the following solutions for the temperature field

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle T(z,t)=T_{0}+(T_{0}-T_{s})\cdot \text{Erf}\left(\frac{z}{2\sqrt{D_{s}t}}\right)\quad \text{for}~z\leqslant 0\quad \text{and} & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle T(z,t)=T_{0}+\frac{e_{s}}{e_{i}}(T_{0}-T_{s})\cdot \text{Erf}\left(\frac{z}{2\sqrt{D_{i}t}}\right)\quad \text{for}~0\leqslant z\leqslant h(t), & \displaystyle\end{eqnarray}$$

where $e_{s,i}=\sqrt{\unicode[STIX]{x1D706}_{s,i}\unicode[STIX]{x1D70C}_{s,i}C_{ps,i}}$ are the effusivities of the substrate and the solid, and $T_{0}$ the contact temperature at the solid–substrate interface (a constant in time in this self-similar framework). The effusivity of a material is the physical quantity that shows both its heat capacity and its ability to diffuse it. $T_{0}$ is an integration constant to be determined by the boundary conditions. It corresponds to the temperature at the solid–substrate interface and the self-similar behaviour indicates that this temperature is a constant in time. The error function, $\text{Erf}$ is defined here by

(2.9) $$\begin{eqnarray}\text{Erf}(x)=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{x}\text{e}^{-\unicode[STIX]{x1D709}^{2}}\,\text{d}\unicode[STIX]{x1D709}.\end{eqnarray}$$

Then, by imposing the Stefan condition at $z=h(t)$ (2.3), we obtain the following transcendental equation

(2.10) $$\begin{eqnarray}St=\frac{\sqrt{\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FD}}}{2}\text{e}^{\unicode[STIX]{x1D6FD}/4}\left(\frac{e_{i}}{e_{s}}+\text{Erf}\left(\frac{\sqrt{\unicode[STIX]{x1D6FD}}}{2}\right)\right)\end{eqnarray}$$

that links the Stefan number, $St$ , with the ratio of the diffusion coefficients,  $\unicode[STIX]{x1D6FD}$ , defined, respectively, by

(2.11a,b ) $$\begin{eqnarray}St=\frac{C_{pi}(T_{m}-T_{s})}{L}\quad \text{and}\quad \unicode[STIX]{x1D6FD}=\frac{D_{eff}}{D_{i}}.\end{eqnarray}$$

This implicit relation (2.10) has to be solved numerically and the values of $D_{eff}$ used further on are those obtained numerically. From the relation (2.10) it is easy to deduce the asymptotic behaviours for small and large Stefan numbers

(2.12a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}\sim \frac{4e_{s}^{2}}{\unicode[STIX]{x03C0}e_{i}^{2}}St^{2}\quad \text{for}~St\ll 1,\quad \text{and}\quad \unicode[STIX]{x1D6FD}\sim 4\ln (St)\quad \text{for}~St\gg 1.\end{eqnarray}$$

An interesting physical quantity to compute is the dimensionless solid–substrate interface temperature $\bar{T}_{0}$ yielding

(2.13) $$\begin{eqnarray}\overline{T}_{0}=\frac{T_{0}-T_{s}}{T_{m}-T_{s}}=\frac{1}{1+{\displaystyle \frac{e_{s}}{e_{i}}}\text{Erf}\left({\displaystyle \frac{\sqrt{\unicode[STIX]{x1D6FD}}}{2}}\right)}\end{eqnarray}$$

such that $\overline{T}_{0}$ varies between 0 and 1 with $\overline{T}_{0}=0$ for $T_{0}=T_{s}$ and $\overline{T}_{0}=1$ for $T_{0}=T_{m}$ (see figure 1). Then from the asymptotic relations between $\unicode[STIX]{x1D6FD}$ and $St$ (2.12), we find that $T_{0}\rightarrow T_{m}$ for $St\ll 1$ , while for $St\gg 1$ it converges to an intermediate temperature

(2.14) $$\begin{eqnarray}\lim _{St\rightarrow \infty }T_{0}=T_{s}+(T_{m}-T_{s})\frac{1}{1+{\displaystyle \frac{e_{s}}{e_{i}}}},\end{eqnarray}$$

which corresponds to the interface temperature when two infinite media are suddenly put in contact (de Ruiter et al. Reference de Ruiter, Soto and Varanasi2018).

Figure 2.

Results of the model for the different substrates (copper: red, steel: orange, marble: light blue) and ice (dark blue) that is for comparison with the Stefan problem with infinite ice and infinite liquid water. (a,b) Dimensionless temperature profiles $\overline{T}(z,t)=(T(z,t)-T_{s})/(T_{m}-T_{s})$ as functions of the variable $z/\sqrt{D_{eff}t}$ , obtained from (2.8) and (2.7) for two characteristic Stefan numbers $St=250$  (a) and $St=0.25$  (b). This rescaling accounts well for the two self-similar variables valid in each domain (ice and substrate). (c) Dimensionless contact temperature, obtained from (2.13). Regardless of the substrate, it goes from 1 for $St\rightarrow 0$ to the asymptotical values $(1+e_{s}/e_{i})^{-1}$ (indicated on the figure) when $St\rightarrow \infty$ , reminiscent of the contact temperature of two infinite bodies initially at different temperatures. (d) The effective diffusion coefficient of the solidification front normalized by the thermal diffusion coefficient of the ice,  $\unicode[STIX]{x1D6FD}=D_{eff}/D_{i}$ , as a function of the Stefan number. The differences between substrates are more important at low Stefan number, and the asymptotic regime at low Stefan number $\unicode[STIX]{x1D6FD}\propto St^{2}$ is indicated (dashed line). Here $\unicode[STIX]{x1D6FD}$ has been computed by solving numerically the implicit equation (2.10).

In order to understand the thermal fields in the substrate and the solid layer, we plot the vertical position normalized by the thickness of the growing ice layer ( $z/\sqrt{D_{eff}t}$ ) as a function of the normalized temperature profiles

(2.15) $$\begin{eqnarray}\overline{T}(z,t)=\frac{T(z,t)-T_{s}}{T_{m}-T_{s}}\end{eqnarray}$$

for two specific values of the Stefan number: a large value, $St=250\gg 1$ (figure 2 a) and a value of order one, $St=0.25\sim 1$ (figure 2 b); and for four pairs of melt/substrate. The melt is water as it is the most common in our experiments and the four substrates have been chosen: copper, steel and marble in order to have a large range of thermal conductivities, and ice in order to compare our result with the classical Stefan problem with no substrate. We observe that the temperature profiles evolve from the substrate temperature $T_{s}$ ( $\overline{T}(z\rightarrow -\infty ,t)=0$ ) up to the solidification temperature $T_{m}$ ( $\overline{T}(h(t),t)=1$ ). In the high Stefan numbers regime (figure 2 a), the temperature profiles exhibit the shape of the error function within both the solid and the substrate domains. The temperature gradient is discontinuous at the contact between the ice and the substrate ( $z=0$ ), because of the discontinuity of the thermal conductivities. For $St=0.25$ (figure 2 b), the temperature profile tends to be linear in the ice while it still exhibits an error function profile in the substrate.

Figures 2(c) and 2(d) present, respectively, the normalized temperature of the solid–substrate interface, $\overline{T}_{0}$ (2.13), and the normalized effective diffusion coefficient $\unicode[STIX]{x1D6FD}=D_{eff}/D_{i}$ (by solving (2.6)), both plotted as functions of the Stefan number, for the same four pairs of melt/substrate as before. Remember that increasing the Stefan number corresponds to a decrease of the substrate temperature.

To better understand the different regimes observed in figure 2, we can investigate the two asymptotic behaviours of our model: the solidification-dominated regime ( $St\ll 1$ ) and the solid-cooling dominated ( $St\gg 1$ ). Let us first consider the latter case of large Stefan numbers ( $St\gg 1$ ) that corresponds to $C_{pi}\unicode[STIX]{x0394}T\gg L$ , meaning that the latent heat $L$ released by the solidification at the solid–liquid interface $z=h(t)$ , is negligible compared to the heat energy released by cooling the solid ( $C_{pi}\unicode[STIX]{x0394}T$ ) from the melting temperature $T_{m}$ to the substrate temperature $T_{s}$ . In this case, the substrate is hardly warmed up by the solidification and does not influence the front propagation. As a consequence, the dynamics is dominated by the self-similar variation of the temperature profiles with time; they exhibit two $\text{Erf}$ functions that join at the solid–substrate interface (figure 2 a). If the diffusive coefficients in the solid and the substrate are the same, which here is the case when the substrate is ice and almost the case when the substrate is marble, then the $\text{Erf}$ shape of the temperature profile is symmetrical with respect to $z=0$ . In other cases (steel and copper), this symmetry of the error function is broken. In this limit $(St\gg 1)$ , the temperature of the solid–substrate interface is the asymptotic limit given by the (2.14) (figure 2 c). Finally, figure 2(d) confirms that in this limit ( $St\gg 1$ ), the front propagation dynamics is not influenced much by the substrate thermal conductivity as all the substrates tend to the classical Stefan problem (dark blue curve) with no substrate.

In the opposite limit ( $St\ll 1$ ), the latent heat released by the solidification is much larger than the heat released by cooling the solid down to  $T_{s}$ , indicating that the heat reaching the substrate comes mainly from the solid–liquid interface. The solid layer is inert in the sense that it only transfers the heat flux coming from the solidification to the substrate; the heat flux stays constant through the solid layer. This explains why the temperature field is linear in the solid (figure 2 c). The slope is selected by the front dynamics and evolves slowly with time, it can be estimated using (2.12) for $St\ll 1$ ( $\unicode[STIX]{x1D6FD}\sim St^{2}$ ) as follows

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{i}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}=\unicode[STIX]{x1D70C}_{i}L\frac{\text{d}h}{\text{d}t}=\frac{1}{2}\unicode[STIX]{x1D70C}_{i}L\sqrt{\frac{D_{eff}}{t}}\sim St\unicode[STIX]{x1D70C}_{i}L\sqrt{\frac{D_{i}}{t}}.\end{eqnarray}$$

This linear profile in the solid joins an $\text{Erf}$ profile in the substrate which goes deeper (relatively to $h(t)$ ) compared to the large Stefan regime. The freezing process is efficient enough to warm the substrate up. This explains why the interface temperature $T_{0}$ increases (figure 2 d) and tends toward $T_{m}$ ( $T_{0}\rightarrow T_{m}$ for $St\rightarrow 0$ (2.13)). Finally, the effective diffusion coefficient is smaller than in the large Stefan regime (figure 2 d). This can be understood easily since decreasing the Stefan number corresponds to warming the substrate up, so the ice grows slower. In this regime we observe that all the curves follow $D_{eff}\sim St^{2}D_{i}$ as shown on (2.12). What is very interesting is that $D_{eff}$ with metal substrates differs by several orders of magnitude from the one with no substrate (dark blue curve). Indeed, as the Stefan number decreases, the heat flux released by cooling the solid down to $T_{s}$ becomes negligible with respect to the latent heat. As a consequence, the solidification front propagation becomes more and more influenced by the substrate’s thermal parameters. In this case, the higher the substrate’s thermal conductivity, the faster the solidification front propagation (see (2.12)). This property demonstrates how important it is to take into account the heat transfer within the substrate.

2.3 Experimental comparison

With the aim of validating our model and its assumptions, we carry out a model experiment which consists of putting in contact a liquid with a substrate at a temperature below the liquid solidification temperature (see figure 3 a). The liquid is initially poured into a transparent PMMA tube of 3.4 cm diameter. We use two different liquids (water and hexadecane) and two different substrates (copper and steel). We also varied the initial temperature of the water (20 and  $0\,^{\circ }\text{C}$ ) and the thickness of the initial water layer (between 4 mm and a few centimetres). Finally, the propagation of the solidification front ( $h(t)$ ) is recorded using a camera, between one and ten frames per second depending on the experiment.

Figure 3.

(a) Schematic of the experimental set-up: a liquid-filled cylinder is set upon a cold substrate. The liquid starts freezing into a solid layer of thickness $h(t)$ , which is monitored using a camera. (b) Growth of the solid phase as a function of time. Here $h(t)$ is plotted in the case of water freezing on a copper plate at two different temperatures (blue and red squares), the liquid being initially at room temperature. The same experiments are repeated with water slightly warmer than $0\,^{\circ }\text{C}$ (triangles), showing no differences. The dynamics of each experiment is fitted by a square root function (full lines), in order to get its diffusion coefficient $D_{exp}$ . Dashed lines represents the dynamics predicted by the classical Stefan problem (without the substrate), for the two different temperatures. (c) Comparison of the experimental freezing dynamics and the theory: the diffusion coefficients $D_{exp}$ and $D_{eff}$ , respectively, experimental and theoretical, are compared for a wide range of parameters. The colour refers to the substrate cooling $\unicode[STIX]{x0394}T=T_{s}-T_{m}$ . Circles (respectively squares) represents water at room temperature freezing on copper (respectively steel). Stars refer to thin thicknesses of water freezing on copper, triangles to $0\,^{\circ }\text{C}$ water freezing on copper. Finally, pentagons are for hexadecane freezing on both metals. The $y=x$ guideline is shown as a dashed line.

Figure 3(b) presents the growth of the ice layer ( $h(t)$ ) with respect to time for water solidifying on a copper substrate. The blue and red squares are for water at $20\,^{\circ }\text{C}$ and two different temperatures of the substrate, respectively, $-38\,^{\circ }\text{C}$ and $-26\,^{\circ }\text{C}$ as indicated on the graph. As expected the freezing front propagates faster when the substrate is colder. On the same graph, the empty triangles show the front propagation with water initially at $0\,^{\circ }\text{C}$ . Their variations coincide with the preceding ones on the whole range of time, indicating that solidifying water at $20\,^{\circ }\text{C}$ or at $0\,^{\circ }\text{C}$ gives rise to the same solidification dynamics. This important observation justifies a posteriori the approximation to consider the liquid layer at $T_{m}$ in the model and to neglect the heat flux in the liquid. Superimposed on these four experimental curves are lines which are the best fit using a square root function of the form

(2.17) $$\begin{eqnarray}h(t)=\sqrt{D_{exp}t},\end{eqnarray}$$

with $D_{exp}$ an experimental diffusion coefficient that quantifies the solidification dynamics. We observe that, as shown in the model, the solidification front follows a diffusive law. Note that in the experiments shown here, the same diffusive law is followed up to 6 mm, suggesting that the 1-D approximation seems valid at least up to an aspect ratio of around one fifth. Finally, we plot on the same graph with dashed line the variation of $h(t)$ given by the classical Stefan problem associated to the two substrate temperatures. In these plots, the substrate is treated as if it were ice and the thermal parameters of the substrate are consequently not taken into account. It appears that, by promoting the heat transfer, the substrate highly increases the front’s propagation velocity.

By fitting the solidification front evolution by (2.17) on each experiment we obtain $D_{exp}$ , the experimental diffusion coefficients. On figure 3(c) we compare them to the effective diffusion coefficient given by our 1-D model $D_{eff}$ , for all of our experiments: water (0, $20\,^{\circ }\text{C}$ and different liquid thicknesses) and hexadecane, on copper and steel (see the different markers on the graph), and for various substrate temperatures between 0 and $-50\,^{\circ }\text{C}$ , indicated by the colour of the markers. All the data gather along the black dashed line that has a slope of 1, indicating that for each experiment $D_{exp}\simeq D_{eff}$ : our 1-D model provides therefore an excellent estimation of the dynamics of the upwards solidification when a liquid is placed on the top of a cold substrate.

3.1 Experimental set-up and qualitative description

The classical drop impact set-up consists of a syringe pump pushing a liquid through a capillary tube from which the drop falls. As the pumping is slow enough, the size of the drop is entirely controlled by the radius of the capillary tube. We used two different tubes of inner diameter $1600$ and $250~\unicode[STIX]{x03BC}\text{m}$ leading to two drop radii: $R=1.9$ and 1.2 mm, yielding, respectively, volumes $30~\unicode[STIX]{x03BC}\text{L}$ and $7~\unicode[STIX]{x03BC}\text{L}$ . The impact velocity $U_{0}$ is controlled by the height of fall  $H$ , which in our case ranges from 15 to 45 cm, so that $U_{0}$ ranges from 1.7 to $3~\text{m}~\text{s}^{-1}$ (following roughly $U_{0}=\sqrt{2gH}$ ).

Our substrates, made of blocks ( $100~\text{mm}\times 100~\text{mm}\times 30$  mm) of different materials (steel, copper and marble), are placed into a bowl and cooled down by pouring a certain amount of liquid nitrogen. The minimal temperature reached in this work is around $-80\,^{\circ }\text{C}$ . Due to the substrate heat capacity and the bowl’s thermal isolation, it takes several hours for the system to warm up to room temperature. The change in the substrate temperature is thus much less than $1\,^{\circ }\text{C}$ during the time of an experiment which is roughly 1 s.

In order to minimize frost formation, the whole system is placed inside a regulated atmosphere chamber which allows us to drastically reduce the humidity (less than 1 % humidity inside the chamber). The substrate temperature $T_{s}$ is measured before each experiment using a surface thermometer. The dynamics of the impact is studied using a high-speed camera, and the height profile of the frozen drop is extracted with a polychromatic confocal sensor (CCS OptimaPlus from STIL Optics) moving along a translation platform.

At room temperature, as a water drop impacts a solid substrate, it spreads, reaches a maximal radius (Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014) and immediately starts to retract back to an equilibrium radius (Bartolo et al. Reference Bartolo, Josserand and Bonn2005; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016). Figure 4(a–h) shows a sequence of snapshots of a water drop impacting a sub-zero substrate. Here again, the drop spreads rapidly and reaches a maximal spreading diameter (a,b), but then the drop does not retract, it is stuck at its maximal diameter (c). This pinning is due to the formation of a thin layer of ice between the substrate and the liquid during the spreading. On images (c–e) the drop is thus made of a thin layer of ice, attached to the substrate and growing vertically, beneath a thicker layer of water where capillary waves are damped. This remaining water layer is not stable due to its high aspect ratio, so that it still needs to retract in order to reach its equilibrium, and we observe on (e) where the liquid layer has actually just started to unpin from the edge. Then, the liquid retracts on the ice layer (f), until it reaches its typical equilibrium contact angle on ice (Knight Reference Knight1971) and forms a spherical-cap-shaped drop (g). Eventually, it completely freezes, yielding a pointy ice drop (Anderson, Worster & Davis Reference Anderson, Worster and Davis1996; Snoeijer & Brunet Reference Snoeijer and Brunet2012; Marin et al. Reference Marin, Enriquez, Brunet, Colinet and Snoeijer2014), on top of an ice pancake (h).

Figure 4.

Snapshots (a–h) show a sequence at different times (indicated on each image) of a drop of water of radius 1.9 mm impacting at velocity $2.6~\text{m}~\text{s}^{-1}$ on an aluminium substrate at temperature $-9\,^{\circ }\text{C}$ . (i) is a height profile of the final frozen splat, measured using optical profilometry. The aspect ratio is 20, so the real splat is much flatter.

Figure 4(i) shows the height profile of the frozen drop (h) obtained by scanning the drop diameter with an optical profilometer. This profile shows clearly two different parts: a quasi-cylindrical plate (the so-called ice pancake), whose thickness is denominated as  $h_{p}$ , on top of which we find an ice pattern, the consequence of the water retraction on ice. We thus define clearly $h_{p}$ as the ice layer thickness of the edge of the structure when it unpins. This structure with an ice pancake in contact with the impacted solid is found any time a drop impacts a sub-zero cold surface. The ice pattern that is eventually formed on the top of the ice pancake can exhibit different shapes that will be the subject of future work. Here we focus on the ice pancake that is, in fact, crucial in all relevant applications, in particular for the stability of the frozen drop. Indeed, as the temperature of the ice varies, thermal expansion/contraction generates elastic stress that can lead to the delamination of the pancake (de Ruiter et al. Reference de Ruiter, Soto and Varanasi2018) or the formation of cracks (Ghabache et al. Reference Ghabache, Josserand and Séon2016), influenced by the thickness  $h_{p}$ . The aim of the present paper is thus dedicated to the quantification of this thickness  $h_{p}$ .

3.2 The stationary contact line (SCL) regime

As we observe on figure 4, after spreading, the contact line seems steady (b–e), before a dewetting transition occurs leading to the retraction of the water film on the ice layer. It is thus expected that the thickness of the ice pancake $h_{p}$ is built by thermal conduction during this time. Consequently, we need first to characterize the duration of this regime, that we call SCL regime (Rivetti et al. Reference Rivetti, Salez, Benzaquen, Raphaël and Bäumchen2015). Figure 5(a) presents the variation of the liquid film radius during an impact, plotted as a function of time, for three different substrate temperatures. The dashed line shows the liquid drop radius evolution at room temperature and the two solid lines at freezing temperature, $-14$ and $-30\,^{\circ }\text{C}$ . Without freezing (dashed line), the drop spreads rapidly, reaches its maximal diameter and almost instantaneously retracts. The behaviour is different when the drop freezes. Indeed, after having reached its maximal diameter (illustrated on the curve with the inserts), the stationary contact line regime is observed: the liquid remains attached to the ice layer close to its maximum radius. There, the liquid radius barely varies with time during approximately 60–70 ms before the retraction of the liquid film starts. This time between the spreading and the retraction regimes defines the SCL time, $\unicode[STIX]{x1D70F}_{SCL}$ , as shown on the curves. We note that the two SCL times, for the two different temperatures, seem to be quite close.

Figure 5.

(a) The dynamics of spreading and retraction of the drop. The radius of the liquid film is plotted as function of time. Inset pictures show the aspect of the film on the different stages. On all graphs, impact velocities are represented by colours (light blue for the slowest to dark purple for the fastest), substrates by symbols (circles for copper, squares for steel and triangles for marble) and drop radii by the filling (full symbols for $R=1.9$  mm, empty symbols for $R=1.2$  mm). (b) Here $\unicode[STIX]{x1D70F}_{SCL}$ as a function of the temperature $\unicode[STIX]{x0394}T$ for different impact velocities. The drop radius is $1.9$  mm and the substrate is steel. (c) Here $\unicode[STIX]{x1D70F}_{SCL}$ as a function of the temperature $\unicode[STIX]{x0394}T$ for different substrates. The impact velocity is $2.6~\text{m}~\text{s}^{-1}$ and the drop radius is 1.9 mm (full symbols). Results for a drop of $R=1.2$  mm impacting on steel at the same velocity are shown with empty squares. (d) Theoretical prediction for the water film thickness $h_{w}$ as a function of its experimental estimation $h_{tot}^{exp}$ . Each point represents a series of experiments with given impact velocity, substrate and drop radius. The error bars contain the variations within each series. The dashed line is a linear fit with the slope set to 1.

Figure 5(b) shows the dependence of the time $\unicode[STIX]{x1D70F}_{SCL}$ with the temperature, for different heights of fall and same drop size (1.9 mm) and substrate material (steel). Instead of the substrate temperature $T_{s}$ , we use $\unicode[STIX]{x0394}T=T_{m}-T_{s}$ , where $T_{m}$ is the freezing temperature of the liquid (here for water $T_{m}=0\,^{\circ }\text{C}$ so that $\unicode[STIX]{x0394}T=-T_{s}$ when using the Celsius scale), which is the relevant temperature here, because we saw that the energy needed to cool the drop down to $T_{m}$ can be neglected. The first observation is that $\unicode[STIX]{x1D70F}_{SCL}$ always reaches a plateau where it is independent of $\unicode[STIX]{x0394}T$ . The value of this plateau strongly depends on the drop impact velocity and thus on the initial spreading of the liquid film: the larger the maximal diameter at impact, the shorter the liquid layer stays stationary atop the ice pancake before retracting. These results suggest therefore that $\unicode[STIX]{x1D70F}_{SCL}$ might be independent of the heat transfer and due solely to the liquid dynamics.

Figure 5(c) presents the variation of $\unicode[STIX]{x1D70F}_{SCL}$ with $\unicode[STIX]{x0394}T$ for our three different substrates, marble (triangle markers), steel (square) and copper (circle), and for two different drop radii on steel substrate (full and empty squares). In each case, the plateau regime is quickly reached confirming this independence of SCL time with $\unicode[STIX]{x0394}T$ . As expected, the drop size, as the drop impact velocity, seems to play a role in the selection of the plateau value of $\unicode[STIX]{x1D70F}_{SCL}$ . What is more surprising here is the variation of the plateau value with the substrate material, $\unicode[STIX]{x1D70F}_{SCL}$ decreases as the substrate thermal conductivity increases. This appears as a contradiction with the previous result that $\unicode[STIX]{x1D70F}_{SCL}$ does not depend on the substrate temperature and, thus, seems to be independent of the thermal parameters. We will come back to these observations in the following.

To explain the existence of this delay $\unicode[STIX]{x1D70F}_{SCL}$ between spreading and retraction when water freezes, we will use an argument of a previous work (Rivetti et al. Reference Rivetti, Salez, Benzaquen, Raphaël and Bäumchen2015) on the relaxation of a contact line pinned at the edge of a polymer film. Similarly, we assume here that the time $\unicode[STIX]{x1D70F}_{SCL}$ is due to the relaxation dynamics of the contact angle $\unicode[STIX]{x1D703}$ formed by the liquid film and the growing ice pancake, and pinned at its edge. The contact line relaxation starts at an angle given by the spreading dynamics (influenced by the drop impact velocity, the substrate wetting angle, and the drop size), and ends when $\unicode[STIX]{x1D703}$ reaches a critical value $\unicode[STIX]{x1D703}^{\ast }$ enabling its unpinning, a required prelude to its receding motion. In the experiment of Rivetti et al. (Reference Rivetti, Salez, Benzaquen, Raphaël and Bäumchen2015), this angle is $\unicode[STIX]{x1D703}^{\ast }=4.5^{\circ }\pm 0.5^{\circ }$ and surprisingly appears to be independent of the liquid and of the film thickness.

Because the liquid film is very thin, its dynamics can be taken in the lubrication regime so that the contact line relaxation is expected to follow a capillaro-viscous relaxation time $t_{w}\propto (h_{0}\unicode[STIX]{x1D702})/\unicode[STIX]{x1D6FE}$ , with the proportional coefficient being a function of $\unicode[STIX]{x1D703}^{\ast }$ , being about $10^{5}$ in their experiments. Note that, as soon as the contact line retracts, they showed that $\unicode[STIX]{x1D703}$ increases to a receding contact angle which stays roughly constant during dewetting.

In order to show that the same relaxation dynamics is at play in our experiments, we compare on figure 5(d) the characteristic thickness $h_{w}$ of the liquid film that corresponds to the relaxation time $\unicode[STIX]{x1D70F}_{SCL}$ following this capillaro-viscous dynamics, i.e. $h_{w}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70F}_{SCL}/(10^{5}\unicode[STIX]{x1D702})$ , and the characteristic thickness of the water film $h_{tot}^{exp}$ formed at the end of the spreading, neglecting ice formation, $h_{tot}^{exp}=V_{tot}/\unicode[STIX]{x03C0}R_{max}^{2}$ with $V_{tot}$ the total volume of the drop and $R_{max}$ the spreading radius. Remarkably, the experimental data gather along a line of slope 1, suggesting that this relaxation scenario of the contact line is correct. Note that, since $h_{tot}^{exp}$ corresponds to the total height of the spread drop, it is the sum of the ice pancake and the water film, so that the intersection of the dotted line with the $x$ -axis gives a consistent estimation of the ice pancake thickness. We point out the fact that the thicknesses of the ice and liquid layer are varying in time, so that $h_{tot}^{exp}$ is just a rough approximation of the liquid film thickness. This might explain the data scattering around the linear prediction.

We have used for comparison the same value of $\unicode[STIX]{x1D703}^{\ast }$ as for the polymer films (leading to the $10^{5}$ prefactor for  $h_{w}$ ), while there is no reason that this value is relevant for water films. In fact, when trying to fit the best line for the data of figure 5(d), we found an angle slightly different but with no significant improvement when comparing with the $\unicode[STIX]{x1D703}^{\ast }$ correlation, so that we have kept this value for the sake of simplicity. The fact that this dewetting angle $\unicode[STIX]{x1D703}^{\ast }$ appears to be almost independent on the substrates, the liquids and the film thicknesses, while our experiments and those of Rivetti et al. (Reference Rivetti, Salez, Benzaquen, Raphaël and Bäumchen2015) concern film with totally different thickness is a very interesting and intriguing result. It might shed light on a universal dewetting mechanism that deserves specific investigation in the future.

Finally, following these results, the stationary contact line time ( $\unicode[STIX]{x1D70F}_{SCL}$ ) varies linearly with the thickness of the liquid pancake at impact ( $h_{tot}^{exp}$ ) and can be roughly expressed by the unique formula

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{SCL}=10^{5}\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D6FE}}\left(70+\frac{V_{tot}}{\unicode[STIX]{x03C0}R_{max}^{2}}\right).\end{eqnarray}$$

In this relation the quantity $R_{max}$ has been evaluated experimentally so far since no accurate modelling is available that takes into account the solidification.

As the spread drop diameter increases with the impact velocity, the liquid pancake is thinner and we can understand the decrease of $\unicode[STIX]{x1D70F}_{SCL}$ when the impact velocity increases, observed on figure 5(b). These results also resolve the contradiction that we draw concerning the dependence of $\unicode[STIX]{x1D70F}_{SCL}$ with the different substrates. The contact line relaxation dynamics suggests indeed that the substrate intervenes through its wetting properties that influences the spreading radius, and therefore the liquid film thickness pinned on the thin ice layer. In this scenario, the thermal properties of the substrate plays no direct role, thus explaining the independence of $\unicode[STIX]{x1D70F}_{SCL}$ with the temperature.

3.3 The ice pancake

At this point, we know the time during which the ice pancake is building ( $\unicode[STIX]{x1D70F}_{SCL}$ ) and we have a validated model for the ice growth dynamics ( $D_{eff}$ ); we now need the experimental measurements of the pancake thickness, $h_{p}$ . Figure 6(a–c) shows the ice thickness $h_{p}$ , deduced from the height profiles of the frozen drops (see figure 4 i), as a function of $\unicode[STIX]{x0394}T$ , respectively, for four different drop impact velocities (a), three different substrates (b) and two different drop sizes (c). These control parameters, and their corresponding markers, are the same as those used in figure 5.

Figure 6.

(a) Evolution of the underlying ice plate thickness $h_{p}$ as a function of the temperature $\unicode[STIX]{x0394}T$ for different impact velocities, on the same substrate (steel) and with the same drop radius (1.9 mm). The colder the substrate is, the thicker the ice layer is. For a given temperature the ice layer is thinner for the highest impact velocities. (b) Evolution of $h_{p}$ along $\unicode[STIX]{x0394}T$ , for different substrates, with the same impact velocity ( $2.6~\text{m}~\text{s}^{-1}$ ) and drop radius (1.9 mm). For a given temperature, the more conductive the substrate, the thicker the ice layer. (c) Evolution of $h_{p}$ with $\unicode[STIX]{x0394}T$ for different drop radius, on the same substrate (steel) and at the same velocity ( $2.6~\text{m}~\text{s}^{-1}$ ). The smaller the drop the less it freezes. (d) Rescaling of the measurements against the model: the measured diffusivity of the solidification front, $h_{p}^{2}/\unicode[STIX]{x1D70F}_{SCL}$ is plotted against the effective diffusion coefficient obtained through the previously described model ((2.6) and (2.10)). The dashed line has a slope 1.

These graphs show that the underlying ice layer becomes thicker when: (i) the substrate is colder (figure 6 a–c); which is expected since more liquid can be frozen during the solidification time that appeared to be mostly constant with $\unicode[STIX]{x0394}T$ (figure 5 b,c). (ii) The substrate heat conductivity increases (figure 6 b); the heat is indeed transferred to the substrate with higher efficiency, so that the freezing front can propagate faster during a solidification time that, again, does not vary much (figure 5 d). (iii) The drop impact velocity is slower (figure 6 a), or the drop size is smaller (figure 6 c), because the spread drop stays freezing longer before retracting (figure 5 b,c).

We can now compare our experimental results to our theoretical model. Figure 6(d) presents the variation of $h_{p}^{2}/\unicode[STIX]{x1D70F}_{SCL}$ , which represents the experimental diffusion coefficient of the solidification front, as a function of the theoretical diffusion coefficient $D_{eff}$ , for all the substrate materials, substrate temperatures, drop sizes and impact velocities investigated experimentally. The substrate temperature ranges from $-5\,^{\circ }\text{C}$ to $-80\,^{\circ }\text{C}$ , which yields a Stefan number ranging from approximately $10^{-2}$ to 1. We observe a nice collapse of all the data into a straight line of slope 1. We can therefore conclude that the ice pancake thickness is well described by the expression

(3.2) $$\begin{eqnarray}h_{p}^{2}\simeq D_{eff}\unicode[STIX]{x1D70F}_{SCL}.\end{eqnarray}$$

Together with the formula obtained for the time $\unicode[STIX]{x1D70F}_{SCL}$ (3.1), it provides the following expression for the ice pancake thickness

(3.3) $$\begin{eqnarray}h_{p}^{2}\simeq D_{eff}\unicode[STIX]{x1D70F}_{SCL}=10^{5}D_{eff}\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D6FE}}\left(70+\frac{V_{tot}}{\unicode[STIX]{x03C0}R_{max}^{2}}\right)\end{eqnarray}$$

that accounts quantitatively for almost all of our experiments. This indicates clearly that our simplified 1-D model describes correctly the dynamics and that the ice pancake formation is controlled by the thermal diffusion during the unpinning time of the remaining liquid film.

We observe that some experiments are not collapsing on the line with the others. They appear particularly for an effective diffusive coefficient greater than $3.10^{-7}~\text{m}^{2}~\text{s}^{-1}$ . This limit of our model probably comes from the hypothesis of semi-infinite water in the freezing model. Indeed, if we estimate the thickness of the water film at impact $h_{tot}^{exp}$ with the drop volume $V_{tot}$ and the spreading radius $R_{max}$ , we find $h_{tot}^{exp}\sim V_{tot}/\unicode[STIX]{x03C0}R_{max}^{2}\sim 200~\unicode[STIX]{x03BC}\text{m}$ , with a drop of volume $V_{tot}=30~\unicode[STIX]{x03BC}\text{L}$ spreading to a radius of $R_{max}=7$  mm. This thickness is smaller than the highest values of  $h_{p}$ , meaning that almost the whole water drop freezes during $\unicode[STIX]{x1D70F}_{SCL}$ and the hypothesis of semi-infinite water does not hold anymore.