Past exam of the mathematics course of the University of Cambridge / 2010 / iii / Paper 72 (Solidification of Fluids)

Write $ϵ=1-{\rho }_s/\rho >0$ and $\Delta v=1/{\rho }_s-1/\rho =ϵ/{\rho }_s$. When the ice sphere grows by volume $4\pi a^2\dot{a}dt$, it consumes liquid volume $({\rho }_s/\rho )4\pi a^2\dot{a}dt$. The difference is expelled from the fixed cavity. Thus mass conservation fixes the outward Darcy flux in the surrounding porous medium:The Darcy law $u_r=-(\Pi /\mu ){\partial }_{r}p$ and $p(\infty )=p_m$ implyMatching at the cavity wall, whose internal pressure is uniform, givesThe solid and liquid pressures are equal when the Gibbs--Thomson relation is neglected; nevertheless the coexistence temperature changes with their common pressure because their specific volumes differ. Along coexistence, equality of chemical potentials givesHence the common-pressure melting-point slope is $d{T}_m/dp=-T_m\Delta v/L$. To first order in the pressure change,The $T_m$ multiplying this coefficient is absolute temperature, approximately $273.15\mathrm{K}$ for water; a temperature difference may equally be measured in kelvin or degrees Celsius. Substituting the numerical value $0$ from the Celsius scale into a thermodynamic denominator would be incorrect.

In the quasi-stationary thermal conduction approximation, neglect advective heat transport and use the stated effective conductivity $k$ outside the sphere. The radial Laplace equation gives the bounded solid temperature and the liquid/porous-region temperature:This treats the exterior as a uniform conducting medium, as required by the printed reduction; distinct cavity and porous-medium conductivities would give an additional thermal resistance. There is no solid-side heat flux in the bounded quasi-stationary solution. The Stefan condition therefore becomesEliminate $T_i$ and $p$ to obtainFor $x=a/R$, $t=t_0\tau$ and $t_0=R^2{\rho }_{s}L/(k\Delta T)$, this isThe pressure follows either from the Darcy flux or the Stefan condition:In particular $T_i-T_{\infty }=\Delta TK/(x+K)$: rising pressure progressively reduces the available undercooling. This is spherical freezing inside a porous cavity, whose fixed hydraulic outlet radius $R$ distinguishes it from freezing directly through a porous matrix.

For the idealized zero-radius start, integration givesA finite seed replaces the right side by its difference from the initial value. The reduction ends when the sphere fills the cavity.

For $K\gg 1$, hydraulic drainage is easy and the pressure feedback is small throughout the cavity. The leading expressions areThis is ordinary conduction-controlled spherical solidification, with a completion time asymptotic to $t_0/2$.

For $K\ll 1$, distinguish the initial crossover from the pressure-dominated regime. At $x\ll K$, equivalently $\tau \ll K^2$, the same square-root growth applies and $p-p_m\simeq p^{\ast }x/K$. At $K\ll x\le 1$, equivalently $K^2\ll \tau \lesssim 1/K$,The interface approaches $T_{\infty }$, the pressure approaches the finite coexistence value $p_m+p^{\ast }$, and growth becomes hydraulically limited. The leading completion time is $t_0/(3K)$. The cube-root approximation is not valid at the instant of nucleation, where the initial square-root regime must be retained.

Use the one-sided thermal conduction model encoded by the requested dispersion relation: neglect solid-side perturbation heat flux, take equal phase densities, and let the normal point from solid into liquid. Write the isotropic linear kinetic undercooling law using the solidification attachment mobility $\mathcal{M}>0$:Here $\gamma$ is solid-liquid surface energy, $\mathcal{K}$ is the sum of principal curvatures, and $T_m$ is absolute temperature. Positive curvature convex into liquid lowers the equilibrium freezing temperature through the Gibbs--Thomson relation. The planar base state has $\mathcal{K}=0$.

Set $\Delta T=T_m-T_{\infty }>0$ and defineThe paper's $S$ is the latent-to-sensible heat ratio, the reciprocal of the common Stefan number convention. Below, velocities, lengths and growth rates are dimensionless in these scales. In the travelling coordinate $\xi =X-Vt$, the liquid heat equation gives ${\theta }_0^{\prime \prime }+V{\theta }_0^{\prime }=0$, while the bounded solid steady field is constant. Matching a decaying liquid profile to the Stefan condition $SV=-{\theta }_0^{\prime }(0)$ givesThe kinetic undercooling condition is ${\theta }_i=1-V$. ThusThe speed is constant for $S<1$, determined by $S$ and the attachment coefficient in the specified scales. Thermal diffusivity sets the exponential decay length $\kappa /U$, not the selected speed. Including solid-side perturbation conduction would change the stability equation below.

For linear stability analysis, corrugate the interface to $\xi =\zeta e^{\sigma t+i\alpha y}$ and write the liquid perturbation as $A_1{e}^{-p\xi }{e}^{\sigma t+i\alpha y}$. The liquid heat equation requiresFor this graph the linearized curvature is ${\alpha }^2\zeta e^{\sigma t+i\alpha y}$. Evaluate the temperature on the displaced interface: its perturbation is $A_1+{\theta }_0^{\prime }(0)\zeta =A_1-SV\zeta$. The kinetic undercooling law therefore givesLikewise the derivative evaluated at the displaced surface is ${\theta }_0^{\prime }(0)+{\theta }_0^{\prime \prime }(0)\zeta -p{A}_1$, with ${\theta }_0^{\prime \prime }(0)=S{V}^2$. The Stefan condition at perturbation order readsElimination proves the one-sided kinetic solidification dispersion:The destabilizing term comes from a protrusion encountering colder liquid; capillarity and kinetic relaxation oppose it.

At $S=0$, $V=1$ and the base liquid has no thermal contrast generated by latent heat. The interface obeys $\boxed{\sigma =-\Gamma {\alpha }^2}$: positive surface energy damps all nonzero modes, while $\Gamma =0$ is neutral at this order. Here the thermal perturbation amplitude vanishes, so the auxiliary decay exponent is immaterial to that pure kinetic relaxation. At $S=1$, $V=0$ and the formal limit givesFor $\Gamma >0$ it has a negative real branch $-{\alpha }^2<\sigma <0$ and no growing branch; at $\Gamma =0$ the interface branch is neutral. Strictly, $V=0$ is a singular endpoint of the travelling profile: no finite-width stationary exponential connects ${\theta }_i=1$ to the far-field value zero. These stability statements describe the limit of the travelling-front model.

For $0<S<1$, assume $\sigma =O({\alpha }^2)$ at small wavenumber. ExpandSubstituting into the dispersion relation cancels $S\sigma /V$ on both sides and givesAt the capillary stability threshold for a kinetic freezing front $\Gamma =S/V$, continuing the same expansion gives $\sigma =-2S{\alpha }^4/V^3+O({\alpha }^6)$: the nonzero long waves are damped even at threshold.

For large wavenumber, multiply the dispersion relation by $p$ and eliminate $\sigma$ using the thermal equation. The exact algebraic form isIt prevents the uncontrolled replacement $p\simeq \alpha$ when $\sigma$ itself is of order ${\alpha }^2$. In the usual small-capillarity range the leading results areFor the second result, set $p=\sqrt{1-\Gamma }\alpha +d+O({\alpha }^{-1})$ in the cubic; its quadratic-order balance gives $2d-V=S\Gamma /(1-\Gamma )$. Thus positive capillarity stabilizes sufficiently short waves. If also $\Gamma \ll 1$ and $|\sigma |\ll {\alpha }^2$, the useful reduced approximation isits finite terms should not be treated as a uniform expansion for arbitrary fixed $\Gamma$.

For completeness, large dimensionless capillarity has different balances. At $\Gamma =1$, $p\sim S^{1/3}{\alpha }^{2/3}$ and $\sigma =-{\alpha }^2+S^{2/3}{\alpha }^{4/3}+o({\alpha }^{4/3})$. For $1<\Gamma <1+2S/V$, $p\to p_{\ast }\equiv S/(\Gamma -1)>V/2$, so $\sigma =-{\alpha }^2+p_{\ast }^2-V{p}_{\ast }+O({\alpha }^{-2})$. At $\Gamma =1+2S/V$, the real branch reaches $p=V/2$: for $S<1/3$ its growth rate is $-{\alpha }^2-V^2/4+O({\alpha }^{-4})$, and for $S=1/3$ that endpoint is exact. For $S>1/3$ at this boundary, or for $\Gamma >1+2S/V$, the selected real-$p$ normal mode ceases to exist at sufficiently large wavenumber; the full diffusion problem includes a thermal continuous spectrum and cannot be assigned a real isolated interface growth rate by this ansatz there. Indeed, writing the bulk thermal perturbation as ${\theta }_1=e^{-V\xi /2}u$ reduces its heat equation to $u_t=u_{\xi \xi }-({\alpha }^2+V^2/4)u$, whose half-line diffusion modes have decay rates $-{\alpha }^2-V^2/4-k^2$ for $k\ge 0$. The thermal branch endpoint is $p=V/2$. Loss of a damped isolated interface mode is not a morphological instability.

The approximate instability condition is already clear from the long-wave coefficient: $\Gamma <S/V$, equivalently $S>\Gamma /(1+\Gamma )$. In fact it is exact for growing modes of the present one-sided model. At $\sigma =0$, let $p_0=[V+\sqrt{V^2+4{\alpha }^2}]/2$, so ${\alpha }^2=p_0(p_0-V)$. The right side of the dispersion relation becomesFor real $\sigma \ge 0$, the residual $F=\sigma (1+S/p)-SV+\Gamma {\alpha }^2+S{V}^2/p$ is strictly increasing, sinceTherefore a unique positive root exists precisely when $\Gamma <SV/p_0^2$. For $0<\Gamma <S/V$, the unstable band isThere are no unstable oscillatory roots hidden by this real-root test. Indeed, if $p=a+ib$ and $\mathrm{Re}\sigma \ge 0$, the thermal equation implies $a\ge V$. For $b\ne 0$, the imaginary parts of the two equations would require $2a-V=S(V-\Gamma {\alpha }^2)/|p+S{|}^2\le SV<V$, contradicting $2a-V\ge V$. Thus growing modes have real $p$ and real $\sigma$.

The sketch shows zero capillarity approaching the positive plateau $SV$, weak capillarity giving a positive band followed by short-wave damping, and threshold/stronger capillarity giving damping at every nonzero wavenumber. The zero-wavenumber mode is a translation. As $S\to 1$ with fixed positive $\Gamma$, both the unstable bandwidth and the positive growth rates shrink to zero; the small-wavenumber expansion is not uniform in that limit.

Use a one-sided top-hat plume model with Boussinesq approximation and neglect the solid/liquid density difference in the interfacial volume source. Define the positive ablation speed $U=-\dot{a}$, volume flux per unit wall span $Q=bw$, and kinematic momentum flux $J=b{w}^2$. Let ambient entrainment per height be $E=\epsilon w$, with entrainment coefficient $\epsilon$, and take $\rho ={\rho }_{\infty }[1+{\beta }_C(C-C_{\infty })]$ for the salinity contribution to density. Conservation of volume and vertical momentum, with negligible wall stress and negligible vertical momentum of entering fluid, givesOnly the ocean-facing edge entrains. Using two edges would change the entrainment coefficient in these formulas.

To derive the heat balance, choose liquid specific enthalpy $h_l=c_{p}T+L$ and solid specific enthalpy $h_s=c_p{T}_s$. A horizontal slab receives ambient volume $Edz$ with enthalpy $h_{\infty }=c_p{T}_{\infty }+L$, and solid mass-equivalent volume $Udz$ with enthalpy $h_s$. ThusSubtract $(c_p{T}_{\infty }+L)Q^{\prime }$ and divide by $c_p$:Similarly, the incoming salt flux is $E{C}_{\infty }+U{C}_s$, so $(QC{)}^{\prime }=E{C}_{\infty }+U{C}_s$. Subtracting $C_{\infty }{Q}^{\prime }$ provesBoth anomalies are negative for fresh, cold melt entering the ocean. These balances account for the heat needed to melt the entering ice, rather than treating it as already-liquid water at $T_s$.

Now use $L/c_p\gg T_{\infty }-T_s$, the large latent-to-sensible heat ratio convention of the paper. Put $\Delta C=C_{\infty }-C_s>0$. The stated small plume-deficit limits allow us to replace plume $T,C$ by their ambient values in the diffusive interfacial fluxes. The Stefan condition and interfacial salt balance are thenFor the compositional Péclet number $P\equiv P{e}_c=U\delta /D$, solve the second equation and impose the liquidus:Substitution in the heat balance gives the diffusive glacier-ablation interface lawHere $\Gamma$ is the liquidus slope, unrelated to the surface-energy parameter of the previous question. If $A_T=LD/(c_p\kappa )$, $\Theta =T_{\infty }-T_m$ and $B_T=\Gamma \Delta C$, the physical root isIt is positive when $\Theta +B_T>0$, exactly the condition that the ocean is warmer than its own liquidus. Since $\delta$ and ambient properties are fixed, $U$ is constant with height in this approximation.

For the melting-driven saline wall plume, measure height $Z=z-z_0$ from a virtual origin. The integrated anomaly balances areIn the specified small-deficit regime the meltwater volume is negligible compared with entrainment. Indeed, the solution below gives $(C_{\infty }-C)/\Delta C=4U/(3\epsilon w)$, so a small compositional plume deficit implies $U\ll \epsilon w$. Thus retain $Q^{\prime }=\epsilon w$ for the leading power laws, although $U$ remains essential in the salt and heat sources. Try $b=BZ$, $w=W{Z}^q$. The salt deficit is proportional to $Z^{-q}$, and the momentum balance has powers $Z^{2q}=Z^{1-q}$, hence $q=1/3$. The volume balance gives $(4/3)BW=\epsilon W$, and the momentum coefficient gives $(5/3)B{W}^2=g{\beta }_{C}U\Delta C/W$. ThereforeThese are entrainment-dominated power laws, not an exact solution with the $U$ term retained in the volume equation. They fail close to the virtual origin, where the formal temperature and salinity deficits diverge. Their buoyancy source is distributed melting along the wall, unlike a line plume with a constant buoyancy flux from a localized source.

The profiles rise from $T_i,C_i$ across the laminar boundary layer to the almost-ambient plume values $T,C$, then approach $T_{\infty },C_{\infty }$ across the entraining edge. The sketch uses the linear-gradient approximation implicit in the prescribed diffusive fluxes; it is not an exact normal advection–diffusion profile at large $P$.

For $P\ll 1$,The last formula retains both first-order terms and is appropriate near the no-ablation threshold. Salt diffusion keeps the interface saline, and freezing-point depression allows ice dissolution in saline water even when $T_{\infty }<T_m$. Latent heat is still needed, but the thermodynamic permission to ablate comes from the depressed liquidus; small $P$ does not mean a negligible energy cost.

For $P\gg 1$, $C_i-C_s\sim \Delta C/P$ and $T_i\sim T_m-B_T/P$. Rapid production of fresh melt dilutes the interface, so its temperature approaches the pure-ice melting point. In the usual thermally driven regime $\Theta >0$ and $\Theta \gg B_T/P$,This is predominantly thermal melting. More generally the retained large-$P$ balance is $A_{T}P\simeq \Theta +B_T/P$. For example, $\Theta =0$ gives $P\sim \sqrt{B_T/A_T}$, which can be large while the driving term is still saline freezing-point depression. Thus the melting/dissolution interpretation must also consider the ocean's temperature relative to $T_m$, not the compositional Péclet number alone.

Take $G,V,D,m,C_0>0$ and $T_E<0$, with the liquidus referenced to the pure-water melting point. The solid-liquid segregation coefficient $k_D=0$ means that newly formed ice contains no salt. A steady sharp front producing only pure ice cannot remove the incoming solute flux. More explicitly, in the liquid the steady solute flux is $J_C=-VC-D{C}_z=-V{C}_0$, whereas in pure solid it is zero. Interfacial salt conservation would require $V{C}_0=0$, a contradiction. An equilibrium mushy layer avoids this by carrying concentrated pore liquid down to the eutectic region, where a salt-bearing solid assemblage can form.

Use equal densities and common downward translation velocity $-V$. Write $q=1-\varphi$ for the liquid fraction; $\varphi$ here is the solid fraction. Local liquidus equilibrium in the mush requiresAbove a stationary roof $h_0$, the liquid advection-diffusion equation is $-V{C}_z=D{C}_{zz}$, with solutionThe marginal equilibrium at a mush-liquid boundary condition is $-m{C}_l^{\prime }(h_0)=G$, hence $C_i=C_0+GD/(mV)$. Matching the temperature and liquidus at the roof givesA positive equilibrium mush requires $(-m{C}_0-T_E)/G>D/V$; otherwise this positive-thickness construction does not apply.

Including solute diffusion in the liquid pore space, the steady solute flux in the mush isThus $q[C+(D/V)C_z]=C_0$. Inserting $C_z=-G/m$ and the value of $h_0$ yieldsThis is the diffusive steadily pulled equilibrium mush. Omitting the diffusive contribution would give a different solid-fraction distribution.

The imposed temperature and local liquidus coincide throughout the mush. In the liquid put $y=V(z-h_0)/D\ge 0$. Their difference iswith equality and matching first derivatives at the roof. The inequality follows because the derivative of $y-1+e^{-y}$ is $1-e^{-y}\ge 0$. Thus the liquid is not constitutionally supercooled. In the phase diagram, the pore-liquid trajectory follows the liquidus from the eutectic composition $C_E=-T_E/m$ to $(C_i,T_i)$, then curves into the liquid region and approaches $C=C_0$ as height increases. Pure-ice crystals have composition zero; the trajectory denotes the liquid composition, not the bulk mush composition.

After pulling stops, the imposed temperature keeps the pore-liquid composition fixed wherever the material remains mush. Solute conservation becomesBecause $C_{zz}=0$, this reduces toThe method of characteristics uses $dz/dt=-D/(z+T_E/G)$ and $d\varphi /dt=0$. HenceThese are the solute-diffusion characteristics in a fixed-temperature mush; their validity is local to the coexisting solid–liquid region.

For an explicit initial-data solution, define $H=-T_E/G$, $c=m{C}_0/G$ and $d=D/V$, so $h_0=H-c-d$. A characteristic through $(z,t)$ has initial footWhile that foot lies in the original mush and no phase boundary has removed the point,This recovers the initial distribution and shows that interior ice fraction increases as salt diffuses upward from colder, more concentrated pore liquid.

If one simply transports the initial zero-solid-fraction level, $f(\eta (h_0,0))=0$ implies the printed formal lawOn the physical concentration branch $h<H$, this is $h(t)=H-\sqrt{(H-h_0{)}^2-2Dt}$, so $\dot{h}=D/(H-h)>0$. It starts moving above $h_0$ immediately.

This formal continuation cannot be the physical mush–liquid roof of a column sealed to salt by the membrane. The no-flux constraint in a fixed-temperature mush is $J_C=-Dq{C}_z=0$ at the membrane. At a smooth zero-fraction equilibrium roof, however, $q=1$ and $C_z=-G/m$, giving $J_C=DG/m\ne 0$. Moreover the formal roof has left the membrane-bounded domain. The printed interface law is incompatible with a salt-impermeable membrane when interpreted as the physical roof. Its derivation above is a characteristic extrapolation; it does not enforce the stated boundary condition. If “impermeable” were intended to block bulk flow but allow salt diffusion, a different boundary condition would need to be specified.

The closed-column physical sketch is consequently different. Salt transported upward accumulates below the membrane and produces a liquid brine cap; the cap composition need not lie on the liquidus at every height. Its lower boundary can have a finite jump in solid fraction, so it need not follow the initial zero-fraction characteristic. Near the cold solid base, zero salt supply from below requires the liquid fraction to vanish, and a pure-ice region develops. In the surviving mush, the characteristic formula increases $\varphi$ relative to its initial value. An intermediate profile therefore has pure ice at the bottom, a remaining mush with elevated solid fraction, and liquid at the top. The plotted intermediate roof is schematic, not claimed to obey the inconsistent formal law.

One can characterize the final state quantitatively if the solid base also admits no salt diffusion, as in the usual pure-solid/no-solid-diffusion closure. In a stationary column, solute flux is constant and is zero at the sealed roof. A finite-width equilibrium mush would have $J_C=DGq/m$, forcing $q=0$ everywhere in that region; it therefore cannot remain mush. The final state is pure ice for $0<z<h_f$ and a uniform-composition liquid cap for $h_f<z<h_0$. At their boundary,The cap is liquid because $T(z)>-m{C}_f=T(h_f)$ above its base. The fixed external thermal field can supply or remove the heat involved; this is a solute-budget determination rather than an isolated enthalpy budget.

The final salt budget of a sealed mush column determines $h_f$. Initiallysince $(1-{\varphi }_0)C_m=C_0[1+d/(c+h_0-z)]$. Conservation requires $C_f(h_0-h_f)=N$. Writing the final cap depth ${\ell }_f=h_0-h_f$, the equation and its unique positive root areThe initial salt integral is less than $C_E{h}_0$, so $0<{\ell }_f<h_0$ and $0<h_f<h_0$; the resulting $C_f<C_E$ remains on the ice–brine branch above the eutectic. If the basal eutectic region is instead allowed to exchange solute with the column, this salt budget must include that flux and the numerical final height is not fixed by the printed data alone. No universal physical $h(t)$ follows from the interior characteristic equation without the cap and basal boundary conditions.