Stefan problem

A Stefan problem is a moving-boundary problem in which a phase boundary moves as latent heat is released or absorbed. Temperature satisfies a heat equation in each phase and the Stefan condition supplies the interfacial energy balance.

A two-phase diffusion model with constant properties has similarity solutions built from the complementary error function in each half-line. Equilibrium gives $T_a=-m{C}_a$ and $C_s(a)=k_D{C}_a$. The Stefan condition and interfacial solute balance determine $\lambda$ and $C_a$; finite diffusion on both sides must be retained when constructing the narrow transition near $\lambda =0$.

For $r=-\lambda >0$ and vanishing solid-to-liquid diffusivity ratio, $H(r)=\sqrt{\pi }r{e}^{r^2}\mathrm{erfc}(-r)$ gives the interfacial liquid concentration. It increases from zero to $C_s$ as $r$ increases. Its derivative is $C_s{H}^{\prime }/(1+H{)}^2$, which decreases from $\sqrt{\pi }{C}_s$ to zero. For a direct proof, write $A=H/r=\sqrt{\pi }{e}^{r^2}\mathrm{erfc}(-r)$. Then $A^{\prime }=2rA+2$, so $A\ge \sqrt{\pi }+2r$. Differentiation gives $2H^{\prime 2}-(1+H)H^{\prime \prime }=2[(1+r^2+2r^4)A^2+r(4r^2-1)A+2r^2-2]>0$: the bracket increases with $A$, and substitution of its lower bound gives only positive polynomial coefficients, including constant $\pi -2$. Thus the concentration is strictly concave as a function of melting speed parameter $r$. Combining this with a linearized thermal balance produces one fold on the melting branch and, with the freezing branch, an interval of three possible interface speeds.

For a one-phase spherical Stefan problem with melting-temperature boundary value $T_s-\Gamma /a$ and supercooling $\Delta T=T_s-T_{\infty }>0$, use the latent-to-sensible heat ratio $S=\mathcal{L}/(c_p\Delta T)$. Large $S$ separates heat diffusion from slow interface motion. The exterior harmonic temperature is $T=T_{\infty }+(\Delta T-\Gamma /a)a/r$. The Stefan condition gives the displayed radius equation with $a_c=\Gamma /\Delta T$. This is the large inverse-Stefan number convention; the common sensible-to-latent Stefan number is $1/S$.

The Gibbs--Thomson relation depresses a spherical interface's equilibrium temperature by $\Gamma /a$. For $a>a_c$, its interface is warmer than the supercooled ambient and heat loss drives growth. For $a<a_c$, its interface is colder than the ambient and heat gain drives melting. The equilibrium radius is unstable because the derivative of the quasistatic radius velocity at $a_c$ is positive.

In the quasistatic spherical model, ${\eta }_0=a_0/a_c\in (0,1)$. Separating the radius equation gives $t_e=(S/\kappa ){\int }_0^{a_0}{a}^2/(a_c-a)da$, which is finite and evaluates to the displayed expression. Near disappearance, $a^3\sim 3\kappa a_c(t_e-t)/S$. The time diverges logarithmically as the initial radius approaches the critical radius from below.

A spherical solid of radius $a$ growing inside a fixed cavity of radius $R$ expels volume at rate $4\pi (1-{\rho }_s/{\rho }_l)a^2\dot{a}$. The Darcy flux outside the cavity gives $p-p_m=\mu (1-{\rho }_s/{\rho }_l)a^2\dot{a}/(\Pi R)$. The common-pressure melting-point slope and quasi-steady spherical thermal conduction give ${\rho }_{s}La\dot{a}=k(T_i-T_{\infty })$. Define $x=a/R$, $\tau =tk(T_m-T_{\infty })/(R^2{\rho }_{s}L)$ and $K=L^2\Pi /[k{T}_m\mu (1/{\rho }_s-1/{\rho }_l{)}^2]$. Elimination yields $(x+x^2/K)x_{\tau }=1$ and the displayed zero-seed integral. The pressure rise is $p^{\ast }x/(x+K)$, with $p^{\ast }=L(T_m-T_{\infty })/[T_m(1/{\rho }_s-1/{\rho }_l)]$. Easy drainage gives square-root growth; low permeability gives an initial square-root regime followed by cube-root growth and pressure saturation. The model stops at $x=1$ and uses a uniform exterior conductivity.

A crystal dendrite is a branched crystal-growth morphology. Diffusive access to undercooling or supersaturation favours protruding tips, while interfacial capillarity regularizes fine scales and crystal interface anisotropy selects directions. Dendritic growth of a snowflake is a vapour-growth example.

For equal thermal properties and bounded travelling profiles on two deep half-spaces, the solid temperature is $T_m$ and the liquid profile is exponential. The Stefan condition requires $L=c_p\Delta T$. This balances the latent heat of solidification with warming of the initially supercooled liquid.

Retaining perturbation thermal diffusion givesThe displaced interfacial temperature and Stefan condition give the displayed implicit dispersion relation. In the capillarity-free limit its growing branch is exactly $\sigma =V\alpha$, which checks the limits of the small-wavelength thermal approximation for a solidification front.

Approximating perturbation temperatures by harmonic fields decaying as $e^{\pm \alpha \xi }$, while retaining displacement through the base gradient, produces this reduced dispersion relation. The Gibbs--Thomson relation supplies $\Gamma =\gamma T_m/(\rho L)$. This short-wave reduction suppresses perturbation advection and time dependence; its constant term and a threshold at $\alpha \kappa /V=O(1)$ are not controlled next-order results of the full diffusion problem.

The reduced dimensionless growth rate $g(q)=-1+q-2q^3/R_c$ has a positive maximum precisely when $R_c>27/2$. Its maximizing $q=\sqrt{R_c/6}$ follows by differentiation. At equality $q=3/2$, outside a strict $q\gg 1$ short-wave regime. The inequality is an exact criterion for the reduced cubic, not a quantitative onset theorem for full thermal stability of a travelling solidification front.

At zero liquid superheat, dense solid released by viscous overturn of a dense crust can sink and remain solid. For a thin crust on an asteroid of fixed radius $R$, efficient deposition and the equal-density volume approximation give $b^3-b_0^3=3R^2{J}_s(t-t_0)$. If surface gravity scales as $R$, the crustal flux scales as $R^{1/3}$ and the core radius as $R^{7/9}{t}^{1/3}$. The approximation ends when remaining liquid is exhausted or the renewed crust is no longer thin.

A growing ice-filled sphere expels liquid because ice is less dense than water. With porosity $\varphi$, $ϵ=1-{\rho }_s/\rho$, and pressure-induced freezing-point depression $T_i=T_m-\gamma p_i$, the Darcy velocity gives $p_i=\mu \varphi ϵa\dot{a}/\Pi$. Quasi-steady spherical thermal conduction gives $\varphi {\rho }_{s}La\dot{a}=k(T_i-T_{\infty })$. Hence $a\dot{a}=k(T_m-T_{\infty })/[\varphi {\rho }_{s}L(1+\Lambda )]$, where $\Lambda =k\gamma \mu ϵ/(\Pi {\rho }_{s}L)$. The interface approaches the far-field temperature at low permeability of a porous medium, while its pressure tends to a finite equilibrium value.

The Stefan condition equates latent heat generated by motion of a phase boundary to the jump in conductive heat flux across it. With interface speed $V$ and normal directed from solid to liquid, one common sign convention is

Kinetic undercooling makes the speed of a phase boundary proportional to the difference between its equilibrium and actual interface temperatures. A curvature correction raises the required undercooling for protrusions and therefore stabilizes short wavelengths.

The attachment mobility $\mathcal{M}>0$ is the coefficient converting interfacial undercooling into a solidification speed. In the linear isotropic law, $U=\mathcal{M}(T_{\mathrm{eq}}-T_i)$, and capillarity makes $T_{\mathrm{eq}}=T_m-{\gamma }_T\mathcal{K}$. Thus $T_i=T_m-U/\mathcal{M}-{\gamma }_T\mathcal{K}$. This kinetic coefficient refers to a moving phase boundary, rather than an order-parameter relaxation coefficient. The Gibbs--Thomson relation determines ${\gamma }_T=T_m\gamma /(\rho L)$ for isotropic solid-liquid surface energy $\gamma$.

Instantaneous interface kinetics is the limit of negligible kinetic undercooling at a moving solidification interface. Its temperature is set by local liquidus equilibrium and curvature-induced melting-temperature depression, without a further velocity-dependent attachment term. It does not imply that diffusion or heat transport is instantaneous.

A planar front growing into a supercooled melt can amplify corrugations because protrusions enter colder liquid. Interface curvature opposes this effect. Linearization couples the perturbation growth rate to the decaying thermal exponent ahead of the front.

Use undercooling $\Delta T=T_m-T_{\infty }$, solidification attachment mobility $\mathcal{M}$, velocity scale $U_0=\mathcal{M}\Delta T$, length $\ell =\kappa /U_0$ and time $\kappa /U_0^2$. Let $S=L/(c_p\Delta T)$ be the latent-to-sensible heat ratio, $\theta =(T-T_{\infty })/\Delta T$, and $\Gamma ={\gamma }_T/(\Delta T\ell )$. In the one-sided liquid-diffusion model, a travelling planar front has $\theta =S{e}^{-V\xi }$ and $V=1-S$. For a corrugation $\zeta e^{\sigma t+i\alpha y}$, the decaying thermal perturbation is $A_1{e}^{-p\xi }$, with $p^2-Vp=\sigma +{\alpha }^2$ and the branch $p=[V+\sqrt{V^2+4(\sigma +{\alpha }^2)}]/2$. The displaced kinetic condition gives $A_1=(SV-\sigma -\Gamma {\alpha }^2)\zeta$. The displaced Stefan condition gives $S\sigma \zeta =p{A}_1-S{V}^2\zeta$. Eliminating $A_1$ proves the displayed dispersion relation. Neglect of solid-side perturbation heat flux is part of this model.

The cubic follows by eliminating $\sigma =p^2-Vp-{\alpha }^2$ from one-sided kinetic solidification dispersion, with $V=1-S$. For $\Gamma =0$, the interface branch tends to $\sigma =SV-SV/\alpha +O({\alpha }^{-2})$. For fixed $0<\Gamma <1$, substituting $p=\sqrt{1-\Gamma }\alpha +d$ gives $2d-V=S\Gamma /(1-\Gamma )$, hence $\sigma =-\Gamma {\alpha }^2+[S\Gamma /\sqrt{1-\Gamma }]\alpha +O(1)$. At $\Gamma =1$ the balance is $p^3\sim S{\alpha }^2$, giving $\sigma =-{\alpha }^2+S^{2/3}{\alpha }^{4/3}+o({\alpha }^{4/3})$. For $1<\Gamma <1+2S/V$, $p\to S/(\Gamma -1)>V/2$ and $\sigma =-{\alpha }^2+p^2-Vp+o(1)$. Stronger capillarity can remove the real isolated interface mode from the chosen thermal branch; it does not create growth. These distinctions explain why replacing $p$ by $\alpha$ is not a controlled fixed-capillarity limit.

For $0<S<1$ and $V=1-S$, the translational branch of one-sided kinetic solidification dispersion has $\sigma =(S/V-\Gamma ){\alpha }^2+O({\alpha }^4)$. Expand $p=V+(\sigma +{\alpha }^2)/V+O({\alpha }^4)$ and substitute to prove this limit. At zero growth, $p_0=[V+\sqrt{V^2+4{\alpha }^2}]/2$, and the destabilizing right side has sign ${\alpha }^2(SV/p_0^2-\Gamma )$. For nonnegative real $\sigma$, the residual has derivative at least one, so there is one growing root precisely when $\Gamma <SV/p_0^2$. Thus instability exists exactly when $\Gamma <S/V$, with unstable band ${\alpha }^2<p_c(p_c-V)$, $p_c=\sqrt{SV/\Gamma }$ for $\Gamma >0$. At threshold the next term is $\sigma =-2S{\alpha }^4/V^3+O({\alpha }^6)$. Surface energy removes short-wave growth even when long waves remain unstable.

An icicle ripple is a corrugation of the ice-water interface produced while a thin flowing water film freezes. Coupling lubrication theory, surface tension, thermal conduction, and the Stefan condition determines its growth and migration.

For a water film of thickness $h$, lateral freezing speed $V$, ripple wavenumber $\alpha$, and $\Gamma =\gamma /(3\rho g)$, a long-wave model gives the dispersion relationThe fastest-growing mode has ${\alpha }_m=5^{-1/6}(\Gamma h{)}^{-1/3}$ and travels upward at speed $\sqrt{5}V/6$ when the vertical coordinate points downward.

A phase boundary is an interface separating distinct phases of matter. Its temperature, velocity, and flux jumps are constrained by phase equilibrium and conservation laws.

A planar freezing front can push a spherical particle steadily when repulsive sphere-plane dispersion force balances lubrication resistance. The thin parabolic gap gives hydraulic drag $6\pi \mu V{a}^2/d$, and balancing it against $Aa/(6d^2)$ gives the displayed speed. This balance assumes a nonretarded interaction, small gap and negligible inertial and other particle forces.

For a sphere partially below a planar melting isotherm, write $c=\cos {\theta }_0$ and depth $a(\cos \theta -c)$. Local premelting gives film thickness $d^3=K/(\cos \theta -c)$. Conservation in the thin spherical film gives flux $Va\sin \theta /2$ and pressure $p=3\mu V{a}^2(\cos \theta -c{)}^2/K$. Integrating vertical repulsive and hydraulic stresses over the wetted cap yields the displayed speed. The unknown thermodynamic constant $K$ cancels, but a narrow rim where the film ceases to be thin is excluded from the leading approximation.

The solid-liquid segregation coefficient is the ratio of the solute concentration in newly formed solid to that in the adjacent equilibrium liquid. A coefficient of zero means complete solute rejection: newly formed crystals are pure solvent. It does not mean the bulk solute in a mushy layer vanishes, because its pore liquid can remain concentrated.

A crystal's interfacial free energy and attachment kinetics depend on the interface normal $\mathbf{n}$. The orientation-dependent energy $\gamma$ and mobility $M$ select directions even when bulk diffusion is isotropic. This supplements the transport mechanism of crystal dendrites; isotropic transport alone does not specify the number of preferred arms.

Ordinary atmospheric snow crystals are predominantly hexagonal ice Ih. Rotating by sixty degrees around its crystallographic axis leaves the equilibrium orientation dependence unchanged, giving six equivalent basal-plane directions. Temperature and supersaturation determine whether plates, columns or crystal dendrites develop. Dendritic growth of a snowflake needs this crystal symmetry as well as the diffusion instability.

Premelting is the persistence of a liquid layer at a surface or interface below the bulk melting temperature. Molecular disjoining pressure and phase equilibrium can stabilize the film. Its hydraulic connection to a reservoir permits frost heave when water freezes onto the base of loaded ice.

Thermomolecular pressure is the temperature-dependent force per area transmitted through a thin liquid film between a solid and a substrate. For a repulsive ice-film interaction $p_T=A/(6\pi a^3)$, it can drive water through a porous material and support a load. Its thermodynamic relation to undercooling is the Clapeyron pressure relation for equal-density phases.

For equal solid and liquid mass densities, $d\mu =-sdT+dp/\rho$ and $s_l-s_s=L/T_m$ imply the displayed relation from equality of chemical potentials. $T_m$ must be absolute temperature. The phase pressures differ; this should not be confused with differentiating a coexistence curve along a common pressure.

Frost heave is displacement caused by segregated ice growth supplied by liquid flow, rather than solely by the expansion of water on freezing. A premelting film allows liquid access beneath the ice; thermomolecular pressure draws supply against viscous resistance and load.

In an equal-density, no-slip lubrication approximation, repulsive disjoining pressure $A/(6\pi d^3)$ sustains a premelting film of thickness $d$ beneath a rigid disk of radius $a$. With imposed radial undercooling $\Delta T(1-r^2/a^2)$, local thermomolecular pressure is $P_T=\rho L\Delta T(1-r^2/a^2)/T_m$, and $d^3=A/(6\pi P_T)$. A constant vertical lifting rate $U$ requires inward volume flux $q_r=-Ur/2$. Integrating the resulting lubrication pressure and disk force balance gives the displayed rate, where $W$ is the apparent load, $\mu$ is liquid dynamic viscosity, $L$ is specific latent heat, and $T_m$ is the absolute bulk melting temperature. Positive lifting requires $W<\pi a^2\rho L\Delta T/(2T_m)$. For fixed positive load, the rate rises from zero at that threshold towards $A/(6\pi \mu a^2)$ as undercooling increases. That plateau is an asymptotic limit within this constitutive model, not a claim that continuum-film physics persists at arbitrarily low temperatures.

A porous sheet of thickness $d$ and permeability of a porous medium $\Pi$ supplies a thin premelting film. The Darcy flow law equates water supply to ice growth for equal densities. With equal film and ice conductivities, negligible latent heat relative to the through-going conductive flux, and $a\ll h$, phase equilibrium gives $p_T=(\rho L\Delta T/T_m{)}^{3/4}(A/(6\pi h^3){)}^{1/4}$. The increasing gravitational load limits growth. A quasi-steady temperature field alone does not imply the negligible-latent-heat approximation.

For ${\eta }^{\prime }={\eta }^{-3/4}-\eta$, set $u={\eta }^{7/4}$ to get $u^{\prime }=(7/4)(1-u)$. The displayed solution approaches the stable equilibrium $\eta =1$. From the formal zero start it behaves as $(7\tau /4{)}^{4/7}$ at early times and $1-(4/7)e^{-7\tau /4}$ at late times. The zero-start asymptote is an intermediate reduced-model regime, since the physical thin-film approximation fails at vanishing ice thickness.

Combine the full quasi-steady Stefan condition above with $\delta T=T_m{p}_T/(\rho L)$, $a=(A/(6\pi p_T){)}^{1/3}$ and the Darcy flow law to obtain an implicit pressure-growth system. Even when latent heat is negligible, a conductivity ratio $K_s/K_l=r$ changes the leading molecular driving pressure by $r^{3/4}$. Thus a reduced equal-conductivity formula requires an explicit material and rate approximation.

The excess free energy per area of a solid-liquid phase boundary. Curvature changes the equilibrium temperature through the Gibbs--Thomson relation. Crystal anisotropy can make this energy and its associated surface stiffness depend on orientation.

The Neumann solution places a planar phase boundary at $h(t)=2\lambda \sqrt{\kappa t}$. Similarity reduction of the heat equation gives temperature profiles in terms of the error function, while the Stefan condition determines $\lambda$ implicitly.

For equal phase thermal properties, a roof undercooling $\Delta T$ and liquid superheat $\Theta$, the planar Stefan problem has interface $a=2\lambda \sqrt{\kappa t}$. The solid error function profile and liquid complementary-error-function profile give $S\lambda =e^{-{\lambda }^2}[1/\mathrm{erf}\lambda -(\Theta /\Delta T)/\mathrm{erfc}\lambda ]/\sqrt{\pi }$. Here $S=L/(c_p\Delta T)$ is the latent-to-sensible heat ratio. Liquid heat supplied to the interface reduces crust growth.

The small-$\lambda$ limit of superheated two-phase Stefan similarity gives $a\simeq A\sqrt{\kappa t}$, where $A=2/[\sqrt{{\vartheta }^2/\pi +2S}+\vartheta /\sqrt{\pi }]$ and $\vartheta =\Theta /\Delta T$. For fixed superheat at large latent-to-sensible heat ratio, $A\sim \sqrt{2/S}$. For much larger superheat, $A\sim \sqrt{\pi }/\vartheta$. The approximation retains the leading effect of a hot interior and its coefficient need not be exact when $S$ is only order one.

Constitutional supercooling occurs when liquid ahead of a solidification phase boundary lies below its local liquidus. In a saline liquid with $T_L=-mC$ and $m>0$, salt rejection enriches the interface liquid and depresses its liquidus; as the concentration decreases away from the interface, the liquidus rises. If it rises faster than the actual temperature, a supercooled interval appears. This can make a planar solidification phase boundary unstable.

For a saline Stefan problem with common thermal properties, $ϵ=\sqrt{D/\kappa }\ll 1$, and fixed latent-to-sensible heat ratio, set $T_0=-m{C}_0$, ${\theta }_{-}=T_0-T_{-\infty }$ and ${\theta }_{+}=T_{\infty }-T_0>0$. A finite-$\lambda$ similarity solution has $T_i=(T_{\infty }+T_{-\infty })/2+O(ϵ)$. At the onset of constitutional supercooling, $\lambda =O(ϵ)$ and the error in this mean is $O({ϵ}^2)$. The salt-rejection function for a saline Stefan front then gives $m\lambda C_i=({\theta }_{-}-{\theta }_{+})/(2\sqrt{\pi })+O({ϵ}^2)$, whereas $\sqrt{Dt}{T}_{l,x}(a^{+})=ϵ{\theta }_{+}/\sqrt{\pi }+O({ϵ}^2)$. Comparing the actual temperature gradient to the liquidus gradient gives the displayed onset, with supercooling on the side of greater ${\theta }_{-}$. For a freezing front at $a\ge 0$ and $D<\kappa$, the ratio of these two gradients decreases into the liquid, so their interface comparison also decides whether a supercooled interval occurs anywhere ahead. The formula is an asymptotic criterion; finite-$ϵ$ onset requires the full Stefan condition and salt balance.

A saline Stefan problem couples a heat equation, a solute diffusion equation, a Stefan condition, solute conservation, and a concentration-dependent liquidus at a moving solid-liquid boundary. Different thermal and solutal diffusivities produce nested boundary layers.

In a planar no-flow model with pure ice, a stable cold brine below freshwater can support an ice layer whose upper phase boundary freezes freshwater while its lower boundary melts into diluted brine. Heat conduction and salt diffusion obey different similarity scales. Local liquidus equilibrium, salt conservation and two Stefan conditions determine both boundary positions. The heat lost into the brine supports net ice growth even while some basal ice melts.

For an ice layer between freshwater and cold brine, set $ϵ=\sqrt{D/\kappa }$ and use solutal similarity positions $a=2{\lambda }_a\sqrt{Dt}$, $b=2{\lambda }_b\sqrt{Dt}$. A linear ice temperature, two Stefan conditions and interfacial salt conservation determine both negative interface constants. As $ϵ\to 0$ with fixed latent-to-sensible heat ratio $S$, ${\lambda }_b=O(1)$ but ${\lambda }_a=O({ϵ}^{-1})$. Freezing above occurs over the thermal diffusion scale, while dilution-driven melting below has the much shorter solutal scale.

For $S=L/(c_p\Delta T)\gg 1$, $ϵ=\sqrt{D/\kappa }\ll S^{-1}$ and $R=m{C}_0/\Delta T>1$, the leading thickness is the displayed expression. The basal position is $2A\sqrt{Dt}$, with $A$ determined byOnly when $ϵSA\ll 1$ can the basal translation be omitted from the upper position. Then $A{e}^{A^2}\simeq \sqrt{\pi }RS/4$, exposing the logarithmic refinement to the scale ordering.

For stationary liquid below a salt-free solid and an upward melting interface $a(t)=2A\sqrt{Dt}$, the complementary error function salt profile and flux balance $D{C}_z(a)=-\dot{a}{C}_i$ give the displayed concentration ratio. For $A>0$ it lies below one. At large $A$ it behaves as $e^{-A^2}/(2\sqrt{\pi }A)$.

When salt-free ice forms from a saline liquid, the moving phase boundary excludes the dissolved salt. With equal phase mass densities, no bulk flow, liquid salt diffusion coefficient $D$, and interface concentration $C_i$, salt conservation gives the displayed condition, where $x$ increases from solid to liquid. Positive freezing speed $\dot{a}$ therefore requires a negative liquid concentration gradient. During melting the same condition gives dilution. Coupling salt rejection, the liquidus, and the Stefan condition produces the saline Stefan problem.

For salt-free ice and a liquid concentration profile proportional to $\mathrm{erfc}[x/(2\sqrt{Dt})]$, salt rejection at $a=2\lambda \sqrt{Dt}$ gives $C_{i}F(\lambda )=C_0$, where $D$ is the salt diffusion coefficient, and $C_i,C_0$ are the interface and far-field concentrations. This follows from $-D{C}_x(a^{+})=\dot{a}{C}_i$. The function is positive and strictly decreasing for all real $\lambda$, since an integration by parts givesConsequently $C_i=C_0/F(\lambda )$ increases with $\lambda$, equals $C_0$ at zero, and is enriched for freezing and diluted for melting. The positivity proof also prevents spurious negative interface concentrations when using the complementary error function formula.

For a planar ice-ocean interface with no fluid motion, similarity places the boundary at $a(t)=2\lambda \sqrt{Dt}$. The sign of $\lambda$ distinguishes ice growth from ablation, while rejected salt during growth can cause constitutional supercooling in the adjacent ocean.

A mushy layer contains solid crystals and liquid in local phase equilibrium, rather than a single sharp phase boundary. Its solid fraction changes as heat and solute are transported. A concentration-dependent liquidus couples temperature to liquid composition; a connected crystal framework makes the residual liquid behave as a porous-media flow.

A uniform-permeability mushy layer on an impermeable eutectic floor, under the quadratic roof pressure $p_0-E^2{x}^2/2$, has Darcy flux $\mathbf{q}=\alpha (x,-z)$, where $\alpha =\Pi E^2/\mu$. With pulling velocity $-V{\mathbf{e}}_z$, pure primary crystals, solid fraction $\varphi$, liquid concentration $C$ and constant phase properties, its one-dimensional balances areThe solidification product at the eutectic floor retains bulk composition $C_s=(1-\varphi (0))C_E$. Integration gives $C_s=C_0-(\alpha /V){\int }_0^h(C-C_0)dz$: horizontal strain removes solute before the material freezes. A negative predicted $C_s$ is inadmissible, and diagnoses failure of the assumed steady branch rather than a negative-composition product.

Set $\delta =\sqrt{2\kappa /\alpha }$, $R=\alpha h/V$ and $S=L/[c_p(T_0-T_E)]$. In the regular limit $R\gg S\gg 1$, at fixed finite thermal strain number $\alpha h^2/\kappa$ and bounded phase fraction, the heat equation reduces to $\kappa T^{\prime \prime }+\alpha z{T}^{\prime }=0$. Integration gives the displayed error function profile. Matching the roof temperature gradient $G$ yieldsThe product composition is $C_s=C_0-\alpha (T_0-T_E)\delta [1-e^{-(h/\delta {)}^2}]/[mV\sqrt{\pi }\mathrm{erf}(h/\delta )]$. The asymptotic reduction requires admissible phase fractions. If the thermal strain number grows singularly with $R$, the two ordering inequalities alone do not control boundary-layer errors; the full strain-driven steady mush equations retain the latent-heat term.

In a stationary mushy layer with fixed $T=T_E+Gz$, pure solid and local liquidus $C=-T/m$, solute conservation is ${\partial }_t[(1-\varphi )C]=D{\partial }_z[(1-\varphi )C_z]$. Since $C_z=-G/m$ is constant, it reduces to $(z+T_E/G){\varphi }_t=D{\varphi }_z$. The method of characteristics makes the displayed $\eta$ constant and gives $\varphi =f(\eta )$ within the region that remains mush. For negative $T_E$, the characteristics move upward where the equilibrium concentration is positive. Extrapolating a zero-solid-fraction level is not sufficient to determine a boundary when additional flux constraints create a liquid cap or a solid-fraction jump.

For nonzero imposed thermal gradient and linear liquidus equilibrium, $C_z=-G/m\ne 0$ wherever a mixture of crystals and pore liquid persists. A sealed boundary requires zero solute flux, so a local-equilibrium mush trace there would require $1-\varphi =0$, not a zero-solid-fraction roof. In particular a pure-liquid trace with the equilibrium concentration gradient has nonzero diffusive flux. A sealed column must develop a nonequilibrium-in-the-liquid temperature/composition relation, a separate liquid cap, or an appropriate phase-fraction jump; simply transporting the initial zero-fraction level violates the boundary condition.

Assume a column is sealed to salt at both ends, diffusion occurs only in liquid, and its imposed temperature increases strictly upward. A finite-width equilibrium mushy layer cannot be stationary: its solute flux is nonzero unless its liquid fraction vanishes. The final state is pure ice below a height $h_f$ and uniform-composition brine above it. Their interface obeys $C_f=-(T_E+G{h}_f)/m$. The displayed mass conservation equation fixes $h_f$ from the initial salt content. For a diffusive steadily pulled equilibrium mush, let $H=-T_E/G$, $c=m{C}_0/G$, $d=D/V$ and $h_0=H-c-d$. The initial salt content is $N=C_0[h_0+d\log (1+h_0/c)]$. Thus the liquid-cap depth ${\ell }_f=h_0-h_f$ obeys ${\ell }_f({\ell }_f+c+d)=c[h_0+d\log (1+h_0/c)]$, which has one positive root. This conclusion assumes no subsequent exchange with an underlying eutectic reservoir.

A stagnant mushy-layer model has no bulk advection, equal phase mass densities, pure solid, and local equilibrium on the liquidus. With liquid fraction $\varphi$, specific enthalpy $E=c_{p}T+L\varphi$ and bulk solute concentration $B=\varphi C$, neglecting solute diffusion inside the mushy layer gives $T_t+(L/c_p){\varphi }_t=\kappa T_{zz}$ and $(\varphi C{)}_t=0$. The latter fixes composition at each stationary material level; it does not automatically set every level to the initial liquid concentration when solute can diffuse ahead of the growing mush.

For $\mathcal{C}=C_0/\Delta C\gg 1$ and inverse Stefan number $\mathcal{S}=O(\mathcal{C})$, a stagnant mushy-layer model has nearly constant effective heat capacity, giving ${\Omega }^2{\theta }_t=\kappa {\theta }_{zz}$ with ${\Omega }^2=1+\mathcal{S}/\mathcal{C}$. With negligible liquid solutal diffusivity, roof value $\theta =-1$, positive liquid superheat ${\theta }_{\infty }$, and edge $h=2\lambda \sqrt{\kappa t}$, the mush profile is $\theta =-1+\mathrm{erf}(\Omega \eta )/\mathrm{erf}(\Omega \lambda )$. Matching heat flux requires $\Omega e^{{\lambda }^2}\mathrm{erfc}\lambda ={\theta }_{\infty }{e}^{{\Omega }^2{\lambda }^2}\mathrm{erf}(\Omega \lambda )$. The square on $\Omega$ in the exponential is essential. The porosity $\mathcal{C}/(\mathcal{C}-\theta )$ records the leading composition balance; the linear thermal profile is asymptotic, not exact at finite $\mathcal{C}$.

At a sharp mush–liquid edge with linear liquidus $T_L=T_m-mC$, marginal equilibrium makes the liquid temperature gradient tangent to the local liquidus: $T_z=-m{C}_z$. Together with $T=T_L$ at the edge and $T\ge T_L$ in the exterior liquid, it supplies the usual absence-of-constitutional supercooling closure. When solute diffusion is neglected inside the mushy layer but retained outside it, interfacial solute conservation can require a small solid fraction jump; imposing unit porosity as an additional finite-diffusivity condition generally overdetermines this sharp-edge reduction.

The solid fraction is the volume fraction occupied by solid in a solid–liquid mixture. The complementary liquid fraction is $1-\phi$.

The liquid fraction is the volume fraction occupied by liquid in a solid–liquid mixture. For a saturated rigid mushy layer, it is the porosity available to Darcy flux.

For constant phase mass densities ${\rho }_s,{\rho }_l$, a fixed solid framework has mixture density $\bar{\rho }={\rho }_s\phi +{\rho }_l(1-\phi )$. Mass conservation giveswhere $\mathbf{u}$ is Darcy flux. Equal phase densities give a divergence-free flux. Unequal densities require expansion or contraction flow during freezing or melting.

Directional solidification advances a freezing region preferentially along an imposed temperature gradient. Pulling material through fixed heat exchangers can make the mushy layer stationary in the apparatus frame even while solid and liquid pass through it.

The frozen-temperature approximation prescribes the external thermal field during interface perturbations, rather than solving a perturbation heat equation. An interface displacement samples that prescribed field at its new location. Local liquidus and Gibbs-Thomson effect conditions still apply, and the solute field may retain full time-dependent diffusion.

For pure solid rejecting all solute, pulling speed $V>0$, solutal diffusivity $D$ and finite solute inventory per unit area $N$, a stationary planar front has $C(z)=C_I{e}^{-Vz/D}$ in its liquid half-space. The advection-diffusion equation gives this profile, and ${\int }_0^{\infty }Cdz=N$ fixes $C_I=VN/D$. A liquidus $T_I=-m{C}_I$ fixes its position in an imposed temperature gradient. The pure-solid branch requires $C_I<C_E$, where $C_E$ is the eutectic composition.

For an interface perturbation $\eta e^{ikx+st}$ with instantaneous interface kinetics, take the liquid concentration perturbation proportional to $e^{-qz}$. Here $\Gamma >0$ is the coefficient of curvature-induced melting-temperature depression. The advection-diffusion equation and interfacial mass conservation giveFor a growing normal mode the unique decaying liquid branch has $q=(V+\sqrt{V^2+4D(s+D{k}^2)})/(2D)$. Writing $b=(G+\Gamma k^2)/(m{C}_I)$, its rate is $s=(D/2)(V/D-b)(\sqrt{b^2+4k^2}-b)$. Instability requires $G+\Gamma k^2<mV{C}_I/D$. The Gibbs-Thomson effect suppresses sufficiently short waves; a quasistatic solute diffusion to solute diffusion changes the rate and must be stated separately.