1. Introduction
In his poem De Rerum Natura, Lucretius asks: ‘don't you see, besides, how drops of water falling down against the stones at last bore through the stones?’. This sentence, dating back to the first-century BC, is a testament to the long research interest for the problem of drop impacts on a substrate, which is always of topicality nowadays (Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Blanken et al. Reference Blanken, Saleem, Thoraval and Antonini2021; Cheng, Sun & Gordillo Reference Cheng, Sun and Gordillo2022). A better understanding of the droplet dynamics and maximum spreading diameter after impact is motivated by the wide range of industrial and natural applications such as, amongst others, spray deposition (Pasandideh-Fard et al. Reference Pasandideh-Fard, Pershin, Chandra and Mostaghimi2002), aerosol generation (Joung & Buie Reference Joung and Buie2015) or raindrop erosion (Zhao et al. Reference Zhao, Zhang, Tjugito and Cheng2015). This led to the elaboration of models describing single drop impacts in the capillary and viscous limits (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010), or in the transition between these two asymptotic regimes (Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016), in an isothermal context (i.e. when no thermal effects are involved).
The particular configuration of drop impact involving phase change also received a significant attention due to its relevance for three-dimensional or inkjet printing (Wang et al. Reference Wang, Tsai, Wu and Hwang2016; Lohse Reference Lohse2022), spray coating or cooling processes (Shukla & Kumar Reference Shukla and Kumar2015; Breitenbach, Roisman & Tropea Reference Breitenbach, Roisman and Tropea2018), or aircraft icing problematics (Baumert et al. Reference Baumert, Bansmer, Trontin and Villedieu2018), for instance. Several studies focused on droplets impinging on heated walls (Moita, Moreira & Roisman Reference Moita, Moreira and Roisman2010; Quéré Reference Quéré2013; Liang & Mudawar Reference Liang and Mudawar2017) and identified different regimes for the spreading dynamics. In particular, Chandra & Avedisian (Reference Chandra and Avedisian1991) explored experimentally the case of 
$n$-heptane droplets impinging at a low impact velocity on a stainless steel substrate whose surface temperature could be varied from 
$24\,^{\circ }{\rm C}$ to 
$250\,^{\circ }{\rm C}$, encompassing both the liquid boiling point and the Leidenfrost point. This allowed these authors to describe the spreading process below and above the Leidenfrost point. In the second situation, the impacting droplet levitates above the substrate, due to the formation of a vapour layer under the expanding liquid film. From there, Tran et al. (Reference Tran, Staat, Prosperetti, Sun and Lohse2012) investigated water drop impacts on hot surfaces and provided a comprehensive phase diagram highlighting the existence of three regimes: contact boiling, gentle film boiling and spraying film boiling. Staat et al. (Reference Staat, Tran, Geerdink, Riboux, Sun, Gordillo and Lohse2015) studied the impact of ethanol droplets on a hot surface, varying both the Weber number and the substrate temperature, to determine the transition towards splashing and the dynamic Leidenfrost point (onset of the Leidenfrost effect). They evidenced a strong dependency of the splashing threshold with the substrate temperature. The transition regime between contact boiling and film boiling has been investigated by Shirota et al. (Reference Shirota, van Limbeek, Sun, Prosperetti and Lohse2016) using total internal reflection imaging. Another experimental contribution from Antonini et al. (Reference Antonini, Bernagozzi, Jung, Poulikakos and Marengo2013) revealed that the Leidenfrost effect, superhydrophobicity and sublimation of the substrate have a similar influence on the dynamics of an impinging drop, with droplet rebound being observed in each situation. Liu, Cai & Tsai (Reference Liu, Cai and Tsai2020) recently studied drop impacts on heated nanostructures and highlighted that hot nanotextures can enhance jetting and splashing during the impact process.
Although the case of drop impact on a hot wall has been the subject of important scientific literature, fewer experimental studies have been dedicated to situations featuring liquid solidification or substrate melting. In a seminal contribution, Madejski (Reference Madejski1976) derived a theoretical analysis of the spreading dynamics of a liquid droplet impinging onto a solid substrate cold enough to trigger solidification, based on energy conservation, alongside experiments of metal drop impacts performed on different substrates. Numerous studies have since been dedicated to the case of a metal drop impinging on a cold substrate, providing estimates for the maximum spreading diameter (Pasandideh-Fard et al. Reference Pasandideh-Fard, Bhola, Chandra and Mostaghimi1998; Gielen et al. Reference Gielen, de Ruiter, Koldeweij, Lohse, Snoeijer and Gelderblom2020), or evidencing the intriguing self-peeling phenomenon occurring for a cold enough surface temperature (de Ruiter, Soto & Varanasi Reference de Ruiter, Soto and Varanasi2018). When studying the outcome of water droplet impacts on a cold substrate, Ghabache, Josserand & Séon (Reference Ghabache, Josserand and Séon2016) observed different crack patterns developing in the resulting frozen puddle, depending on the surface temperature, and proposed a model to estimate the thresholds towards the fragmentation and the hierarchical regimes. Schremb, Roisman & Tropea (Reference Schremb, Roisman and Tropea2018) studied the impact of supercooled water droplets on a smooth ice target and developed an analytical framework to describe the lamella thinning as well as the final ice thickness. Recently, Thiévenaz, Séon & Josserand (Reference Thiévenaz, Séon and Josserand2020) studied the influence of solidification (or freezing) on the maximum spreading of a water droplet impacting on a cold surface, at rather large impact velocities. These authors proposed a model of effective viscosity that allowed them to successfully describe their experiments. If these studies focused on the sole case of solidification during a drop impact, recent experimental works investigated the thermodynamic configuration in which the droplet is able to melt the solid surface it impacts (Jin, Zhang & Yang Reference Jin, Zhang and Yang2017; Ju et al. Reference Ju, Yang, Yi and Jin2019; Lolla et al. Reference Lolla, Ahmadi, Park, Fugaro and Boreyko2022), but a model describing the effect of substrate melting on the impact outcome is currently missing.
Therefore, a unified description of the influence of phase change on the maximum spreading diameter resulting from a drop impact remains elusive. In particular, although these situations are conceptually close, no modelling describes the effects of substrate melting and liquid solidification on the impact outcome in a common framework. These aspects motivated the present experimental study, which aims at investigating the case of a temperature-controlled water droplet falling onto a cold substrate, made either of ice or of cold brass, in order to understand the effects of substrate melting and liquid solidification on the maximum spreading diameter.
2. Experimental set-up and methods
The experimental set-up designed to this end is schematized in figure 1(a). It consists of a 
$2.2\ {\rm {m}}$ high vertical beam, that holds an aluminium block that can be set at an adjustable height. A vertical needle of outer diameter 
$1.83\ {\rm {mm}}$ passes through the block, which was hollowed out to host a heating cartridge placed in contact with the needle and connected to a generator. At the bottom of the beam, a cooling unit can be used as is or to generate an ice layer: it is made of a thin brass plate, cooled by a Peltier heat sink that is in turn connected to a cold bath operating with a mixture of water and ethylene glycol. Drop impacts can then be made directly onto cold brass, to study the role of liquid solidification during the impingement process, or onto an ice layer, which is produced by depositing a certain amount of water on the plate before freezing it rapidly with the help of the Peltier modulus. When conducting experiments with the brass substrate, a perspex plate is placed on top of it to limit the formation of frost, without altering much its surface temperature due to the low thermal effusivity of plastic. This protection is removed just before performing a drop impact. A thermocouple, placed on top of the substrate, allows us to adjust the heat flux imposed by the Peltier so that the brass or ice layer is set at an initial surface temperature 
$T_s$. At the beginning of an experiment, the needle is positioned so that its tip is located at a distance 
$H$ from the substrate (with a 
$4$ mm accuracy). Then, the liquid is gently pushed into the needle using a syringe pump and a pendant drop is formed. The droplet temperature is set and maintained at a controlled value 
$T_d$ using a regulation loop, based on a tension generator connected to the heating cartridge and to two thermocouples that are respectively placed inside the needle and in contact with the cartridge. As long as the temperature in the needle is lower than 
$T_d$, the cartridge is powered and heats the liquid up, but it is switched off when the measured temperature becomes greater than or equal to 
$T_d$. This regulation system allows us to quickly obtain a very stable initial temperature for the pendant drop, which is measured with a 
$\pm 2\,^{\circ }{\rm C}$ accuracy. The syringe pump is then activated anew to slowly inject liquid inside the drop so that it eventually reaches its critical volume and detaches from the needle under the action of gravity. As a result, a droplet of initial diameter 
$D_0$ starts its fall over the vertical distance 
$H$ before impacting and spreading over the ice or brass substrate. This process is recorded from above by a high-speed camera, which operates at 5000 fps.
Figure 1. (a) Schematic representation of the experimental set-up. The initial diameter of the liquid droplet is 
$D_0$ and the distance between its original position and the substrate is 
$H$. (b–e) Photographs of a transparent water droplet spreading on melting tinted ice, with 
$T_s=-5\,^{\circ }{\rm C}$, 
$T_d=25\,^{\circ }{\rm C}$ and 
$H=0.4$ m, at times (b) 
$t=0$, (c) 
$t=3.5$ ms, (d) 
$t=10$ ms and (e) 
$t=35$ ms after impact. ( f–i) Pictures of a tinted water droplet spreading on cold brass, with 
$T_s=-32.2\,^{\circ }{\rm C}$, 
$T_d=19.3\,^{\circ }{\rm C}$ and 
$H=1.5$ m, at times ( f) 
$t=0$, (g) 
$t=1$ ms, (h) 
$t=3.6$ ms and (i) 
$t=35$ ms after impact. In (b–e) and ( f–i), melting and solidification are evidenced by the increase or decrease of the brightness with time, respectively, which reveals that more and more ice melts (respectively, an increasing amount of water solidifies). In (d) the red line is the extracted contour of the liquid film when it reaches its maximum radial extent. Circles of diameter 
$D_{min}$ (yellow dashed line) and 
$D_{max}$ (white dash-dotted line) corresponding to the averaged minimum and maximum droplet diameter, respectively, are also represented. Horizontal bars in (b, f) give the scales for each corresponding image sequence.
In most experiments reported in the present study, the liquid used is pure water, whose initial temperature has been varied between 
$18\,^{\circ }{\rm C}$ and 
$80\,^{\circ }{\rm C}$, while the substrate temperature 
$T_s$ ranged between 
$-33\,^{\circ }{\rm C}$ and 
$-2\,^{\circ }{\rm C}$ for both the ice or brass surface. The height 
$H$ has been explored in the range 
$[ 0.04\unicode{x2013}2.2 ]$ m so that the resultant impact velocity 
$U$ of the water droplet, which is evaluated by a home-made code accounting for the air resistance, is varied between 
$0.9\ {\rm {m}}\ {\rm {s}}^{-1}$ and 
$5.9\ {\rm {m}}\ {\rm {s}}^{-1}$.
A complementary set of experiments of dimethyl sulfoxide (DMSO) drop impacts on cold brass have also been performed, as this fluid has significantly different thermal properties than water and also a reduced surface tension. For these particular tests, the droplet temperature was kept constant (at 
$T_d=25\,^{\circ }{\rm C}$), whereas three substrate temperatures (
$T_s=-30\,^{\circ }{\rm C}$, 
$T_s=10\,^{\circ }{\rm C}$ and 
$T_s=30\,^{\circ }{\rm C}$) and four initial heights (
${H=10}$ cm, 
$H=50$ cm, 
$H=128$ cm and 
$H=170$ cm) have been investigated. The explored surface temperatures have been chosen so as to have one temperature above the DMSO freezing point (which is 
$T_f \simeq 18.6\,^{\circ }{\rm C}$) and two below. As a result, these experiments solely involve no phase change or liquid solidification.
For each initial condition, experiments are repeated three times to ensure reproducibility. Using standard correlations for the convective heat transfer of a sphere moving in a fluid (Yuge Reference Yuge1960), the maximum temperature drop during the fall of a water droplet can be estimated to be less than 
$2\,^{\circ }{\rm C}$ at worst, i.e. when 
$T_d=80\,^{\circ }{\rm C}$ and 
$H=2.2$ m, as described in Appendix B. In the case of experiments with DMSO, for which 
$T_d$ is very close to the ambient air temperature, the estimated change is below 
$1\,^{\circ }{\rm C}$. As a result, we choose to neglect this effect so that the droplet temperature at the time of impact is considered to be equal to 
$T_d$. Finally, the initial droplet diameter 
$D_0$ has been measured for all experiments from the last image showing the drop before collision with the substrate: it is found to be 
$D_0 \simeq 4.0 \pm 0.15\ {\rm {mm}}$ for water and 
$D_0 \simeq 3.2 \pm 0.15\ {\rm {mm}}$ for DMSO.
3. Qualitative and quantitative results
In the present experiments, when a liquid droplet impacts the substrate, it starts spreading radially on top of it. This dynamics happens on the characteristic kinetic time scale 
$D_0/U$, of the order of milliseconds in the present experiments, until the drop reaches its maximum spreading diameter. After this moment, there is no significant retraction of the contact line on the substrate.
An important aspect is to determine whether the ice effectively melts on the same time scale when a hot water droplet impacts on its surface, or if the liquid (water or DMSO) solidifies when impinging a cold brass substrate. To this end, two kinds of preliminary experiments are performed with water droplets, where either the solid or liquid phase is dyed with fluorescein. The spreading is illuminated from above using a UV-blue light (with a wavelength of 
$470\ {\rm {nm}}$) and the camera lens is covered by a green optical filter (with cutting wavelength of 
$495\ {\rm {nm}}$), so that only fluorescent regions appear bright on the obtained images. In the first situation, an ice layer is dyed during its formation, whereas the impinging drop remains translucent. As fluorescence does not happen when fluorescein molecules are diluted in solid water due to a self-quenching phenomenon (Huerre et al. Reference Huerre, Monier, Séon and Josserand2021), its detection is a signature of melted water originating from the ice layer. The pictures in figures 1(b)–1(e), corresponding to an experiment for which 
$T_s=-5\,^{\circ }{\rm C}$, 
$T_d=25\,^{\circ }{\rm C}$ and 
$H=0.4$ m, show (b) the drop prior to the collision with the substrate, (c) the liquid film during the spreading process, (d) the moment the maximum diameter is reached and (e) the final footprint left by the impact. If the initial droplet is almost not visible due to its transparency, the liquid phase becomes increasingly luminous as time goes by: this demonstrates the melting of the ice during the whole impact process. Conversely, in the second case corresponding to a water drop impact onto cold brass, only the impinging liquid has been dyed. Figure 1( f–i) presents a typical image sequence of this configuration for an experiment where 
$T_s=-32.2\,^{\circ }{\rm C}$, 
$T_d=19.3\,^{\circ }{\rm C}$ and 
$H=1.5$ m, with ( f) the tinted droplet just before the first contact with the substrate and (i) the resulting imprint long after the spreading phase. The fact that the intensity in the liquid decreases with time, with some part of the expanding droplet becoming increasingly dark, reveals that solidification of the liquid layer is at play. It should be underlined that these phenomena are also visible during spreading: this suggests that the two typical times of radial expansion and of phase change are of the same order, so that there is no scale separation between the two processes. Another comment arises from the images of figure 1: the centre of the spreading water droplet pictured in figure 1(h), photographed at the moment the maximum diameter has been reached for the corresponding experiment, remains fluorescent (hence, some liquid remains) and not completely dark (hence, solid) as in figure 1(i). This shows that the arrest criterion does not correspond to the moment the liquid–solid interface is reached by the droplet free surface.
In order to study quantitatively the impact outcome, we performed drop impact experiments on ice or cold brass, in which only the liquid (water or DMSO) has been tinted with fluorescein and the parameters 
$T_d$, 
$T_s$ and 
$H$ have been varied systematically in the ranges indicated in § 2. In what follows, we specifically focus on the maximum diameter reached by the spreading liquid film that, under the presence of solidification or substrate melting, has received little attention so far in the literature. It should be mentioned, though, that the transient dynamics has already been addressed by several previous studies (Jin et al. Reference Jin, Zhang and Yang2017; Ju et al. Reference Ju, Yang, Yi and Jin2019; Thiévenaz et al. Reference Thiévenaz, Séon and Josserand2020), showing that the spreading dynamics itself was only smoothly affected by the phase change.
The final contour of the contact line, obtained when the liquid film reaches its maximum radial extent, can then be extracted as illustrated by the red solid line in figure 1(d). Image processing allows us to locate the positions of the local maxima (i.e. the tip of the digitations) and minima (located between two fingers), relative to the centre of mass of the contour (
$O$). From the corresponding radial distances to 
$O$, the minimum and maximum droplet diameters, noted respectively as 
$D_{min}$ and 
$D_{max}$, are defined as the averaged positions of the local minima and maxima, respectively. From all our experiments, we observe a linear relationship between 
$D_{max}$ with 
$D_{min}$ regardless of the droplet or substrate temperatures, which reads 
$D_{max}=1.07D_{min}$ for water and 
$D_{max}=1.04D_{min}$ for DMSO. As a result, the more the spreading, the more the fingers’ elongation. The proportionality between the two diameters is an intriguing result, which suggests, for instance, that for water, the typical size of the digitations 
$D_{max}-D_{min}$ is about 
$7\%$ of the spreading diameter 
$D_{min}$ and that 
$D_{max}$ can be described in a similar way as 
$D_{min}$. As a result, the maximum spreading ratio is defined as 
$\beta _m \equiv D_{min}/D_0$, with 
$D_0$ the initial diameter of the droplet.
In figure 2, 
$\beta _m$ is presented as a function of the impact velocity 
$U$ for experiments involving water with a fixed (a–c) droplet temperature 
$T_d$ or (d–f) substrate temperature 
$T_s$, which are indicated above each plot. The markers correspond to experiments involving water, and for which (
$\lozenge$) solidification or (
${\bigcirc}$) fusion (substrate melting) occurs. As it was not always straightforward to distinguish which kind of phase change was at play for a given experiment, the marker to assign is determined by the sign of the moving ice–liquid front position, which is calculated from the results presented in § 4. In all cases, the spreading ratio clearly increases with 
$U$. Furthermore, for a fixed value of 
$T_d$, increasing 
$T_s$ results in larger 
$\beta _m$: the higher the substrate temperature, the larger the spreading ratio. This is illustrated for ice in (b,c) where 
$T_d=50\,^{\circ }{\rm C}$ and 
$T_d=80\,^{\circ }{\rm C}$, respectively: in these situations, experiments conducted at 
$T_s=-2\,^{\circ }{\rm C}$ (light grey) are significantly above those performed at 
$T_s=-25\,^{\circ }{\rm C}$ (blue). To a lesser extent, at a given value of 
$T_s$, 
$\beta _m$ is larger when 
$T_d$ is increased. This is especially visible in ( f) for 
$T_s=-2\,^{\circ }{\rm C}$, where a gentle order exists with the value of the droplet temperature 
$T_d$. These observations agree with the results gathered by Jin et al. (Reference Jin, Zhang and Yang2017) and Ju et al. (Reference Ju, Yang, Yi and Jin2019) for drop impacts on ice, and by Thiévenaz et al. (Reference Thiévenaz, Séon and Josserand2020) for droplet impingement on a cold metal substrate. Experiments with DMSO present a similar evolution as the water drop impact tests presented in figure 2(a). They are not included here as they feature a slightly larger initial droplet temperature (
$T_d=25^{\circ }{\rm C}$).
Figure 2. Evolution of the maximum spreading ratio 
$\beta _m \equiv D_{min}/D_0$ as a function of the impact velocity 
$U$, for water drop impacts on (a) brass and (b–f) on ice. In (a–c) the droplet temperature is fixed at (a) 
$T_d=19\,^{\circ }{\rm C}$, (b) 
$T_d=50\,^{\circ }{\rm C}$ and (c) 
$T_d=80\,^{\circ }{\rm C}$, while the colourbar denotes the substrate temperature 
$T_s$. Contrariwise, in (d–f) the substrate temperature is fixed at (d) 
$T_s=-25\,^{\circ }{\rm C}$, (e) 
$T_s=-15\,^{\circ }{\rm C}$ and ( f) 
$T_s=-2\,^{\circ }{\rm C}$, while the markers’ colours represent this time the droplet temperature 
$T_d$. The symbols correspond to experiments with water where (
$\lozenge$) solidification or (
${\bigcirc}$) fusion occurs. The nature of the substrate is indicated above each plot.
4. Discussion and modelling
Building upon previous studies dedicated to drop impacts (Madejski Reference Madejski1976; Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Lagubeau et al. Reference Lagubeau, Fontelos, Josserand, Maurel, Pagneux and Petitjeans2012; Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016), two dimensionless numbers can be defined to describe the outcome of the impact process. On one hand, the Weber number 
${We} \equiv \rho U^2 D_0 / \gamma$ compares inertia to capillarity, with 
$\rho$ and 
$\gamma$ the density and surface tension of the drop, respectively. The values for 
$\rho$ and 
$\gamma$ involved in the expression of 
${We}$ are taken here at 
$T_d$ using standard correlations (see Appendix A). This choice was made as both the initial kinetic energy of the droplet and the surface energy once the maximum diameter has been reached are expected to involve liquid volume (respectively, surface) set at that temperature. Indeed, one can show that the thermal boundary layer in the liquid is much smaller than the liquid film thickness when the maximal diameter is reached. On the other hand, the Reynolds number, defined here as 
${Re} \equiv U D_0/\nu _f$, with 
$\nu _f$ the kinematic viscosity, compares inertia to viscous effects. This time, 
$\nu _f$ is evaluated at the melting point 
$T_f$ since, in the spreading dynamics, this is the viscosity close to the substrate (thus, near 
$T_f$) that is relevant. Although this might seem surprising, at first glance, as the fluid properties for 
${We}$ have been evaluated at 
$T_d$, taking the kinematic viscosity at the initial temperature of the droplet, 
$T_d$, resulted in a significantly enhanced scattering of our results, which clearly suggests that the melting point is more relevant to describe the typical temperature of the dissipative layer for our experiments.
Except for a few tests that have a small falling distance 
$H$, most of the impact velocities in the present experiments are greater than or equal to 
$2\ {\rm {m}}\ {\rm {s}}^{{-1}}$. For these data 
$230 < {We} < 2000$ and 
$4500 < {Re} < 13\,100$: as 
${Re} \gg 100$, the impact outcome is thus expected to be closer to the inertial-viscous regime provided by Madejski (Reference Madejski1976) than to the inertial-capillary regime, and hence, our data to be relatively well parameterised by the Reynolds number. The evolution of 
$\beta _m$ with 
${Re}$ is illustrated in figure 3(a). Overall, the spreading ratio is found to increase with the Reynolds number, but it can be noted that the data are scattered in this representation. Indeed, an order appears with 
$T_s$, which is visible, for instance, for water drop impacts in which ice melts (
${\bigcirc}$), as illustrated by the inset of figure 3(a). In addition, for a given Reynolds number, experiments featuring water solidification (
$\lozenge$) have a systematically smaller value for 
$\beta _m$ than data involving ice melting (
${\bigcirc}$). Last but not least, DMSO drop impacts (
), featuring either no phase change or liquid solidification, lie systematically above experiments with water solidification (
$\lozenge$).
Figure 3. (a) Maximum spreading ratio, 
$\beta _m$, as a function of the Reynolds number, 
${Re}$. The colourbar indicates the substrate temperature, while the symbols denote experiments involving (
$\lozenge$) water droplets experiencing freezing, (
${\bigcirc}$) water droplets causing ice melting and (
) DMSO droplets. (b–c) Schematic views of the ice (b) growth or (c) melting during film spreading. Here 
$\delta _\nu$ corresponds to the size of the viscous boundary layer, 
$h$ is the position of the substrate and 
$\delta _{eff}$ the size of the effective boundary layer; 
$u_r$ is the radial velocity field.
The poor collapse of the data in this representation is, in fact, expected as the spreading dynamics is affected by the presence of phase change. In the case of isothermal drop impacts belonging to the inertial-viscous regime, it has been shown that the arrest criterion corresponds to the moment the viscous boundary layer reaches the free surface of the expanding liquid film (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010). Building upon this, Thiévenaz et al. (Reference Thiévenaz, Séon and Josserand2020) evidenced the fact that, for experiments involving solidification, spreading appears to stop when the sum of the growing ice thickness and the developing viscous boundary layer reaches the free-surface elevation of the expanding droplet. Further possible evidence of such an arrest criterion may be found in the study performed by Pasandideh-Fard et al. (Reference Pasandideh-Fard, Bhola, Chandra and Mostaghimi1998), who studied tin drop impacts onto a cold stainless steel substrate. Indeed, in the numerical simulations conducted by these authors, it can be observed that spreading stops when the solidified layer at the centre of the splat approaches the free surface, whereas other regions of the spreading film remain in the liquid state.
These results from previous studies shed light on the relevance to predict the moment the viscous boundary layer reaches the free surface of the liquid film when solidification or substrate melting occurs. In such a situation, the ice grows or melts at the base of the expanding liquid film, as illustrated in figures 3(b)–3(c), thus changing the position of the solid surface on which the viscous boundary layer of thickness 
$\delta _{\nu }$ develops, thereby modifying the value of 
$\beta _m$. The position of the moving interface can, at first order, be modelled by solving the classical Stefan problem, if one assumes that the influence of advection is negligible. The importance of this latter effect is known to depend on the Prandtl number 
${Pr}$, which compares the typical sizes of the thermal and viscous boundary layers (Roisman Reference Roisman2010). In the case of water close to its melting point, this number is quite large (
${Pr} \sim 14$ at 
$0.01\,^{\circ }{\rm C}$) and so advection can safely be neglected in the absence of a phase change (Moita et al. Reference Moita, Moreira and Roisman2010). We assume that the presence of solidification or melting does not change this result much. Then, the position of the moving interface can be estimated by solving the heat diffusion equations in all phases (liquid, ice and possible substrate) alongside the so-called ‘Stefan condition’, which states that the velocity of the phase change front is directly related to the thermal flux difference at the boundary. The calculation leads to the derivation of a self-similar solution for the temperature field and the liquid–ice front position (see Appendix C for more details). To put it in a nutshell, under these assumptions the front of the ice layer 
$h$ follows a diffusive law of the form 
$h(t)= s \sqrt {\alpha _{eff}t}$, where 
$s=-1$ in the case of melting (respectively, 
$s=1$ for solidification), and 
$\alpha _{eff}$ is an effective thermal diffusivity (
$\alpha _{eff} \geqslant 0$). By introducing 
$\chi =s\sqrt {\alpha _{eff}/\alpha _i}$, with 
$\alpha _i$ the ice thermal diffusivity, 
$\alpha _{\rm {eff}}$ and 
$s$ are found numerically by solving the transcendental equation on 
$\chi$,
\begin{equation} \frac{\chi \sqrt{\rm \pi}}{2 {St}} = \frac{{\rm e}^{-\chi^2/4}}{r_i/r_s+\operatorname{erf} \left( \chi/2 \right)} + \frac{r_d}{r_i} \frac{{\rm e}^{-\chi^2/(4 \omega_d)}}{1-\operatorname{erf} \left( \chi/(2\sqrt{\omega_d}) \right)} \frac{T_f-T_d}{T_f-T_s}, \end{equation}
with 
$r_d$, 
$r_i$ and 
$r_s$ the thermal effusivities of the liquid, ice and possible substrate, respectively; 
$\omega _d=\alpha _d/\alpha _i$, with 
$\alpha _d$ the liquid thermal diffusivity; and 
${St}=c_{p,i} (T_f-T_s) / \mathcal {L}_f$ the Stefan number, with 
$c_{p,i}$ the ice thermal capacity and 
$\mathcal {L}_f$ the latent heat of fusion (for the definitions of these quantities, see also Appendix A). For drop impacts on ice, it should be noted that 
$r_s=r_i$. For each experiment, the value for 
$s$ obtained when solving (4.1) indicates whether freezing (
$s=1$) or melting (
$s=-1$) occurred, so that the symbols used in figures 2, 3(a) and 5 are chosen accordingly.
From there, following the approaches developed by Eggers et al. (Reference Eggers, Fontelos, Josserand and Zaleski2010), Roisman (Reference Roisman2010) and later by Thiévenaz et al. (Reference Thiévenaz, Séon and Josserand2020), it is possible to estimate the size of the viscous boundary layer relative to the initial position of the substrate, which is expected to eventually dictate the arrest. This is done by considering the 
$r$ component of the axisymmetric Navier–Stokes equations, in the Prandtl boundary layer framework for an incompressible flow:
\begin{equation} \partial_t u_r + u_r \partial_r u_r + u_z \partial_z u_r = \nu_f \partial_z^2 u_r. \end{equation}
Here 
$u_r$ and 
$u_z$ are the radial and vertical components of the velocity field, respectively, and 
$\partial _a$ stands for partial differentiation with respect to variable 
$a$. In the inviscid case and in the absence of phase change, using the streamfunction 
$\psi$ defined from 
$u_r \equiv - \partial _z \psi /r$ and 
$u_z \equiv \partial _r \psi /r$, the solution describing the impact can be taken as 
$\psi =-r^2z/t$, corresponding to a time decreasing arrest point flow with 
$u_r=r/t$ and 
$u_z=-2z/t$. Then, in the situation of a viscous flow subjected to solidification or melting, since both the viscous boundary layer, growing from 
$z=h(t)$, and the solid–liquid front position 
$h(t)$ follow a diffusive-in-time evolution, we can consider the following ansatz for the streamfunction:
\begin{equation} \psi \equiv \sqrt{\nu_f} \frac{r^2}{\sqrt{t}} f ( \zeta ). \end{equation}
Here 
$\zeta =[z-h(t)]/\sqrt {\nu _f t}$ is the self-similar variable and 
$f$ an unknown function of 
$\zeta$. As 
$u_r=-(r/t) f^\prime (\zeta )$, 
$f^\prime$ provides an insightful description of the shape of the boundary layer. Following Thiévenaz et al. (Reference Thiévenaz, Séon and Josserand2020), we inject the expression of 
$\psi$ into (4.2), which leads to
\begin{equation} f''' ={-}f' - \frac{1}{2} \left( \zeta + s \sqrt{\frac{\alpha_{\rm{eff}}}{\nu_f}} \right) f'' - f'^2 + 2ff''. \end{equation}
The boundary conditions are a zero velocity at the solid–liquid interface 
$\zeta =0$, and the recovery of the inviscid profile at infinity: this translates into 
$f(0)=0$, 
$f'(0)=0$ and 
$f'(+\infty )=-1$, respectively. The resolution of (4.4) for these boundary conditions is achieved numerically, by using a shooting method algorithm. The evolution of 
$\zeta$ as a function of 
$-f^\prime$ (i.e. in a ‘velocity profile’-like representation) is illustrated in figure 4(a) for several representative values of 
$\sigma \equiv s \sqrt {\alpha _{eff}/\nu _f}$. For each case, 
$-f^\prime$ increases from 
$0$ at the contact with the substrate to 
$1$ for 
$\zeta \simeq 2$, where the inviscid flow solution is thereby recovered. Furthermore, the curves for different values of 
$\sigma$ depart from each other, with those for large values of this parameter reaching the asymptotic behaviour earlier, meaning that the viscous boundary layer in this case is reduced in size when compared with lower 
$\sigma$. Therefore, this shows that the typical form factor of the viscous boundary layer, which can roughly be estimated as 
$\xi \simeq -1/f^{\prime \prime }(0)$, is a function of 
$\sigma$. In other terms, the present modelling predicts a coupling between the flow and the phase change dynamics. This fact can be verified in figure 4(b), where 
$\xi$ decreases with 
$\sigma$ in a weakly nonlinear manner. In the present experiments, 
$\xi$ ranges from 0.86 (for 
$\sigma \simeq 0.35$) to 1.04 (for 
$\sigma \simeq -0.2$). These values are highlighted in figure 4(b) by the dash-dotted and dashed lines, respectively.
Figure 4. (a) Evolution of the self-similar variable 
$\zeta =[ z - h(t) ]/\sqrt {\nu _f t}$ as a function of 
$-f^\prime$ (‘velocity profile’ representation). Each curve corresponds to a given value of 
$\sigma \equiv s \sqrt {\alpha _{eff}/\nu _f}$. (b) Evolution of the form factor 
$\xi =-1/f''(0)$ of the viscous boundary layer as a function of 
$\sigma$. The dashed and dash-dotted lines highlight the upper and lower limits for 
$\xi$ covered in the present experiments, respectively.
It should be mentioned that a zero velocity condition has been imposed at the phase change front in the above analysis, although a volume-change flow actually exists at the liquid–ice interface during solidification or melting. However, a rough estimate of the induced velocity 
$v_{pc}$ gives 
$v_{pc} \simeq (\Delta \rho /\rho _d) ({\mathrm {d}h}/{\mathrm {d}t})$, with 
$\Delta \rho = \rho _d - \rho _i$ (
$\rho _d$ and 
$\rho _i$ being the liquid and ice densities taken at the melting point, respectively). Given that 
$\Delta \rho < \rho _d$ and that 
$\alpha _{eff} \ll \nu _f$, 
$v_{pc}$ can reasonably be neglected in comparison to other velocities such as, for instance, 
$u_z$ evaluated at 
$z=\sqrt {\nu _f t}$, which explains the choice to take 
$f(0)=f^\prime (0)=0$.
From this analysis, it then becomes possible to evaluate the vertical height 
$\delta _{\rm {eff}}$ reached by the viscous boundary layer compared with the initial substrate position from 
$\xi = [\delta _{eff}-h(t)]/\sqrt {\nu _f t}$. This yields
\begin{equation} \delta_{eff} = \xi \sqrt{\nu_f t} + s \sqrt{\alpha_{eff}t}. \end{equation}
Noticeably, this height displays an overall diffusive-like behaviour. Therefore, we introduce an effective water kinematic viscosity, 
$\nu _{eff}$, which is defined as 
$\delta _{eff} \equiv \xi \sqrt {\nu _{eff}t}$, so that
\begin{equation} \nu_{eff}= \left( \sqrt{\nu_f } + \frac{s}{\xi} \sqrt{\alpha_{eff}} \right)^2. \end{equation}
In the case of solidification, one obtains 
$\nu _{eff} > \nu _f$, which means that the viscous boundary layer will reach the liquid free surface sooner than for an isothermal drop impact. Freezing thus enhances dissipation and reduces the spreading diameter. Conversely, when substrate melting occurs 
$\nu _{eff} < \nu _f$: the boundary layer will meet the free surface later than for the isothermal case, so that dissipation appears to be reduced while spreading is favoured. For some experiments, the effective kinematic viscosity that is predicted differs significantly from 
$\nu _f$: for instance, for drop impacts on ice where 
$T_s=-2\,^{\circ }{\rm C}$ and 
$T_d=80\,^{\circ }{\rm C}$, one obtains 
$\nu _{eff} \simeq 0.65 \nu _f$. We stress that, in this model, the two cases of liquid solidification and substrate melting during the impact of a droplet onto its solid phase are encompassed into the same framework, which generalizes the approach followed by Thiévenaz et al. (Reference Thiévenaz, Séon and Josserand2020).
As a result of this analysis, we define an effective Reynolds number as 
${Re_{eff}} \equiv U D_0/\nu _{eff}$, i.e. based on the effective viscosity 
$\nu _{eff}$ that is evaluated by means of (4.6). Thus, 
${Re_{eff}}$ takes into account the influence of phase change on the development of the viscous boundary layer. The spreading ratio 
$\beta _m$ is shown as a function of 
${Re_{eff}}$ in figure 5(a). This representation reveals a collapse of our experimental data for water onto a master curve, regardless of the nature of the initial substrate (brass or ice) and the dynamics of the ice–water interface (melting or solidification). This shows that the effective Reynolds number better captures the physics at play than the Reynolds number, as highlighted by a comparison with figure 3(a).
Figure 5. (a) Evolution of 
$\beta _m$ with the effective Reynolds number 
${Re_{eff}} \equiv U D_0/\nu _{eff}$, with 
$\nu _{eff}$ the effective kinematic viscosity defined in (4.6). The solid line indicates 
$\beta _m=0.82\, {Re_{eff}} ^{1/5}$. (b) Plot of 
$\beta _m {Re_{eff}}^{-1/5}$ as a function of the impact parameter 
${P_{eff}} \equiv {We} \, {Re_{eff}}^{-2/5}$. The data (
$\Delta$) as well as the universal law (solid black line) obtained by Laan et al. (Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014) for isothermal drop impacts are also reported. The colourbar indicates the substrate temperature, whereas the symbols correspond to experiments involving (
$\lozenge$) water droplets freezing, (
${\bigcirc}$) water droplets causing ice melting and (
) DMSO droplets.
Nevertheless, a significant number of the present experiments are not compatible with an inertial-viscous scaling of the form 
$\beta _m \propto {Re_{eff}}^{1/5}$, as highlighted by the comparison between the measured values of 
$\beta _m$ and the solid line reported in figure 5(a). In addition, data at low 
${Re_{eff}}$ appear to be slightly more scattered, and experiments corresponding to DMSO drop impacts (
) remain above those featuring water droplets freezing (
$\lozenge$). Such a behaviour is reminiscent of the transition from the inertial-viscous regime to the inertial-capillary regime, which has been thoroughly discussed in previous studies (Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016). In the capillary limit, the initial kinetic energy of the drop, which scales as 
$\rho U^2 D^3_0$, is completely converted into surface energy that scales as 
$\gamma D^2_{min}$ (with 
$\rho$ and 
$\gamma$ the density and surface tension of the drop evaluated at 
$T_d$, respectively). As a result, one obtains the scaling 
$\beta _m \sim {We}^{1/2}$, with 
${We} = \rho U^2 D_0 / \gamma$ the Weber number (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010). Contrariwise, in the viscous regime, the kinetic energy is balanced by viscous dissipation, which leads this time to the scaling 
$\beta _m \sim {Re}^{1/5}$ (Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010). To bridge between these two asymptotic scenarios, a universal rescaling has been proposed by Laan et al. (Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014) in the context of isothermal drop impacts, in which 
$\beta _m \textit{Re}^{-1/5}$ is a function of a sole impact parameter, 
${P} \equiv {We} \, {Re}^{-2/5}$. Adopting this approach, and using 
${Re_{eff}}$ instead of 
${Re}$, we plot in figure 5(b) the evolution of the rescaled spreading ratio 
$\beta _m {Re_{eff}}^{-1/5}$ as a function of the impact parameter 
${P_{eff}} \equiv {We} \, {Re_{eff}}^{-2/5}$ for all our experiments. In addition, the data from Laan et al. (Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014) (
$\Delta$) corresponding to isothermal drop impacts are reported in figure 5(b) with 
${Re_{eff}}={Re}$ and 
${P_{eff}}={P}$. Very noticeably, drop impacts involving water solidification (
$\lozenge$), DMSO (
), as well as most experiments with water featuring substrate melting (
${\bigcirc}$) superimpose with the data of Laan et al. (Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014), and are captured by the universal empirical law
\begin{equation} \beta_m {Re_{eff}}^{{-}1/5} = \frac{\sqrt{{P_{eff}}}}{1.24 + \sqrt{{P_{eff}}}}, \end{equation}
evidenced by these authors for the isothermal case (solid black line). The typical deviation of these experiments from (4.7) is less than 10%, similar to the dispersion of the original data from Laan et al. (Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014). However, a closer inspection reveals that water drop impacts on an ice substrate at 
$T_s=-2\,^{\circ }{\rm C}$ (grey circles) slightly deviate from relation (4.7). As these experiments belong to the transition region (
${P_{eff}} \sim 10$), and as their initial substrate temperature is close to the melting point, this suggests that wettability effects could start to play a role here (Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016). Nevertheless, as the wetting of water on ice is still a subject of active research, it is not straightforward to conclude on that aspect within the present analysis.
5. Conclusion
In the present investigation, experiments of water drop impacts onto ice and cold brass were performed, in which both liquid and substrate temperatures were varied, alongside with the falling height, in order to reach a deeper understanding of the influence of melting and solidification on the impact outcome. The maximum spreading ratio is found to increase with both temperatures as well as with the impact velocity, and the typical size of the corrugations, when present, is proportional to the final radial extent of the main liquid film. Phase change results in a modification of the viscous boundary layer, thereby affecting the overall viscous dissipation occurring within the spreading droplet. Modelling this effect through the use of an effective viscosity allows us to capture the physics at play, and to relate it to a universal law developed for the isothermal configuration. Additional experiments of DMSO drop impacts onto a cold brass substrate also show promising agreement with the proposed modelling, which suggests that the approach can be generalized to other fluids. These results pave the way for a comparison with experiments of molten metal drop impacts on a cold substrate or in the presence of evaporation, which could further validate or enrich the approach followed here. A detailed study of the effects of frost on drop impacts, using a controlled humidity set-up, would also be needed to reach a better understanding of the environmental situation. Another configuration of interest, for practical applications as well as to extend the results from the present work, would be to investigate the maximum spreading diameter following non-isothermal drop impacts in the absence of phase change. Indeed, in this scenario, varying the temperature of the initial droplet or the substrate will change the value of the contact temperature, and is thus expected to affect the behaviour of both the thermal and viscous boundary layers. This should, in turn, modify the maximum spreading diameter of the liquid film.
Acknowledgements
The authors warmly thank C. Frot and A. Garcia for their help in the elaboration of the experimental set-up.
Funding
This work was partially supported by Agence de l'Innovation de Défense (AID) – via Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité (CIEDS) – (project 2021 – ICING).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author, upon reasonable request.
Author contributions
W.S. and R.G. contributed equally to the present study.
Appendix A. Thermophysical properties of the solid and liquid phases
Table 1 reports the thermal properties of the different solids and liquids involved in the present study. The values of densities, specific heats and thermal conductivities indicated are those used when solving the Stefan problem (see Appendix C) to obtain (4.1) in the main text. From these quantities, one can readily obtain the thermal diffusivity 
$\alpha _n=k_n/(\rho _n c_{p,n})$ and effusivity 
$r_n=\sqrt {k_n \rho _n c_{p,n}}$, where subscript 
$n$ corresponds either to 
$s$, 
$i$ or 
$d$ depending on whether the substrate, ice or liquid is to be considered, respectively.
Table 1. Thermal properties of the solid and liquid phases involved in the present study. Here 
$T_f$ corresponds to the melting point, 
$\mathcal {L}_f$ to the latent heat of fusion, and 
$\rho$, 
$c_p$, and 
$k$ are the density, specific heat and thermal diffusivity, respectively (with their associated subscripts 
$s$, 
$i$ or 
$d$ denoting the substrate, ice or liquid, respectively.
As stated in the main text, some quantities are also evaluated at the melting point or at the droplet temperature, 
$T_d$. To this end, and for water, standard correlations have been used to estimate densities (Pátek et al. Reference Pátek, Hrubý, Klomfar, Součková and Harvey2009), viscosities (Pátek et al. Reference Pátek, Hrubý, Klomfar, Součková and Harvey2009; Dehaoui, Issenmann & Caupin Reference Dehaoui, Issenmann and Caupin2015) and surface tensions (Pátek, Součková & Klomfar Reference Pátek, Součková and Klomfar2016). For DMSO, for which 
$T_d=25\,^{\circ }{\rm C}$, droplet density was taken as 
$\rho =1095.5\ {\rm {kg}}\ {\rm {m}}^{{-3}}$, kinematic viscosity as 
$\nu _f=1.83\times 10^{-6}\ {\rm {m}}^2\ {\rm {s}}^{{-1}}$ and surface tension as 
$\gamma =43.53\ {\rm {mN}} {\rm {m}}^{-1}$.
Appendix B. Temperature loss during the droplet's fall
For the purpose of this investigation, it is important to estimate the temperature loss of the droplet over the course of its fall. To do so, assuming the droplet temperature 
$T_w$ to be uniform spatially, we observe that its variation with time reads
\begin{equation} m c_{p,d} \frac{{\rm d}T_w}{{\rm d}t}={\rm \pi} C(t) D^2_0 \left( T_w(t) - T_a \right) \leqslant {\rm \pi}C(t)D^2_0 \left( T_{d} - T_a \right), \end{equation}
where 
$D_0$ is the initial diameter of the drop, whose shape is assumed not to change with time, 
$m$ is its mass, 
$T_d$ and 
$T_a$ are the initial temperatures of the liquid droplet and the air, respectively, 
$c_{p,d}$ is the specific heat of the liquid considered and 
$C$ is the convective heat exchange coefficient. The right-hand side of (B1) constitutes an upper bound of the actual temperature loss. Then, one may note that 
$C(t)$ can be expressed as 
${Nu}\, k_a/D_0$, with 
$k_a$ the thermal conductivity of the air and 
${Nu}$ the Nusselt number. Integrating (B1) over time, and making the conservative assumption of a pure free fall thereby leads to the following inequality for the global temperature loss 
$\Delta T$:
\begin{equation} \Delta T \leqslant \frac{6 \left( T_d - T_a \right) k_a}{\rho c_{p,d} D^2_0} \int_0^{t_f} {Nu}(t) \, {\rm d}{t}. \end{equation}
Here 
$t_f=\sqrt {2H/g}$ corresponds to the free-fall time over a vertical distance 
$H$ (with 
$g$ the gravitational acceleration) and 
$\rho$ to the liquid density at 
$T_w=T_d$. From there, it is possible to use classical correlations to express the Nusselt number, such as those established by Yuge (Reference Yuge1960) for a sphere in forced convection, yielding
$$\begin{gather} {Nu} = 2+0.493\,{Re}^{0.5}_a\ \left( 10 < {Re}_a < 1.8 \times 10^3 \right), \end{gather}$$
$$\begin{gather}{Nu} = 2+0.3\,{Re}^{0.5664}_a\ \left( 1.8 \times 10^3 < {Re}_a < 1.5 \times 10^5 \right), \end{gather}$$
where 
${Re}_a=U(t)D_0/\nu _a$ is the Reynolds number associated to the air motion, with 
$U(t)$ and 
$\nu _a$ the droplet velocity at time 
$t$ and the air kinematic viscosity, respectively.
Finally, by defining 
${Re}_c=1.8 \times 10^3$ and 
$t_c=\nu _a {Re}_c/(gD_0)$, one gets the following upper bound for the temperature loss 
$\Delta T$ during the droplet's fall:
\begin{align} \Delta T &\leqslant \frac{6 \left( T_d - T_a \right) k_a}{\rho_w c_{p,d} D^2_0} \left\{ 2\, t_f + \frac{0.986}{3} \left( \frac{gD_0}{\nu_a} \right)^{0.5} t^{3/2}_c \right.\nonumber\\ &\quad +\left. \frac{0.3}{1.5664} \left( \frac{gD_0}{\nu_a} \right)^{0.5664} \left(t^{1.5664}_f - t^{1.5664}_c \right) \right\}. \end{align}
If, then, the left-hand side of (B5) is evaluated for water and for the various initial parameters explored in the present investigation, it is found that the maximal temperature loss (occurring for 
$T_d=80\,^{\circ }{\rm C}$ and 
$H=2.2$ m) is less than 
$2\,^{\circ }{\rm C}$. In the case of DMSO drop impacts, for which 
$T_d$ is very close to 
$T_a$, the estimated change is below 
$1\,^{\circ }{\rm C}$. To evaluate (B5), the following values for 
$T_a$ have been taken: 
$T_a=25\,^{\circ }{\rm C}$ in the case of water and 
$T_a=18\,^{\circ }{\rm C}$ for DMSO. Consequently, we choose to neglect this temperature change over the course of the present study for the sake of simplicity.
Appendix C. The three phase Stefan problem
The Stefan problem consists of solving the heat diffusion equation in several phases, in the presence of a moving interface of phase change (Rubinshteĭn Reference Rubinshteĭn1971). In the present study we are interested in the one-dimensional formulation of this problem involving either three phases for solidification (the liquid state, its corresponding ice and the solid substrate) or two for melting (the liquid state and its ice). The first situation is detailed below, and is strongly inspired by the self-similar analysis performed by Thiévenaz, Séon & Josserand (Reference Thiévenaz, Séon and Josserand2019). The model equations and boundary conditions are summarized in figure 6 for the three phase Stefan problem involving the liquid phase (droplet), the growing ice layer and the solid substrate. The properties for each phase are respectively denoted by the subscript letters 
$d$, 
$i$ and 
$s$. Assuming that advection can be neglected in the problem (see main text), the equations to be solved in the different phases are
$$\begin{gather} \partial_t T = \alpha_d \partial_z^2 T \quad z > h(t), \end{gather}$$
$$\begin{gather}\partial_t T = \alpha_i \partial_z^2 T \quad 0 < z < h(t), \end{gather}$$
$$\begin{gather}\partial_t T = \alpha_s \partial_z^2 T \quad z < 0, \end{gather}$$
where 
$T$ is the temperature field, and 
$\alpha _d$, 
$\alpha _i$ and 
$\alpha _s$ are the thermal diffusivities of the different phases (
$\alpha _n=k_n/(\rho _n c_{p,n})$, where subscript 
$n$ corresponds either to 
$s$, 
$i$ or 
$d$ depending on whether the substrate, ice or liquid is to be considered, respectively). In addition, the relevant boundary conditions to apply read
$$\begin{gather} \lim_{z \to +\infty} T = T_d, \end{gather}$$
$$\begin{gather}T(h^-) = T(h^+) = T_f, \end{gather}$$
$$\begin{gather}\rho_i \mathcal{L}_f {d}_t h = k_i \left. \partial_z T \right|_{h^-} - k_d \left. \partial_z T \right|_{h^+}, \end{gather}$$
$$\begin{gather}T(0^-) = T(0^+), \end{gather}$$
$$\begin{gather}k_s \left. \partial_z T \right|_{0^-} = k_i \left. \partial_z T \right|_{0^+}, \end{gather}$$
$$\begin{gather}\lim_{z \to -\infty} T = T_s, \end{gather}$$
with 
$T_d$ and 
$T_s$ the droplet and substrate temperatures, respectively, 
$T_f$ the melting point and 
$h$ the position of the moving liquid–ice interface. In addition, 
$k_s$, 
$k_i$ and 
$k_d$ stand for the thermal conductivities of the substrate, ice and liquid phase, respectively, 
$\rho _i$ is the ice density and 
$\mathcal {L}_f$ is the latent heat of fusion. In particular, condition (C5) ensures that the temperature at the liquid–ice interface corresponds to the melting point, (C8) establishes the equality of thermal fluxes at the ice–substrate interface (
$z=0$), and expression (C6) is the so-called Stefan condition. It expresses that the growth of the liquid–ice interface is dictated by the local difference of thermal fluxes (at 
$z=h$).
Figure 6. Model equations and boundary conditions for the three phase Stefan problem involving the liquid phase (expanding droplet), the growing ice layer and the solid substrate.
One may then apply the transformations 
$z \equiv H \, \bar {z}$, 
$h \equiv H \, \bar {h}$, 
$t \equiv (H^2/\alpha _i) \, \bar {t}$ and 
$T \equiv T_s + (T_f-T_s) \, \bar {T}$ to (C1)–(C3) and (C4)–(C9) to obtain the dimensionless problem
$$\begin{gather} \partial_{\bar{t}} \bar{T} = \omega_d \partial_{\bar{z}}^2 \bar{T} \quad \bar{z} > \bar{h}(\bar{t}), \end{gather}$$
$$\begin{gather}\partial_{\bar{t}} \bar{T} = \partial_{\bar{z}}^2 \bar{T} \quad 0 < \bar{z} < \bar{h}(\bar{t}), \end{gather}$$
$$\begin{gather}\partial_{\bar{t}} \bar{T} = \omega_s \partial_{\bar{z}}^2 \bar{T} \quad \bar{z} < 0, \end{gather}$$
with 
$\omega _d=\alpha _d/\alpha _i$, 
$\omega _s=\alpha _s/\alpha _i$, and with the associated boundary conditions
$$\begin{gather} \lim_{\bar{z} \to +\infty} \bar{T} = \bar{T}_d = \frac{T_d-T_s}{T_f-T_s}, \end{gather}$$
$$\begin{gather}\bar{T}(\bar{h}^-) = \bar{T}(\bar{h}^+) = 1, \end{gather}$$
$$\begin{gather}\frac{1}{St} {d}_{\bar{t}} \bar{h} = \left. \partial_{\bar{z}} \bar{T} \right|_{\bar{h}^-} - \left. \kappa \partial_{\bar{z}} \bar{T} \right|_{\bar{h}^+}, \end{gather}$$
$$\begin{gather}\bar{T}(0^-) = \bar{T}(0^+), \end{gather}$$
$$\begin{gather}\left.k_s \partial_{\bar{z}} \bar{T} \right|_{0^-} = \left. k_i \partial_{\bar{z}} \bar{T} \right|_{0^+}, \end{gather}$$
$$\begin{gather}\lim_{\bar{z} \to -\infty} \bar{T} = 0, \end{gather}$$
with 
${St}=c_{p,i} (T_f-T_s) / \mathcal {L}_f$ the Stefan number, whose definition involves the heat capacity 
$c_{p,i}=k_i/(\rho _i \alpha _i)$ of the ice, and 
$\kappa = k_d / k_i$. The Stefan number is a dimensionless quantity that compares the energy needed to decrease the ice temperature from 
$T_f$ to 
$T_s$ to the energy released by phase change, both energies being considered per unit mass.
It is possible to obtain a self-similar solution for (C10)–(C12) that satisfies boundary conditions (C13)–(C18) by conducting a similar analysis to that presented in Thiévenaz (Reference Thiévenaz2019). By doing so, a self-similar temperature field can be found only if it solely depends on the variable 
$\eta = \bar {z}/\sqrt {\bar {t}}$, whereas the liquid–ice front position 
$\bar {h}$ is proportional to 
$\sqrt {\bar {t}}$. This solution reads
$$\begin{gather} \bar{T} (\eta) = \bar{T}_d + \frac{\bar{T}_d-1}{1-\operatorname{erf} \left( \chi/(2 \sqrt{\omega_d}) \right)} \left( \operatorname{erf} \left( \eta/(2 \sqrt{\omega_d}) \right) - 1 \right) \quad\bar{z} > \bar{h}(\bar{t}), \end{gather}$$
$$\begin{gather}\bar{T} (\eta) = \frac{1}{r_i/r_s + \operatorname{erf} \left( \chi / 2 \right)} \left( \frac{r_i}{r_s} + \operatorname{erf} \left( \eta / 2 \right) \right) \quad 0 < \bar{z} < \bar{h}(\bar{t}), \end{gather}$$
$$\begin{gather}\bar{T} (\eta) = \frac{r_i/r_s}{r_i/r_s + \operatorname{erf} \left( \chi / 2 \right)} \left( 1 + \operatorname{erf} \left( \eta / \left( 2 \sqrt{\omega_s} \right) \right) \right) \quad \bar{z} < 0, \end{gather}$$
where 
$r_i=k_i/\sqrt {\alpha _i}$ and 
$r_s=k_s/\sqrt {\alpha _s}$ are the thermal effusivities of the ice and substrate, respectively, while 
$\chi$ is related to the position of the liquid–solid moving interface through 
$\bar {h}(\bar {t}) = \chi \sqrt {\bar {t}}$ and is the solution to the transcendental equation
\begin{equation} \frac{\chi \sqrt{\rm \pi}}{2 \,{St}} = \frac{{\rm e}^{-\chi^2/4}}{r_i/r_s+\operatorname{erf} \left( \chi/2 \right)} + \frac{r_d}{r_i} \frac{{\rm e}^{-\chi^2/(4 \omega_d)}}{1-\operatorname{erf} \left( \chi/(2\sqrt{\omega_d}) \right)} \frac{T_f-T_d}{T_f-T_s}, \end{equation}
with 
$r_d=k_d/\sqrt {\alpha _d}$ the thermal effusivity of water.
