Confinement induced control of similarity solutions in premelting dynamics and other thin film problems
Satyajit Pramanik John S. Wettlaufer
physics.flu-dyn















[]satyajit.math16@gmail.com
Nordita, Royal Institute of Technology and Stockholm University, Stockholm, Sweden




[]john.wettlaufer@yale.edu
Yale University, New Haven, Connecticut, USA
Mathematical Institute, University of Oxford, Oxford, UK
Nordita, Royal Institute of Technology and Stockholm University, Stockholm, Sweden












We study the combined effects of nonlocal elasticity and confinement induced ordering on the dynamics of thermomolecular pressure gradient driven premelted films bound by an elastic membrane. The confinement induced ordering is modeled using a film thickness dependent viscosity. When there is no confinement induced ordering, we recover the similarity solution for the evolution of the elastic membrane, which exhibits an infinite sequence of oscillations. However, when the confinement induced viscosity is comparable to the bulk viscosity, the numerical solutions of the full system reveal the conditions under which the oscillations and similarity solutions vanish. Implications of our results for general thermomechanical dynamics, frost heave observations and cryogenic cell preservation are discussed. 
Finally, through its influence on the viscosity, the confinement effect implicitly introduces a new universal length scale into the volume flux.  Thus, there are a host of thin film problems, from droplet breakup to wetting/dewetting dynamics, whose properties (similarity solutions, regularization, and compact support) will change under the action of the confinement effect. Therefore, our study suggests revisiting the mathematical structure and experimental implications of a wide range of problems within the framework of the confinement effect.








Confinement induced control of similarity solutions in premelting dynamics and other thin film problems
    John S. Wettlaufer
    December 30, 2023
=======================================================================================================








§ INTRODUCTION
 
The free surfaces of a wide range of materials (e.g., rare gas solids, metals, quantum and molecular solids) are observed to have a thermodynamically stable thin liquid film that disjoins the solid from its vapor phase, another grain of solid, or another material entirely. This extension of the equilibrium domain of the melt phase into the solid region of the bulk phase diagram is known as interfacial premelting, the detailed manifestation of which depends on the dominant polarization and surface forces in the system <cit.>. The same general phenomena, whereby volume-volume interactions that are extensive with the system volume compete with those that are extensive with interfacial area, are at play in, amongst many other settings, wetting, biomaterials, ceramics, colloids, and tribology <cit.>. Here we focus on premelting dynamics, in a setting in which the melt phase appears between the solid and a flexible wall and flows from high to low temperatures under the influence of thermomolecular pressure gradients <cit.>. These thin film flows have been studied theoretically and experimentally on ice <cit.>, argon <cit.>, benzene <cit.> and helium <cit.>. Importantly, because the liquid phase is adjacent to its own solid, the region near the latter is endowed with solid-like ordering over a length scale that is system specific. It is this “confinement-induced” ordering that we examine in the context of premelting dynamics, extending a recent study in which the membrane imposed a local elastic hoop stress <cit.>. 

A host of effects are responsible for confined liquids exhibiting properties that differ from their bulk counterparts <cit.>. For example, molecular dynamics simulations of confined nanometric liquids show wall normal fluid layering as a function of temperature, the nature of the wall-fluid and the fluid-fluid interactions, and, for temperatures close to the melting point, transitions between shear thinning and shear thickening are observed <cit.>. Interfacially premelted water on ice surfaces exhibits mobility attenuation several orders of magnitude larger than that of bulk water <cit.>, which is of relevance, for example to the kinetic friction of ice surfaces <cit.>. Recent experiments on the dynamics of premelted water on ice surfaces determined the ratio of surface tension-to-shear viscosity, from which it was inferred that the shear viscosity of premelted water is ten times that of bulk supercooled water <cit.>. This observation is a central motivation of our work. 

Experiments show that the shear viscosity of thin liquid films of different materials are considerably greater than their bulk values <cit.>. For example, the effective viscosity of nanometric thin water films confined between different surfaces increases from the bulk value by almost seven orders of magnitude <cit.>. The increase in effective viscosity of molecularly thin films of octamethylcyclotetrasiloxane  as the films thinned was shown to depend on film thickness as a power law <cit.>. Confinement-induced enhancement of the coefficient of shear viscosity has also been reported for dodecane <cit.> and alkanes <cit.>, and depends upon the quenching rate and the degree of perfection of the mica surfaces. Contrary to these systems, aqueous salt solutions can evidently retain their bulk properties at nanometer scales <cit.>. 

Despite the host of implications of the confinement effects of liquids on the dynamical consequences of premelted liquid films, there has been a single study along these lines. In a recent paper, motivated by the experiments discussed above, we treated the premelted liquid viscosity as a power law function of the film thickness <cit.>. We studied the influence of confinement effects in premelting dynamics on ice surfaces in elastic capillary tubes with a local hoop stress coupling, a setting that had previously been examined without confinement effects <cit.>. 

Our purpose here is to explore the coupled effects of nonlocal elastic deformation of a membrane confining a thin premelted film on the surface of ice, and confinement enhanced viscosity, in a geometry studied previously in the absence of confinement effects <cit.>. In the absence of confinement effects one finds a family of fourth order parabolic partial differential equations (PDEs) that depend on the nature of the interactions responsible for the existence of the premelted film, and there is an associated family of ordinary differential equations (ODEs) with a family of similarity solutions <cit.>. Here, amongst other things, we find that there exists a finite range of the power-law exponents of the confinement induced viscosity model that impose a length scale on the problem and thereby lead to the breakdown of the similarity solutions.  This has general consequences for a wide range of thin film problems.

Of particular relevance are the consequences for continued flow of liquid at low temperatures where confinement effects should be at play and where no quantitative experiments have been undertaken. Settings where this may play out include the flow of unfrozen water in terrestrial, planetary and engineering environments <cit.>, and during biological freezing <cit.>. For example, cardiological treatment is often facilitated by lowering temperatures to lower metabolic rates and during cryogenic preservation cells are preserved below freezing point of water <cit.>. Thus, mechanical damage to cells during freeze-thaw cycling is a fundamental challenge <cit.>, for which many approaches are taken. For example, optimal cryogenic cell preservation may be facilitated by introduction of cryoprotectants, such as ice-binding proteins <cit.>, and controlled rates of solidification <cit.>, all of which will be influenced by the slow flow of confined premelted water abutting compliant walls. 

Finally, we note that similar higher-order diffusion equations govern the dynamics of viscous currents down an incline <cit.>, and a spate of other thin-film settings <cit.> including models of Hele-Shaw flows <cit.>, viscous gravity currents between a rigid surface and a deformable elastic sheet <cit.> and many others too numerous to allow for an exhaustive list. Motivated by these problems, there is a literature on the existence, regularity, and the finite time blow up of nonlinear, degenerate, fourth-order parabolic differential equations (e.g., <cit.> and references therein). Mathematically our study falls under the general rubric of this class of equations, but physically it is a unique thin film system, being driven as it is by thermomolecular pressure gradients and phase changes.  Importantly, because of the new length scale introduced by the confinement effect, its inclusion will influence all of these problems in a fundamental manner. 

In Section <ref>, we briefly review the model and the non-dimensional formulation. In Section <ref>, we discuss the similarity solutions in the two limiting cases of viscosity dominiated by either the bulk or the confinement induced values. Section <ref> examines the membrane deformation for different parameter values, and the numerical results are compared with experiments in Section <ref>, before concluding. 

 



§ PREMELTING DYNAMICS AND THERMODYNAMICS
 
In Figure <ref> we show a schematic of the experimental setting that motivates our study. The interface between a single ice crystal and a flexible polymer membrane is subject to a temperature gradient parallel to the interface <cit.>. The membrane/ice interface is disjoined by a premelted layer of water, which flows from high to low temperature. Because at each temperature the premelted liquid film has a thermodynamically required thickness, determined by Eq. (<ref>) below, mass conservation insures that a gradient in the liquid volume flux will lead to solidification and hence the deformation of the height of the membrane, h̃(x̃, t̃), measured relative to its initial reference height. Our choice of coordinates places the bulk solid-liquid interface at x̃ = 0, where x̃ increases with decreasing temperature. Because the deformation of the membrane height is much smaller than the disk radius, a one dimensional model is appropriate and was first developed to analyze experiments using the bulk liquid viscosity <cit.>. We revisit this model, but with a viscosity that depends on the thickness of the premelted film. 

The free energy of the ice/film/membrane system is given by the grand potential 

    Ω = -P_l V_l - P_s V_s + ℐ (d̃),

where ℐ(d̃) = [ Δσ f(d̃) + σ_sw] A_i is the effective interfacial free energy, 
d̃ is the film thickness, 
A_i is the interfacial area, and P and V are the pressure and volume.  The difference between the dry and wet interfacial free energies is Δσ = σ_sl + σ_lw - σ_sw,  where the σ's are the solid-liquid (sl), liquid-wall (lw), and solid-wall (lw) values, where σ_sw > σ_sl + σ_lw is a sufficient condition for the interface to premelt <cit.>. When the polarization forces underlying the wetting behavior are power law interactions, then f(d̃) is written as 

    f(d̃) = 1 - (  d_0 /d̃)^n-1,

where n depends on the nature of interactions attracting the liquid to the solid and d_0 in the order of molecular diameter. For example, non-retarded (retarded) van der Waals interactions are represented by n = 3 (n = 4). We minimize Ω at fixed temperature T, close to the bulk melting temperature T_m, total volume and chemical potential to obtain two key results. Firstly <cit.>,  

    d̃ = [ -  (n-1) d_0^n-1Δσ/ρ_l q_m ]^1/n t_r^-1/n≡λ_n t_r^-1/n,

where, t_r = (T_m - T)/T_m is the reduced temperature, ρ_l is the liquid density, and q_m the latent heat of fusion. Secondly, the pressure difference between the liquid film and the bulk solid, 

    P_l - P_s = Δσ (n - 1) d_0^n-1d̃^-n,

which reflects the polarization forces attracting the liquid to the solid. Combining the temperature dependent film thickness (<ref>) and the disjoining pressure (<ref>) we obtain the melt pressure as a function of the temperature 

    P_l = P_s - ρ_l q_m t_r,

independent of the nature of interactions, characterized by n. Clearly, when fixing the solid pressure, P_s, and imposing a temperature gradient parallel to the premelted surface, a pressure gradient drives thin film flow towards lower temperatures <cit.>, as shown in Figure <ref>.

As the temperature decreases equation (<ref>) shows that the film thins. However, as the temperature becomes very low, this does not capture the fact that depending upon crystallographic orientation of the solid the film properties are intermediate between a liquid and a solid (for the ice-water interface, see <cit.>). Therefore, the ordering effect of the solid on the liquid may have the dynamical consequences as the premelted film continually thins at low temperatures. This proximity effect on the liquid film viscosity is modeled as  

    η̃( d̃ ) = η_0 - ( η_0 - η_b )[ 1 - (  d_0 /d̃)^γ],

where γ > 0 such that η̃( d̃ ) →η_b (bulk viscosity) as d̃→∞, and η̃( d̃ ) →η_0 when d̃→ d_0^+, and the continuum approximation breaks down <cit.>. For d̃ = λ_n t_r^-1/n and a constant temperature gradient of 𝒢, the reduced temperature depends on position through

    t_r = 𝒢x̃/T_m.

Thus, the viscosity in (<ref>) can be expressed as a function of the position as

    η̃( x̃ ) = η_0 - ( η_0 - η_b ) [ 1 - { d_0 /λ_n ( 𝒢/ T_m )^1/n}^γx̃^γ/n],

and, as we shall see, this position dependence creates an internal length scale in the problem that leads to a parameter dependent breakdown in the family of similarity solutions that are found when the viscosity is constant <cit.>. 

In the lubrication approximation the volume flux Q, per unit breadth through the film of thickness d̃ given by (<ref>), is

    Q(x̃) = - d̃^3 / 12 η̃(x̃) ∂_x̃ P_l.

The pressure P_s exerted on the solid phase by the membrane is proportional to the membrane curvature κ, which in the small slope approximation is κ≈h̃_x̃x̃, giving 

    P_s = -σ̃h̃_x̃x̃,

where the membrane surface tension is σ̃. Combining (<ref>), (<ref>), and (<ref>) – (<ref>), mass conservation is written as

    ∂_t̃h̃ + Γ_n ∂_x̃[ x̃^-3/n/η̃(x̃)/η_b ( 1 + α̃h̃_x̃x̃x̃) ] = 0,

where 

    Γ_n = λ_n^3 ρ_l q_m /12 η_b(T_m/𝒢)^(3/n-1)  and  α̃ = σT_m/ρ_s q_m 𝒢,

Because the temperature gradient 𝒢 is a constant, at each position x̃ thermodynamic equilibrium requires the film thickness to satisfy (<ref>). Therefore, when the thermomolecular pressure driven flow brings liquid in excess of the equilibrium value, the evolution equation (<ref>) describes how the membrane must deform in response to the solidification of the premelted fluid, as seen schematically in Figure <ref>.



























The quantitative properties of the premelting dynamics with the proximity effect of (<ref>) differ from the bulk viscosity case, but not in the high temperature regime that has thus far been experimentally accessible <cit.>. Hence, our calculations here are intended to motivate experiments at low temperatures and long time scales where the proximity effect becomes important. 



 §.§ Non-dimensionalization

We non-dimensionalize the system using  the following scaling; 


    x = x̃/ X_0 ,    t = t̃/ X_0^2+3/n /Γ_n,    h = h̃/ X_0 , 
       η( x ) = η_0 - η̃/η_0 - η_b  = [ 1 - { d_0 /λ_n ( 𝒢 X_0 / T_m )^1/n}^γ x^γ/n] 
                   ≡[ 1 - d_γ, n x^γ/n],     where  
    d_γ, n = {d_0/λ_n(𝒢 X_0/T_m)^1/n}^γ,


from which we obtain 


    ∂_t h + ∂_x [ x^-3/n/η_1(x)( 1 + α h_xxx) ] = 0, 
       η_1(x) = 1 + (η_r - 1)d_γ,n x^γ/n,


where α = α̃/X_0^2 and η_r = η_0/η_b. 
Here, X_0 is a length scale chosen based on the magnitude of the temperature gradient 𝒢, which determines the temperature drop across the size of the sample (Figure <ref>).  Thus, for a sufficiently large experimental cell, the membrane height at the lowest temperatures (the cell center) will be unaffected by the deformation at the highest temperature, where x=0 <cit.>. This, along with appeal to experiment, leads to the following boundary conditions: 


    ∂_x h = ∂_xx h = 0      (x = 0),     and
       ∂_x h = 0      (x = 1).



The boundary conditions
(<ref>) correspond to the bulk solid-liquid interface, where the slope and curvature of the membrane height vanish. 

Equation (<ref>) is a “no flux condition” for the membrane height at the cold end, far from the warm end, but note there is a fluid flux at the cold end. The initial condition for the membrane height, which is measured relative to a reference height, is 

    h(x, t = 0) = 0.
 
To avoid a singularity at the bulk solid-liquid interface, we impose the following necessary condition on the solution

    ∂_xxx h = -α^-1     (x = 0).
 

In order to obtain a unique membrane height we need to specify the fourth boundary condition. In the original theory, described in <ref> below, wherein the bulk liquid viscosity is used, the resulting family of PDEs (namely (<ref>) with a constant viscosity) admits a similarity solution, leading to a family of ODEs and associated boundary conditions; (<ref>)–(<ref>) <cit.>. Analysis of the evolution equation (<ref>) far from the origin provides the fourth boundary condition, (<ref>), which is satisfied numerically using a shooting method. Thus, in that study there was no need to solve the PDEs for h(x,t), but the introduction of the proximity effect here effectively imposes an “external” length scale (at x in the interior of the domain) that leads to a parameter dependent breakdown of the similarity solution. Hence, here we must solve the family of PDEs, for which we specify the curvature at the cold end (x = 1). This differs from using an asymptotic analysis to determine the boundary condition, but as discussed in Appendix <ref> the influence of this boundary condition is only important for parameter values for which the magnitude of the oscillations in the membrane are large at very low temperatures. The experiments done thus far <cit.> have been under conditions wherein the majority of the membrane deformation occurs near the warm end. Thus, at the cold end we impose 

    ∂_xx h = 0      (x = 1).

We note, however, for analysis of future experiments across a variety of materials this condition will depend on both the control parameters as well as experimental geometry. 

Analysis of (<ref>) shows that the liquid viscosity across the film is dominated by the bulk value for the power law exponent γ≥γ_ sat and confinement-induced value for γ≤γ_ crit, whereas the bulk and confinement-induced values are comparable for γ_ crit≤γ≤γ_ sat <cit.>. 

The dimensionless viscosity of the film is η_1(x) = 1 + (η_r - 1)d_γ,n x^γ/n≡ 1 + ϵ(γ) x^γ/n.
For the relevant parameter values, one can show that ϵ(γ) is a decreasing function of γ.  Hence, γ_crit (γ_sat) is the largest (smallest) approximate value of γ such that near the warm (cold) end the confinement-induced value of the viscosity of the thin liquid film, ϵ(γ) x^γ/n, is of the same order of magnitude as that of its bulk value η_b.
See Appendix <ref> for estimates of γ_crit and γ_sat.






We rewrite the mass conservation equation (<ref>) as 

    ∂_t h = 𝒟(x) h_xxxx + 𝒰(x)h_xxx + ℱ(x),

where 


    𝒟(x) = -αx^-3/n/η_1(x) < 0, 
       𝒰(x) =  d𝒟/ dx = αx^-(1+3/n)/η_1^2(x)(3+γ)η_1(x) - γ/n > 0,   and
       ℱ(x) = 𝒰(x) /α = x^-(1+3/n)/η_1^2(x)(3+γ)η_1(x) - γ/n > 0,


which is the form solved numerically using standard second order finite difference formulae and the time integration is performed using a semi-implicit Crank-Nicolson (C-N) method <cit.>. We extend the method and stability analysis of the numerical scheme discussed in <cit.> for a second order equation to the fourth order evolution equation (<ref>). We have used a uniform grid spacing Δ x = x_j+1 - x_j, ∀ j ∈{ -2, -1 , 0, ,N, N+1, N+2 } for the spatial variable, such that

    x_-2 < x_-1 < (0 =) x_0 < x_1 <  < x_N-1 < x_N (= 1) < x_N+1 < x_N+2,

and the time stepping is so chosen that the numerical stability is met <cit.>. (See Appendix <ref> for grid independence test.) 














§ SIMILARITY SOLUTIONS
 
For the bulk viscosity case, η_r = 1, a family of similarity solutions, described by a single parameter n, is obtained as


    g(ζ) = h̃α̃( Γ_n α̃t̃)^-3n/(4n+3), 
       ζ = x̃( Γ_n α̃t̃)^-n/(4n+3),


that satisfies the following fourth order family of ordinary differential equations <cit.>: 


    g^''''   = 3/ng^''' + 1/ζ + n/4n+3ζ^(3/n)+1g^' - 3n/4n+3ζ^3/ng, 
    
     g^'   = g^'' = 0      (ζ = 0), 
     
    g^'   = 0      (ζ = ζ_ c), 
     
    g    ∼ζ^-[(3/n)+1]     (ζ→∞),   and
    
    g^'''   = -1      (ζ = 0),


where, ζ_c denotes the cell center. This similarity solution is valid for both η_r = 1 <cit.> and for parameter dependent values of the viscosity exponent from (<ref>), γ≥γ_ sat,  heuristically like the case for a second order system that treats a related problem in a different geometry <cit.>.  

On the other hand, when the viscosity exponent is γ≤γ_ crit <cit.>, we obtain a two parameter (n,  γ) family of similarity solutions,  

    g_γ(ζ_γ) = h̃α̃[ Γ_n α̃t̃/δ]^-3n/(4n+3+γ),

with the similarity variable 

    ζ_γ = x̃[ Γ_n α̃t̃/δ]^-n/(4n+3+γ),

where 

    δ = ( η_r - 1 ) [ d_0/λ_n( 𝒢/T_m)^1/n]^γ,

that are the solutions to 

    g_γ^'''' = 3+γ/ng_γ^''' + 1/ζ_γ + n/4n+3+γζ_γ^(3+n+γ)/ng_γ^' - 3n/4n+3+γζ_γ^(3+γ)/ng_γ,

subject to the conditions, (<ref>), (<ref>), and (<ref>), written in terms of the appropriate dependent (g_γ) and independent (ζ_γ) variables. As ζ_γ→∞, asymptotic analysis of 

(<ref>) reveals that the membrane height decays algebraically as 

    g_γ∼ζ_γ^-(3+n+γ)/n     (ζ_γ→∞).




The system (<ref>)-(<ref>) (or, the version corresponding to the two parameter family of similarity solutions g_γ(ζ_γ)) is solved numerically using a finite difference method with uniform mesh. To validate our numerical method, we reproduce figure 2 of <cit.>, shown here in Figure <ref>(a). Figure <ref>(b) compares five members of the second family of similarity solutions, g_γ(ζ_γ). 



§ RESULTS
 
Here, we investigate the time evolution of the membrane height, described by the family of fourth order parabolic differential equations (<ref>), whose solutions exhibit oscillatory behavior as is common in a wide class of fourth order equations appearing in many settings (See e.g., <cit.>, and Refs. therein).
In this section we present the numerical results for (a) temperature gradients of 𝒢 = 0.92 × 10,   0.92 × 10^2 and 0.92 × 10^3 K/m, (b) molecular diameters of d_0 = 1,   3,   5,   8 and 10 Å, and (c) a wide range of the power law exponent γ in connection with η_r. Unless otherwise stated, the other physical parameters used are listed in Table <ref>. 







 §.§ Effects of the temperature gradient 𝒢
 

Here we discuss the effect of the temperature gradient on the evolution of the membrane height. Note that both the dimensionless parameters α and  d_γ, n depend on 𝒢. For the case of bulk viscosity model, premelting dynamics is independent of d_γ, n; the only dimensionless parameter that captures the effect of the temperature gradient is α∝𝒢^-1. Figure <ref> shows the membrane height, h(x, t), for the bulk viscosity case and a large (α = 0.0669) and small (α = 0.6686) temperature gradient.  As α decreases, deformation is localized near the warm end, but the qualitative properties of the deformation are unaffected. Namely, as α decreases, because d̃^3 ∝ t_r^-3/n∝( α^-1x)^-3/n, the film thins rapidly towards the lower temperatures. For example, in the case of nonretarded van der Waals interactions (n = 3), a decrease in α by a factor of ten, reduces the film thickness more than 50% everywhere in the domain. Thus, mass conservation dictates that more premelted liquid freezes at higher temperatures, as seen in Figure <ref>.

Similar qualitative effects of a steady temperature gradient on the deformation of an elastic capillary wall with a local hoop stress are seen <cit.>, because the main physical effects are at play; as the temperature gradient increases, the gradient in the liquid volume flux is localized at higher temperatures.






 §.§ Effects of short range cut-off, d_0
 
We of course know the size of a molecule, but depending on the crystallographic orientation, the ordering effect on the liquid differs <cit.>.  Hence, we view the short range cut-off, d_0, as a proxy for crystallographic orientation and its ordering influence. Recall from (<ref>) that 

    d_γ, n = {d_0/λ_n(𝒢 X_0/T_m)^1/n}^γ,

and hence for fixed values of λ_n,  𝒢,   X_0,   T_m and the exponents n and γ, we see that as d_0 increases both d_γ, n and η_1(x) increase (see Eq. (<ref>)). Thus, viscous resistance (volume flux) in the thin liquid film increases (decreases) more rapidly parallel to the temperature gradient as d_0, and hence the ordering influence, increases.  Therefore, the solidification and concomitant deformation, is localized near the warm end, qualitatively like increasing the temperature gradient. For the bulk viscosity model (η_r = 1) premelting  dynamics is independent of d_0. However, d_0 has a significant influence on the solidification in the proximity viscosity model. 

Figure <ref> shows the membrane evolution, h(x, t), for η_r = 10^5, α = 0.6686, and different values of d_0. For each d_0, we compute γ_ crit, γ_ sat and use γ_ crit < γ < γ_ sat in figure <ref>. The corresponding values of the dimensionless parameter d_γ, n are d_6, 3 = 1.2188 × 10^-5, 1.6719 × 10^-2. 






 §.§ The influence of the proximity effect
 
In this section we discuss the influence of the power law exponent, γ, which governs the strength of the decay of the proximity effect, on the evolution of the membrane. We vary γ over a wide range and study the numerical solutions for η_r = 10^3, 10^5 and 10^7. We recover the similarity solutions (<ref>) and (<ref>) when γ≤γ_ crit and γ≥γ_ sat, respectively. For a given value of α this effect is captured via d_γ, n. Here we consider α = 0.0669 and non retarded van der Waals interaction (n = 3).

In the bulk viscosity model the principal positive curvature high temperature maximum leads to an oscillation as follows. Firstly, membrane elasticity acts to suppress the curvature, the response to which is to draw in more premelted fluid. Secondly, the fluid drawn in from the low temperature region creates a negative curvature that elasticity acts to suppress driving fluid towards both the low and high temperature directions. Finally, the fluid driven towards the low temperature direction creates a new peak, and the process repeats itself <cit.>. Here, the same basic oscillation in the membrane height appears when γ≤γ_ crit, and so we focus  on the membrane evolution for γ∈ (γ_ crit, γ_ sat). 





In Figures <ref> we show the membrane evolution for η_r = 10^7, α = 0.0669 and d_γ, 3 = 4.51 × 10^-2  (γ = 5.1), 1.1031 × 10^-4  (γ = 15). Comparing Figure <ref>(a) with Figure <ref>(a) we see that the maximum deformation of the membrane is one order of magnitude larger for the bulk viscosity case than that for γ = 5.1 ≈γ_ crit. On the other hand, for γ = 19 < γ_ sat, the maximum deformation is nearly the same as that for the bulk viscosity model, which is consistent with the analysis used to predict γ_ sat for a similar problem <cit.>. However, for γ = 15, in Figure <ref>(b) we observe two differences from the bulk viscosity model. Firstly, the maximum membrane deformation h_ max is larger than in the bulk viscosity case and secondly, membrane oscillations are suppressed. The interaction between the elasticity induced curvature and the proximity effect on the dynamics of the membrane height can be seen by examining the membrane evolution in a region surrounding the first minimum and the second maximum.  This is shown for d_γ, 3 = 3.7175 × 10^-4  (γ = 13) in Figure <ref>, with the other parameter values the same as in Figure <ref>. It is seen that the oscillations vanish for a parameter dependent time interval; the same general behavior is found for a wide range of temperature gradients (α), molecular diameters (d_0) and viscosity ratios (η_r), as discussed for γ∈ (γ_ crit, γ_ sat) in Section <ref> and Appendix <ref>. 























 


  §.§.§ Reentrant oscillations of h
 
As discussed in section <ref>, the existence of an infinite number of peaks in the membrane height with an exponentially decaying amplitude is found for the case of bulk viscosity (e.g., Figure <ref>), and such oscillations are influenced by the confinement effects modeled by (<ref>). 



Thus, when the confinement-induced and bulk viscosity are comparable, the associated internal length scale imposed on the system influences the membrane oscillations (Figures <ref>(b) and <ref>). 

For such oscillations to persist, the membrane curvature must change sign, and as seen in Figures <ref>(b) and <ref>, the film thickness dependent viscosity is responsible for the vanishing of the first trough and hence the oscillations. The evolution equation (<ref>) shows how the flux operator, written as 𝒟(x) in (<ref>), controls the dynamics of the membrane height. The viscosity model (Eq. (<ref>)) shows that the proximity effects, η_1 = 1 + (η_r - 1) d_γ, n x^γ/n∼ (η_r - 1) d_γ, n x^γ/n, dominate near the cold end (x = 1), whereas the bulk viscosity effects, η_1 = 1 + (η_r - 1) d_γ, n x^γ/n∼ 1, dominate near the warm end (x = 0). In terms of the flux operator, at the cold end we have -𝒟(x) ∼ x^-(3+γ)/n and at the warm end we have -𝒟(x) ∼ x^-3/n, and thus the similarity solution in the bulk viscosity region grows faster, as h ∼ t^3n/(4n + 3), than the solution in the proximity-viscosity dominated region, for which h ∼ t^3n/(4n + 3 + γ). In the crossover between these two scaling behaviors  the membrane height dynamics is intermediate and depends on η_r and d_γ, n (i.e., γ,   n,  λ_n,   d_0,  𝒢,   T_m). Namely, for γ∈ (γ_ crit, γ_ sat), -𝒟(x) exhibits both ∼ x^-3/n and ∼ x^-(3+γ)/n scalings locally within the spatial domain, giving different temporal scalings for the evolution of h(x,t). 

Consider the situation wherein the first minimum of h(x,t) is located in the bulk viscosity-dominated region, but the second maximum is located within the crossover region or beyond. In such a situation, the first minimum grows faster than the second maximum, which is overcome by the minimum until at a time t^†, h(x, t^†) decreases monotonically with x. At later times, when the deformation of the membrane enters the proximity viscosity-dominated region, the oscillatory structure in h(x, t) reappears. When circumstances are such that there is a rapid change in the viscosity across the first minimum the fluid viscosity increases rapidly with x.  However, as elasticity reduces the positive curvature created by the primary maximum, it draws fluid in from both the high and low temperatures regions. However, because the viscous friction in the proximity dominated low temperature region is larger than that in the bulk viscosity high temperature region, more fluid is drawn from the higher temperature lower x region. Therefore, the magnitude of the negative curvature in the low temperatures end is reduced, suppressing the the first minimum. Eventually, the primary peak spans both the bulk viscosity and crossover regions and continued oscillations can be described along similar lines.  





§ DISCUSSION AND CONCLUSIONS
 
The key dynamics in the class of conservation laws that our study is a part of are embodied in the structure of the flux term.  The key distinctions between our work and the majority of such conservation laws governed by one or another manifestation of lubrication theory are that (1) our physical system has a thin film whose thickness that is a known function of position, but the evolution equation is for the elastic surface with which it is in contact and (2) we deal with two families of conservation laws, depending on the nature of the interactions controlling the film thickness and those influencing the ordering of the film as it thins.   

The structure of the flux is given by 𝒟(x), as expressed by (<ref>), and is plotted in Figure <ref> along with the limits for (i) the bulk dominated viscosity, ∝ x^-3/n, and (ii) the proximity dominated viscosity, ∝ x^-(3 + γ)/n. We examine 𝒟(x) with η_r = 10^7, and small,  intermediate, and large temperature gradients in Figures <ref>(a)–(c) for different values of d_γ, n. For example, for α = 0.0669 and d_γ, n = 3.7175 × 10^-4  (γ = 13, n = 3) we have



    - 𝒟(x) ∼ x^-3/n,   for  0 ≤ x ≲ 0.09   (≡ x_l), 
           - 𝒟(x) ∼ x^-(3+γ)/n,   for  (x_r ≡)   0.25 ≲ x ≤ 1,

and for the intermediate range, x_l < x < x_r, the flux decays with x (see Figure <ref>(a)). Based upon 
these flux scaling behaviors, combined with the results from Figures <ref> and <ref> (also Figure <ref>), we make the following observations. 

Firstly, for two points x_1 and x_2 such that x_1 ∈ (0, x_l) and x_2 ∈ (x_l, x_r), h(x_1, t) grows faster than h(x_2, t). Now, if h_1^0 ≡ h(x_1, t_0) < h(x_2, t_0) ≡ h_2^0 at time t = t_0, then the differential growth of h(x,t) at these two points at some later time t_f > t_0 is such that h_1^f = h(x_1, t_f) > h(x_2, t_f) = h_2^f. This also holds true for x_3 ∈ (x_l, x_r) and x_4 ∈ (x_r, 1). 

Secondly, when the primary high temperature maximum of h spreads over 0 ≤ x ≤ x^†, where x^† > x_r, then, provided x^† is at a sufficiently warm position, the second crest and first trough reappear. 

Thirdly, there are times t_1 and t_2 such that ∀ t ∈ (t_1, t_2), h(x, t) is a monotonic function of x.  (We have not obtained analytical approximations of t_1 and t_2.) For fixed values of 

η_r and α,
the crossover region moves towards the cold (warm) end as γ increases (decreases) in the interval (γ_ crit, γ_ sat). Therefore, both t_1 and t_2 increase with γ∈ (γ_ crit, γ_ sat). 


For example, Figure <ref>(b) shows that the crossover region covers almost 70% of the domain for α = 0.6686 and d_9, 3 = 4.2224 × 10^-6. Therefore, in this case we anticipate a large time window (t_1, t_2). Whereas for α = 6.6861 and d_7, 3 = 3.0658 × 10^-7, the crossover region spans more than 80% of the domain. In this case, a second crest located within the crossover region reappears for a short intermediate time interval, showing that the existence and the length of the time interval for h(x, t) to be a monotonic function of x depends nonlinearly on d_γ, n. On the other hand, for α = 6.6861 and d_4, 3 = 1.8968 × 10^-4 (and the same parameters as above), the crossover region is confined to a boundary layer near the warm end (x ≲ 0.01, Figure <ref>(c)). Therefore, we expect oscillations in h(x, t) for ostensibly all t > 0. (Note that ( t_2 - t_1 ) should be compared for different γ∈ (γ_ crit, γ_ sat) for fixed values of α, η_r, and n as the prefactor of the similarity solutions depends on all of these parameters.) 





Finally, we compare the numerical solutions of our model to experiments <cit.>. We choose d_0 = 3 Å (the molecular diameter of water is 2.75 Å) and experiments on nanometric water films confined between different surfaces show that η_r = 10^7 <cit.>. Thus, using these parameters and the experimental temperature gradient and membrane tension <cit.>, we perform numerical simulations with different γ and both van der Waals (n = 3) and electrostatic interactions (n = 2). For van der Waals interactions, the agreement between theory and experiment is best with γ = 5.9, whereas for the electrostatic interactions the agreement is best with γ = 6.7, capturing nearly the entire profile with an overshoot to the maximum deformation for the case of electrostatic interactions (Figure <ref>). Overall, whereas one might view there to be a better agreement with the data relative to the bulk viscosity model (see figure 4 of <cit.>), here we have an additional parameter, γ. Next, we compare the numerical solutions with the measurements of the maximum membrane deformation. As shown in Figure <ref>(b), a least squares fit of the experimental data <cit.> agrees well with our numerical solution with a superior agreement for van der Waals interactions (note again the overshoot in the case of electrostatic interactions), although both are within experimental error. Importantly, we have provided a framework to experimentally test the proximity effect in other materials and to examine the power law treatment of the proximity effect embodied in the exponent γ. 

Our results have implications for frost heave in the freezing of soils, cells and tissues, cancer treatment, food science, fish biology and botany  <cit.>, amongst other phase change settings. In particular, our results show that deformation continues even in very low temperatures regions or if the liquid domain is highly restricted, such as in the case of confined bodies subject to temperature gradients, which may rupture (e.g., see <cit.> and references therein). Moreover, if the temperature gradient is very shallow and the liquid film is thin, deformation can still be large, and our framework provides constraints for designing optimal biological cryogenic preservation <cit.>. Finally, among many other biological applications, the processes discussed here are closely related to phase changes and the ordering of water in membranes, which are active areas of inquiry <cit.>. 




 §.§ Generalizations and Implications
 
We have shown here that for γ_crit⩽γ⩽γ_sat the mobility term in the flux operator is not scale invariant.  In consequence, because the mobility term depends on the position the similarity solution structure breaks down.  This influence of the proximity effect has implication for a wide class of higher order diffusion-like equations arising in thin film flows, which we can write as

    H_t + [ 1/η(H) F (H, H_x, H_xx, H_xxx) ]_x = 0,

where F is a real-valued, nonlinear function of the film thickness H(x, t) and its spatial derivatives up to third oder.   Clearly, the functional form of F and the film-thickness dependent effective viscosity η(H) determine the similarity solution structure, as we have shown here as a function of γ.
We write the confinement effect, as given in equation (<ref>), as
η(H) = η_b + (η_0 - η_b) ( H_0/H)^γ, where H_0 is a structural length scale of the fluid.  Therefore, 


  * For η(H) ≈η_b, there exists an outer solution to equation (<ref>) away from the region H ≈ H_0 whenever 

    (η_0 - η_b) ( H_0/H)^γ≪η_b,         H/H_0≫ (η_r -1 )^1/γ > 1.


  * For η(H) ≈ (η_0 - η_b) ( H_0/H)^γ, there exists an inner solution to equation (<ref>) near the region H ≈ H_0 whenever 

    (η_0 - η_b) ( H_0/H)^γ≫η_b,          1 < H/H_0≪ (η_r -1 )^1/γ.


  * The transition from the outer solution to the inner solution occurs on a length scale satisfying 

    H/H_0∼ (η_r -1 )^1/γ⩾ 1         η_r ⩾ 2,   ∀γ > 0.


This length scale is universal for any thin film problem embodied by (<ref>) under the influence of this proximity effect.  Thus, this introduces new length scales to problems such as for example droplet breakup, contact line dynamics and wetting and dewetting wherein the finite-time singularity structure of the flux operator will be regularized at larger than traditional (e.g., van der Waals) length scales <cit.>, as has been seen in recent experiments on moving contact lines in complex fluids <cit.>.    


A few recent examples of relevance include the inclusion of interface effects in a continuum-based study of transport phenomena in a microfluidic channel by Vo and Kim <cit.>, a classical unsteady Stokes flow analysis of oscillating sphere near a wall by Fouxon and Leshansky <cit.>, which would differ at close range where the confinement effect may be operative, and the study of droplet breakup by Deblais et al. <cit.>.  In this latter combined experimental and numerical study, the authors show inclusion of a small amount of viscosity is sufficient to induce breakdown of the universality of the breakup dynamics predicted by inviscid theory.  Whereas as a general proposition, this is not surprising (namely inviscid theory should in some limit break down), it is a compelling question to examine the relevant length scales of break down and to test the concepts of confinement described here by systematically changing the fluid structure. 





































§ LOW TEMPERATURE BOUNDARY CONDITION 

As discussed in <ref> the original theory wherein the bulk viscosity was used the system of ODEs that resulted from the similarity solution given by (<ref>)–(<ref>) was solved <cit.>. Analysis of the evolution equation (<ref>) at low temperatures far from the origin lead to the fourth boundary condition, (<ref>), which is satisfied by a shooting method when solving the system numerically. Here we imposed (<ref>). In Figure <ref> we show that there is no change in the membrane height near the warm end and minor changes in the membrane height near the cold end. 






§ APPROXIMATION OF Γ_CRIT AND Γ_SAT

Here, we briefly explain the estimates of the transition values of the parameter γ; γ_crit and γ_sat. 
The value of γ_crit is associated with ϵ(γ_crit) ≫ 1, such that η̃(x̃)/η_b ≈ 1 for 0 ⩽ x ⩽ε (≪ 1) and η̃(x̃)/η_b ∼𝐎 (x^γ/n) for ε < x ⩽ 1.  Thus, 

    (η_r - 1) [ d_0/λ_n( 𝒢 X_0/T_m)^1/nε^1/n]^γ_crit≈ 1,     so that
          ln(η_r - 1) + γ_crit[ ln( d_0/λ_n) + 1/nln( 𝒢 X_0/T_m) + 1/nlnε] ≈ln 1 = 0,     and hence
    
          γ_crit≈ - ln(η_r - 1) /ln( d_0/λ_n) + 1/n[ ln( 𝒢 X_0/T_m) + lnε]  .

The value of γ_sat is associated with ϵ(γ_sat) ≈ 1, viz., 

    (η_r - 1) [ d_0/λ_n( 𝒢 X_0/T_m)^1/n]^γ_sat≈ 1,     so that
          ln(η_r - 1) + γ_sat[ ln( d_0/λ_n) + 1/nln( 𝒢 X_0/T_m) ] ≈ln 1 = 0,     and hence
    
          γ_sat≈ - ln(η_r - 1) /ln( d_0/λ_n) + 1/nln( 𝒢 X_0/T_m) .
 

















































































§ CODE VALIDATION
 



We use the five-point stencil formulae for the boundary conditions (<ref>), (<ref>), and (<ref>) and obtain the following algebraic equations: 


    h_-2 - 8 h_-1 + 8 h_1 - h_2 = 0, 
     
    -     h_-2 +16 h_-1 - 30 h_0 + 16 h_1 - h_2 = 0, 
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
        h_N-2 - 8 h_N-1 + 8 h_N+1 - h_N+1 = 0, 
     
    -     h_N-2 + 16 h_N-1 - 30 h_N + 16 h_N+1 - h_N+2 = 0.


The constraint of avoiding the singularity at x = 0 is checked at every time step using the second order centered difference approximation of the necessary condition (<ref>)

    h_2 - 2h_1 + 2h_-1 - h_-2 = -2 Δ x^3 α^-1.







We validate our code by reproducing results from Ref. <cit.>. We solve the partial differential equations corresponding to their constant bulk viscosity model using standard second order finite difference formulae and the time integration is performed using a semi-implicit Crank-Nicolson (C-N) method <cit.>. Additionally, we compute the corresponding similarity solution from the numerically computed membrane deformation, h(x,t), using (<ref>) and (<ref>), as well as solving the corresponding family of ordinary differential equations (<ref>)-(<ref>). The solutions obtained from these two methods are plotted in Figure <ref>, thereby validating our C-N method. 

 

We define the spatial and temporal discretization errors as 

    ‖ E_h(t) ‖ = ( ∑_i = 1^N| h_Δ x(x_i,t) - h_0.5 Δ x(x_2i,t) |^2 )^1/2,   and
          ‖ E_h(x) ‖ = ( ∑_i = 1^N| h_Δ t_1(x_i,n_1Δ t_1) - h_Δ t_2(x_i,n_2 Δ t_2) |^2 )^1/2,

respectively. Here, h_Δ x represents the numerical solution obtained with a spatial grid size Δ x = 1/N, whereas h_0.5 Δ x is the numerical solution corresponding to the spatial grid size 0.5 Δ x, and n_1, n_2 ∈ℕ, such that n_1 Δ t_1 = n_2 Δ t_2. 
Numerical results are spatially (temporally) grid independent if, ‖ E_h(t) ‖ < ϵ_1 uniformly for all t > 0 and a specified tolerance ϵ_1 > 0 (‖ E_h(x) ‖ < ϵ_2, uniformly for all x ∈Ω and a specified tolerance ϵ_2 > 0). For a convergent numerical solution, the optimal spatial and temporal discretization are obtained for N = 512 (or 1024 depending upon various parameters, e.g., η_r,   d_0,  λ_n,  γ,  𝒢) and a dimensional value Δ t = 10 s, such that the global discretization error, 

    ϵ(Δ x, Δ t) = ‖ E_h(t) ‖ + ‖ E_h(x) ‖ ,

is minimized. Figure <ref> shows grid independence of Δ x and Δ t. 



§ ADDITIONAL RESULTS OF REENTRANT OSCILLATIONS OF H 
 
Here we present additional results of parameter dependent reentrant oscillations of the membrane height h̃(x̃, t̃). 

































§ ACKNOWLEDGMENTS

The authors acknowledge the support of the Swedish Research Council Grant No. 638-2013-9243, conversations with their colleagues in Stockholm and New Haven, and J.S.W. acknowledges a Royal Society Wolfson Research Merit Award for support. 













58
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Wilen et al.(1995)Wilen,
  Wettlaufer, Elbaum, and Schick]Wilenetal1995
  author author L. A. Wilen, author J. S. Wettlaufer, author M. Elbaum,
   and author M. Schick,  @noop  journal journal Phys. Rev. B  volume 52, pages 12,426 (year 1995)NoStop
[Dash et al.(2006)Dash,
  Rempel, and Wettlaufer]Dash2006
  author author J. G. Dash, author A. W. Rempel,   and author J. S. Wettlaufer,  @noop  journal journal Rev. Mod. Phys.  volume 78 (year 2006)NoStop
[Wettlaufer and Worster(2006)]Wettlaufer2006
  author author J. S. Wettlaufer and author M. G. Worster, @noop  journal journal Annu.
  Rev. Fl. Mech. volume 38, pages 427
  (year 2006)NoStop
[Schick(1990)]SchickLesHouches
  author author M. Schick, in @noop  booktitle Les Houches
  Session XLVIII (publisher Elsevier, address
  Amsterdam, year 1990) pp. pages
  415–497NoStop
[Safran(1994)]Safranbook
  author author S. A. Safran, @noop  title Statistical
  Thermodynamics of Surfaces, Interfaces, and Membranes (publisher Addison-Wesley, address Reading, Massachusetts,  year 1994)NoStop
[French et al.(2010)French,
  Parsegian, Podgornik, Rajter,
  Jagota, Luo, Asthagiri,
  Chaudhury, Chiang, Granick,
  Kalinin, Kardar, Kjellander,
  Langreth, Lewis, Lustig,
  Wesolowski, Wettlaufer, Ching, Finnis, Houlihan, von
  Lilienfeld, van Oss, and Zemb]French:2010
  author author R. H. French, author V. A. Parsegian, author R. Podgornik, author R. F. Rajter, author A. Jagota,
  author J. Luo, author
  D. Asthagiri, author
  M. K. Chaudhury, author
  Y.-m. Chiang, author
  S. Granick, author S. Kalinin, author M. Kardar, author R. Kjellander, author D. C.  Langreth, author J. Lewis, author S. Lustig,
  author D. Wesolowski, author J. S. Wettlaufer, author W.-Y. Ching, author
  M. Finnis, author F. Houlihan, author O. A.  von Lilienfeld, author C. J.  van Oss,  and author T. Zemb, 10.1103/RevModPhys.82.1887
  journal journal Rev. Mod. Phys. volume 82, pages 1887 (year
  2010)NoStop
[Israelachvili(2011)]Jacobbook
  author author J. N. Israelachvili, @noop  title Intermolecular
  and Surface Forces, edition 3rd ed. (publisher
  Academic Press, address New York, NY, year
  2011)NoStop
[Wilen and Dash(1995)]Wilen1995
  author author L. A. Wilen and author J. G. Dash, 10.1103/PhysRevLett.74.5076 journal
  journal Phys. Rev. Lett. volume 74,  pages 5076 (year 1995)NoStop
[Wettlaufer and Worster(1995)]Wettlaufer1995
  author author J. S. Wettlaufer and author M. G. Worster, @noop  journal journal Phys.
  Rev. E volume 51, pages 4679
  (year 1995)NoStop
[Wettlaufer et al.(1996)Wettlaufer, Worster, Wilen, and  Dash]Wettlaufer1996
  author author J. S. Wettlaufer, author M. G. Worster, author L. A. Wilen,
   and author J. G. Dash,  @noop  journal journal Phys. Rev. Lett.  volume 76, pages 3602 (year
  1996)NoStop
[Rempel et al.(2004)Rempel,
  Wettlaufer, and Worster]Rempel2004
  author author A. W. Rempel, author J. S. Wettlaufer,  and author M. G. Worster, @noop  journal journal J.
  Fluid Mech. volume 498, pages 227
  (year 2004)NoStop
[Zhu et al.(2000)Zhu,
  Vilches, Dash, Sing, and  Wettlaufer]Zhu2000
  author author D.-M. Zhu, author O. E. Vilches,
  author J. G. Dash, author B. Sing,  and author
  J. S. Wettlaufer, 10.1103/PhysRevLett.85.4908 journal journal
  Phys. Rev. Lett. volume 85, pages
  4908 (year 2000)NoStop
[Taber(1929)]Taber1929
  author author S. Taber, @noop  journal journal J.
  Geol. volume 37, pages 428 (year 1929)NoStop
[Mizusaki and Hiroi(1995)]Mizusaki1995
  author author T. Mizusaki and author M. Hiroi, @noop  journal journal Physica
  B volume 210, pages 403 (year 1995)NoStop
[Pramanik and Wettlaufer(2017)]Pramanik2017
  author author S. Pramanik and author J. S. Wettlaufer, 10.1103/PhysRevE.96.052801 journal journal Phys. Rev. E volume
  96, pages 052801 (year 2017)NoStop
[Evans(1990)]Evans1990
  author author R. Evans, @noop  journal journal J.
  Phys. Cond. Matt. volume 2, pages
  8989 (year 1990)NoStop
[Christenson(2001)]Hugo2001
  author author H. K. Christenson, @noop  journal journal
  J. Phys. Cond. Matt. volume 13, pages
  R95 (year 2001)NoStop
[Koplik and Banavar(1995)]Koplik1995
  author author J. Koplik and author J. R. Banavar, @noop  journal journal Annu.
  Rev. Fl. Mech. volume 27, pages 257
  (year 1995)NoStop
[Pittenger et al.(2001)Pittenger, Fain, Cochran, Donev, Robertson, Szuchmacher, and  Overney]Pittenger2001
  author author B. Pittenger, author S. C. Fain,
  author M. J. Cochran, author J. M. K. Donev, author
  B. E. Robertson, author
  A. Szuchmacher,  and author
  R. M. Overney, @noop  journal journal Phys. Rev. B volume
  63, pages 134102 (year 2001)NoStop
[Döppenschmidt et al.(1998)Döppenschmidt, Kappl, and Butt]Butt1998
  author author A. Döppenschmidt, author M. Kappl,  and author H.-J. Butt, @noop  journal journal J. Phys.
  Chem. B volume 102, pages 7813
  (year 1998)NoStop
[Dash(2003)]Dash2003
  author author J. G. Dash, @noop  journal journal Scripta
  Mat. volume 49, pages 1003 (year 2003)NoStop
[Murata et al.(2015)Murata,
  Asakawa, Nagashima, Furukawa, and Sazaki]Murata2015
  author author K. Murata, author H. Asakawa,
  author K. Nagashima, author Y. Furukawa,  and author G. Sazaki, @noop 
  journal journal Phys. Rev. Lett. volume 115, pages 256103 (year
  2015)NoStop
[Granick(1991)]Granick1991
  author author S. Granick, @noop  journal journal
  Science volume 253, pages 1374
  (year 1991)NoStop
[Dhinojwala and Granick(1997)]Dhinojwala1997
  author author A. Dhinojwala and author S. Granick, 10.1021/ja9632318 journal
  journal J. Am. Chem. Soc. volume
  119, pages 241 (year 1997)NoStop
[Raviv and Klein(2002)]Raviv2002
  author author U. Raviv and author J. Klein,  @noop  journal journal Science volume 297, pages 1540 (year
  2002)NoStop
[Zhu and Granick(2003)]Zhu2003
  author author Y. Zhu and author S. Granick,  10.1021/la035155+ journal journal Langmuir volume 19, pages
  8148 (year 2003)NoStop
[Zhu and Granick(2004)]Zhu2004
  author author Y. Zhu and author S. Granick,  @noop  journal journal Phys. Rev. Lett.  volume 93, pages 096101 (year 2004)NoStop
[Major et al.(2006)Major,
  Houston, McGrath, Siepmann,  and Zhu]Major2006
  author author R. C. Major, author J. E. Houston,
  author M. J. McGrath, author J. I. Siepmann,  and author X.-Y. Zhu, @noop 
  journal journal Phys. Rev. Lett. volume 96, pages 177803 (year
  2006)NoStop
[Bureau(2010)]Bureau2010
  author author L. Bureau, 10.1103/PhysRevLett.104.218302 journal journal Phys. Rev. Lett. volume 104, pages 218302 (year
  2010)NoStop
[Kessler and Werner(2003)]Kessler2003
  author author M. A. Kessler and author B. T. Werner, 10.1126/science.1077309 journal
  journal Science volume 299, pages 380 (year 2003)NoStop
[Andersland and Ladanyi(2004)]Andersland2004
  author author O. B. Andersland and author B. Ladanyi, @noop  title An Indroduction to
  Frozen Ground Engineering (publisher Chapman and Hall,  address New York, NY, year 2004)NoStop
[Rempel(2007)]Rempel2007
  author author A. W. Rempel, 10.1029/2006JF000525 journal
  journal Journal of Geophysical Research: Earth Surface  volume 112, pages F02S21 (year 2007)NoStop
[Peterson and Krantz(2008)]Peterson2008
  author author R. A. Peterson and author W. B. Krantz, @noop  journal journal
  Journal of Geophysical Research: Biogeosciences volume 113, pages G03S04 (year
  2008)NoStop
[Tamppari et al.(2012)Tamppari, Anderson, Archer, Douglas, Kounaves, McKay, Ming, Moore, Quinn, Smith,
  Stroble, and Zent]Tamppari:2012
  author author L. K. Tamppari, author R. M. Anderson, author J. Archer,
  P. D., author S. Douglas, author S. P. Kounaves, author C. P. McKay,
  author D. W. Ming, author Q. Moore, author
  J. E. Quinn, author
  P. H. Smith, author
  S. Stroble,  and author
  A. P. Zent, @noop  journal journal Antarct. Sci. volume
  24, pages 211 (year 2012)NoStop
[Rubinsky(2003)]Rubinsky2003
  author author B. Rubinsky, 10.1023/A:1024734003814 journal
  journal Heart Failure Reviews volume
  8, pages 277 (year 2003)NoStop
[Hawthorn and Rolfe(2016)]Hawthorn2016
  author author J. Hawthorn and author E. Rolfe, @noop  title Low Temperature Biology
  of Foodstuffs: Recent Advances in Food Science, Recent advances in food
  science (publisher Elsevier Science, year
  2016)NoStop
[Ishiguro and Rubinsky(1994)]Ishiguro1994
  author author H. Ishiguro and author B. Rubinsky, @noop  journal journal
  Cryobiology volume 31, pages 483 
  (year 1994)NoStop
[Bar-Dolev et al.(2012)Bar-Dolev, Celik, Wettlaufer, Davies, and Braslavsky]Bar-Dolev:2012
  author author M. Bar-Dolev, author Y. Celik,
  author J. S. Wettlaufer,
  author P. L. Davies,  and  author I. Braslavsky, 10.1098/rsif.2012.0388 journal journal
  J. R. Soc. Interface volume 9, pages
  3249 (year 2012)NoStop
[Rubinsky and Pegg(1988)]Rubinsky1988
  author author B. Rubinsky and author D. E. Pegg, 10.1098/rspb.1988.0053 journal
  journal Proceedings of the Royal Society of London B:
  Biological Sciences volume 234, pages
  343 (year 1988)NoStop
[Huppert(1982)]Huppert1982
  author author H. E. Huppert, @noop  journal journal
  Nature volume 300, pages 427
  (year 1982)NoStop
[Troian et al.(1989)Troian,
  Herbolzheimer, Safran, and Joanny]Troian1989
  author author S. M. Troian, author E. Herbolzheimer, author S. A. Safran,  and author J. F. Joanny, @noop  journal journal EPL
  (Europhysics Letters) volume 10, pages
  25 (year 1989)NoStop
[Oron et al.(1997)Oron,
  Davis, and Bankoff]Oron:1997
  author author A. Oron, author S. H. Davis,   and author S. G. Bankoff,  @noop  journal journal Rev. Mod. Phys.  volume 69, pages 931 (year
  1997)NoStop
[Craster and Matar(2009)]Craster:2009
  author author R. V. Craster and author O. K. Matar, 10.1103/RevModPhys.81.1131 journal
  journal Rev. Mod. Phys. volume 81,  pages 1131 (year 2009)NoStop
[Constantin et al.(1993)Constantin, Dupont, Goldstein,
  Kadanoff, Shelley, and Zhou]Constantin:1993
  author author P. Constantin, author T. F. Dupont, author R. E. Goldstein, author L. P. Kadanoff, author M. J. Shelley,  and author S. M. Zhou, @noop  journal journal Phys.
  Rev. E volume 47, pages 4169
  (year 1993)NoStop
[Hewitt et al.(2015)Hewitt,
  Balmforth, and De Bruyn]Hewitt2015
  author author I. J. Hewitt, author N. J. Balmforth,  and author J. R. De Bruyn, @noop  journal journal
  European Journal of Applied Mathematics volume
  26, pages 1 (year 2015)NoStop
[Bertozzi and Pugh(2000)]Bertozzi2000
  author author A. L. Bertozzi and author M. C. Pugh, @noop  journal journal Indiana
  University Mathematics Journal volume 49, pages 1323 (year 2000)NoStop
[Nada and Furukawa(1997)]Nada:1997
  author author H. Nada and author Y. Furukawa, @noop  journal journal J.
  Phys. Chem. B volume 101, pages 6163
  (year 1997)NoStop
[Nada and Furukawa(2000)]Nada:2000
  author author H. Nada and author Y. Furukawa, @noop  journal journal
  Surf. Sci. volume 446, pages 1
  (year 2000)NoStop
[Nada(2016)]Nada:2016
  author author H. Nada, 10.1063/1.4973000 journal journal J. Chem. Phys. volume 145 (year 2016), 10.1063/1.4973000NoStop
[LeVeque(2007)]LeVeque2007a
  author author R. LeVeque, @noop  title Finite Difference
  Methods for Ordinary and Partial Differential Equations (publisher Society for Industrial and Applied Mathematics, year 2007)NoStop
[Bertozzi(1998)]Bertozzi:1998
  author author A. L. Bertozzi, @noop  journal journal Not.
  Am. Math. Soc. volume 45, pages 689
  (year 1998)NoStop
[Chapman et al.(2013)Chapman, Trinh, and Witelski]Chapman:2013fk
  author author S. J. Chapman, author P. H. Trinh,
   and author T. P. Witelski,  10.1137/120872012 journal journal SIAM J. Appl. Math. volume 73, pages 232 (year 2013)NoStop
[Vlahou and Worster(2010)]Vlahou2010
  author author I. Vlahou and author M. G. Worster, @noop  journal journal J.
  Glaciol. volume 56, pages 271
  (year 2010)NoStop
[Parsegian(2005)]Parsegian2005
  author author V. A. Parsegian, @noop  title Van der Waals Forces:
  A Handbook for Biologists, Chemists, Engineers, and Physicists (publisher Cambridge University Press, year
  2005)NoStop
[Eggers(2018)]EggersPRFluids2018
  author author J. Eggers, 10.1103/PhysRevFluids.3.110503 journal journal Phys. Rev. Fluids volume 3, pages 110503 (year
  2018)NoStop
[Nanjundiah et al.(2019)Nanjundiah, Kurian, Kaur, Singla, and Dhinojwala]ContactLine
  author author K. Nanjundiah, author A. Kurian,
  author S. Kaur, author
  S. Singla,  and author
  A. Dhinojwala, @noop  journal journal Phys. Rev. Lett. volume 122 (year 2019)NoStop
[Vo and Kim(2016)]VoSciRep2016
  author author T. Q. Vo and author B. Kim,  10.1038/srep33881 journal journal Sci. Rep. volume 6, pages
  33881 (year 2016)NoStop
[Fouxon and Leshansky(2018)]FouxonPhysRevE2018
  author author I. Fouxon and author A. Leshansky, 10.1103/PhysRevE.98.063108 journal journal Phys. Rev. E volume
  98, pages 063108 (year 2018)NoStop
[Deblais et al.(2018)Deblais, Herrada, Hauner, Velikov, van Roon, Kellay, Eggers, and Bonn]DeblaisPhysRevLett2018
  author author A. Deblais, author M. A. Herrada, author I. Hauner,
  author K. P. Velikov, author T. van Roon, author
  H. Kellay, author J. Eggers,  and author D. Bonn, 10.1103/PhysRevLett.121.254501
  journal journal Phys. Rev. Lett. volume 121, pages 254501 (year
  2018)NoStop


