The net energy effect of introducing a liquid layer between a solid (or vapor) and ice may be described by the Landau free energy Ω V W = A 0 A H / 12 π d 2 , (1) where A 0 is the surface area, and d the liquid layer thickness. We have introduced the Hamaker constant A H ∼ 1 0 − 20 − 1 0 − 19 J (albeit with an unconventional sign to keep A H positive).

Adding the free energy of the ice–water interface σ A ( r ) , where A ( r ) is the area of this interface and σ is the ice–water surface energy per unit area, to the energy of the pre-melted layer given in Equation 1 yields the total free energy Ω = σ A + A H A 0 12 π R − r 2 + Ω 0 , (2) where the bulk free energy Ω 0 is independent of the interface contributions. Because in general Ω = − P V , the combined potential for both the liquid and ice is Ω 0 = − P i V i − P w V 0 − V i , (3) where the ice pressure P i and water pressure P w will in general differ.

At the bulk melting temperature T m = 273 K, there will be no change in Ω 0 under a change in V i when μ and T are kept fixed, so, using the fact that Ω 0 = E − T S − μ N , we can write 0 = d Ω 0 = d E − T m d S − μ d N (4) where E , S , and N are the total ice–water energy, entropy, and molecule number, respectively. The heat needed to melt a volume − d V i of ice is ρ i λ d V i , where ρ i is the ice mass density and λ is the latent heat per unit mass. Using this, the above equation may also be written − ρ i λ d V i = T m d S = d E − μ d N . (5) Because ρ i and λ change very little over a modest temperature variation, we may also get d Ω 0 in the case where T ≠ T m by writing d Ω 0 = d E − μ d N − T T m T m d S ≈ − 1 − T T m ρ i λ d V i , (6) which indicates that the free energy change due to an ice volume increase is negative below the bulk freezing point. Integrating from V i = 0 , where Ω o = − P w ( μ , T ) V 0 , yields Ω 0 ≈ − 1 − T T m ρ i λ V i − P w V 0 , (7) which, when inserted in Equation 2, gives Ω = σ A + A H A 0 12 π R − r 2 − 1 − T T m ρ i λ V i − P w V 0 . (8) The change in this energy as r is increased from r = 0 is Δ Ω = σ A − 1 − T T m ρ i λ V i + A H 12 π A 0 R − r 2 − A 0 R 2 . (9)

The equilibrium value of the ice radius is given by the global minimum of Δ Ω , which, for sufficiently small R values, will be at r = 0 , that is, for the complete liquid state. Above this critical pore size, the minimum will be at r m ≲ R , a value that is given by the equilibrium thickness of the surface melted layer. As may be noted from Figure 3, this minimum does not change much with the pore size R . When Ω ( r = 0 ) > Ω ( r m ) , there will still be ice. The condition for complete melting is that Ω ( r = 0 ) < Ω ( r m ) , which yields the freezing-point depression.

Taking r = R , A H = 0 gives the free energy change in passing from a liquid to a fully frozen pore Δ Ω = σ A − 1 − T T m ρ i λ V i , (10) where A = 4 π R 2 and V i = ( 4 / 3 ) π R 3 . The condition Δ Ω = 0 implies that a pore of radius R will freeze at a temperature T = T F given by T F R = T m 1 − 3 σ ρ i λ R . (11) This is the standard expression for the Gibbs–Thomson effect. For a cylindrical pore, the geometrical factor of 3 must be replaced by 2. In the following, we shall use the value 3. Note that, due to the tendency of the surface tension to minimize the interface area, these smooth geometrical shapes will also be relevant in more complex pore geometries.

Thus far, we have ignored the time it takes for a metastable state to be replaced by the equilibrium state, implicitly assuming that the system has had time to reach the overall minimum state for the free energy. This is in general not the case as some metastable states may be very long-lived, a phenomenon that is quantified in classical nucleation theory [17, 18], which is based on the probability that a free energy barrier is traversed by the thermal activation energy k B T . Moreover, the stability against melting may be very different from the stability against the reverse process of freezing. It is generally much more difficult to superheat a solid than to supercool a liquid [12, 13]. Superheated crystalline solids have only been observed in some rather singular cases where the heated region is along a single crystal plane or the crystals are confined inside a non-melting matrix [19, 20]. The existence of supercooled liquids, on the other hand, only requires the absence of nucleation sites.

Assuming our pre-melted surface layer of water, there is no extra energy cost (nucleation barrier) in forming a new liquid–ice surface during the melting process. Yet, there will be a nucleation barrier that must be crossed during melting when the temperature is T ≈ T F . The reason for this is that when melting happens around T ≈ T F < T m , the bulk free energy must increase, while the surface energy is decreased. Thus, as r decreases from a value around R in Figure 3, the free energy initially increases. As a result, there is a free energy barrier against both melting and freezing. This fact implies the possible existence of solid ice that is superheated relative to its depressed freezing point T F .

Using nucleation theory, it is possible to estimate the lifetime of these metastable states as ∝ exp ( ( Δ Ω ( r max ) − Δ Ω ( R ) ) / ( k B T ) ) , where the free energy Δ Ω ( r ) is shown in Figure 3 and takes its maximum at r = r max . Requiring that the lifetime be within a realistic range, it is possible to show that the melting temperature must be increased above T F by an amount that corresponds to a reduction of the freezing-point depression T m − T F by ∼ 20% for pore sizes greater than ∼ 1 nm.

It may be shown that nucleation barriers are significantly more influential during freezing (supercooled liquid). In this case, however, nucleation pathways other than ice forming as a spherical crystal are likely to dominate, as has been shown for the case where ice nucleates in pockets or corner geometries [3, 4, 11].

In the following, we will consider melting on the basis that metastable states may be neglected, although there is a nucleation barrier to be passed both for the melting and freezing transition in isolated pores. For melting, this assumption implies that there may be quantitative corrections to the depression T m − T F by ∼ 20%, which are ignored.

Having dealt with the equilibrium problem of the freezing-point depression, we now investigate the non-equilibrium effects of this phenomenon in the context of a nanoporous material. We shall consider a melting front, for which the shift in melting temperature is small, and so the shift in the freezing-point depression will not be applied. Note, however, that a freezing front may differ significantly from the melting front through the possible existence of metastable pockets of supercooled liquid.

The pore size distributions may be estimated through nitrogen absorption [21], electron microscopy, or mercury injection experiments and measurements of the heat capacity variations with temperature when there is water present [22]. For silts, clays, and synthetic media made of glass powders [1], they may yield distributions that extend down at least to the nm scale. Freezing and melting of water confined in silica nanopores have been observed down to pore sizes of 3 nm [6].

The distributions may be given in terms of a relative volume fraction per unit length g ( R ) so that ∫ d R g ( R ) = ϕ , the porosity of the medium. Our main assumption is that this distribution may be approximated with a power law above a minimum cut-off length R min , g R = N R − R min β (12) where N = ( β + 1 ) ϕ / ( R max − R min ) β + 1 is the normalization.

Mercury intrusion experiments are challenged by the fact that high injection pressures may crush or deform the smallest pores. Yet, in rigid materials, such as cement, the technique may be used to measure pores down to R min ∼ 1 nm [23]. In order to cover the smaller pore ranges, nitrogen adsorption techniques are often better [21]. Zhao et al. [24] measured pore size distribution for porous sandstone from the Ordos basin by mercury injection, finding g ( R ) -distributions that are well described by R min ∼ 10 nm and β = 1 over 1 to 2 decades in pore sizes. Using N 2 adsorption techniques on porous glass powders, Fujinomori soil, and bentonite clay, Watanabe et al. [1] found g ( R ) -distributions where R min ≈ 1 nm, R max ≈ 3–4 nm, and β = 1–2. Park et al. [25] measured pore-size distributions. Different sediments produced R min values from 1 nm to 100 nm with distributions that could be described by a β ≈ 2 power law over roughly a decade. There is thus a range of natural and synthetic materials that seem to fulfill the assumed pore-size power law distribution over an adequate range of length scales.

In a medium that is described by Equation 12, all the pores are frozen when T = T min = T m − T m r 0 R min , (13) where we have introduced the length r 0 = 3 σ / ( ρ i λ ) ≈ 3.3 nm. Correspondingly, there is an upper temperature T max = T m − T m r 0 R max , (14) where all pore water is melted.

The initial filling fraction of water in the pores C 0 gives the total water (ice or liquid) fraction ϕ C 0 . The fraction of liquid water, C ( T ) , is the fraction contained in the pores that are so small that they have not frozen. These pores have sizes less than r T = 3 σ ρ i λ 1 − T / T m . (15) This means that when T min < T < T max C T = C 0 ∫ R min r T d R g R = N C 0 ∫ R min r T d R R − R min β = N C 0 β + 1 r T − R min β + 1 (16) = N C 0 β + 1 R min β + 1 T − T min T m − T β + 1 , (17)

by use of Equation 12 and Equation 13. Close to the absolute freezing point T min , the above denominator is close to ( T m − T min ) β + 1 , so we shall use C T = B T − T min T min β + 1 , (18) with B = C 0 ϕ R min 2 R max − R min ρ i λ T min 3 σ T m β + 1 (19)

by use of Equation 13 and the definition of N .

Having neglected the thickness of the pre-melted films in the ice-filled pores by setting A H = 0 , we should compare the relative contributions to C ( T ) from these films and the liquid-filled pores. Because there is no film in the liquid-filled pores, we need only take the R > r ( T ) pores into account. We take the film contribution to be given by the film thickness d = R − r ( T ) as Δ C T = C 0 ∫ r T R max d R g R Δ V R V R , (20) where the fraction Δ V / V ≈ 3   d / R is the ratio of the film volume to the pore volume. Then, Δ C T = 3 d N C 0 ∫ r T R max d R R − R min β R , (21) where R max is the upper cut-off for g ( R ) . When β = 1 , this integral is easily evaluated to give Δ C T = 3 d T N C 0 R max − r T − R min ⁡ ln R max r T . (22)

Taking the R max term to dominate in this expression and using Equation 16, the ratio becomes Δ C C ≈ 6 d T R max r T − R min 2 , (23) which may well be larger than one when r ( T ) ≳ R min .

However, as we shall see below, it is the rates of change d C / d T of the volume fractions that are important, not the absolute value of Δ C / C . The film thickness d may be estimated from Equation 9 as the minimum of Δ Ω when σ = 0 . This gives the standard expression [16]. d = A H T m 6 π T m − T ρ i λ 1 / 3 , (24)

where we can use the relatively high value A H = 1 0 − 19 J. Together with the constants given in Table 1, this gives d ≈ 1 nm, while the other relevant length, which appears in Equation 9, is 3 σ / ( ρ i λ ) ≈ 0.3 nm.

Because C ∝ ( T − T min ) 2 and, to leading order Δ C ∝ d ∝ ( T − T m ) − 1 / 3 , we have that d C / d T = ( β + 1 ) C / ( T − T min ) while d Δ C / d T = Δ C / ( 3 ( T m − T ) ) , so that the ratio of the changes in these two quantities due to a temperature change when β = 1 is δ Δ C δ C = T m T − T min d r 0 R max R min 3 . (25) close to the absolute freezing point T min . Here, we have used Equation 13 to substitute ( T m − T min ) / T m = r 0 / R min . The condition that δ Δ C / δ C ≪ 1 may be taken as a condition on the range of pore sizes R max / R min : R max R min ≪ R min 2 d r 0 T − T min T m , (26) or, equivalently, T − T min T m − T min ≫ d R min R max R min . (27) When R min = 30 nm, for instance, and ( T − T min ) / T m = 1 / 273, we get that R max / R min ≪ 10 . It is quite natural that the condition for the domination of pore-versus film fluid is a limited range of pore sizes, as a domination of the large pores, which all carry a film contribution, would leave a smaller fraction of the porosity to be represented by smaller pores.

In other words, when R max / R min ≪ 10 , Δ C changes relatively slowly with T ≈ T min , but C changes significantly. In this case, we may neglect the variations in the film contribution to the overall change in liquid volume fraction. For this reason, we shall only use the C in the following, keeping in mind that it is the fraction of liquid pore water, and not the total fraction of liquid water.

In a 1D setting, the conservation of energy in a slab of thickness d x over a time d t may be written j x − j x + d x Δ A d t = λ ρ Δ A d x d C + c b Δ A d x d T (28) where Δ A is the cross-sectional area, c b is the combined specific heat capacity of the porous medium and the water, and T is the temperature. In Equation 28, the left-hand side is the net energy transfer to the slab, the first term on the right is the energy consumed by melting (latent heat), and the last term is the energy absorbed due to the heat capacities of the water, ice, and the porous medium itself. As d x → 0 , Equation 28 becomes ∂ j ∂ x + λ ρ ∂ C ∂ t + c b ∂ T ∂ t = 0 . (29) To describe the heat flow, we apply the Fourier law, which takes the form j = − κ ∂ T ∂ x , (30) where κ is the bulk thermal conductivity of the porous medium, so inserting this in Equation 29 gives ∂ ∂ x κ ∂ T ∂ x = ρ λ ∂ C ∂ T + c b ∂ T ∂ t (31) where we have used ∂ C / ∂ t = ( ∂ C / ∂ T ) ( ∂ T / ∂ t ) . Generalizing to arbitrary dimension ( ∂ / ∂ x → ∇ ) and replacing T − T min by T min ( C / B ) 1 / ( β + 1 ) using Equation 18 yields the diffusion equation 1 + M ∂ C ∂ t = D 0 ∇ ⋅ C 1 / β + 1 − 1 ∇ C . (32) where D 0 = κ b T min λ ρ β + 1 B 1 / β + 1 = c b T min λ ρ β + 1 B 1 / β + 1 D t , (33) and D t = κ b / c b ≈ 1 mm 2 / s is the average thermal diffusivity of the porous medium, and M = c b T min C 1 / β + 1 − 1 β + 1 ρ i λ B 1 / β + 1 = 3 c b σ T m β + 1 ρ i λ 2 C 0 ϕ 1 / β + 1 R max − R min R min 2 C − γ . (34)

The last expression comes from replacing B by the expression in Equation 19, and γ = β / ( β + 1 ) . Using the material constants in Table 1 gives M = l R max − R min R min 2 C − γ (35) where l ≈ 0.4 nm. So M ≪ 1 when the range of pore sizes is limited and as long as C does not become too small.

We shall proceed to analyze the case where the M -term may be dropped, leaving the equation ∂ C ∂ t = D 0 ∇ ⋅ C − γ ∇ C . (36)

We note at this point that the condition for neglecting the energy needed to change temperature, which is represented by the M -term, coincides with the condition to neglect the contribution of pre-melted films. Both conditions may be fulfilled by media with a limited range of pore sizes above a minimum size R min ≳ 10 nm.

The fact that we have neglected the energy contribution given by the heat capacities means that we have assumed that all the energy is spent melting the ice in the pores. We note in passing that the same assumption is made in treatments of the moving boundary problem associated with melting fronts (the Stefan problem) [26].

The mobile energy density ∝ C , which means that Equation 36 may be read as a statement of energy conservation. We will consider the response to a localized addition of energy that causes a local initial volume V i of melted water. Solving Equation 36 subject to the normalization condition V i = ∫ d V C (37) in d dimensions [27] for a point source initial C ( r , 0 ) gives C r , t = p r / t τ f d t , (38) where f t = 2 − d γ 1 − γ V i γ D 0 t 1 2 − d γ . (39) and p y = γ 2 1 − γ V i γ y 2 + k − 1 γ , (40) where y = r / f ( t ) and k is an integration constant given by the normalization condition. It takes the value of [27]. k = V i 1 − d γ / 2 γ 2 π 1 − γ d 2 Γ 1 γ Γ 1 γ − d 2 2 γ d γ − 2 . (41)

The functional form given in Equation 38 immediately yields the second moment for the concentration profile r rms 2 = ∫ d r r d − 1 r 2 C r , t ∫ d r r d − 1 C r , t ∝ f 2 t ∝ t 2 τ (42) with [27, 28]. Here τ = 1 2 − d γ , (43) or, in terms of β , τ = β + 1 d − d − 2 β + 1 . (44) When the dimension d = 3 , we obtain hyper-ballistic spreading ( τ > 1 ) if 1 / 2 < β < 2 and superdiffusion ( τ > 1 / 2 ) if 0 < β < 2 . A value of β = 1 , which would correspond to a linear initial growth of g ( R ) , gives β = 1 and τ = 2 , which should be compared to the ballistic value τ = 1 and the normal diffusive value of τ = 1 / 2 . A heat pulse will thus spread at an accelerating rate, causing a rapid, long-range front of melting.

The d = 1 case is relevant when a heat pulse spreads downward in the ground. In this case, β = 1 gives τ = 2 / 3 , and β = 2 gives τ = 3 / 4 , both superdiffusive values.

The superdiffusive spread of C and thus T with time follows from the fact that a heat pulse will lower the local melting temperature, thus keeping the remaining ice from receiving more latent heat as the temperature is rising. In contrast, a melting front that propagates through a medium with a single pore size will only spread diffusively, propagating at a speed ∼ t .