Past exam of the mathematics course of the University of Cambridge / 2026 / iii / Paper 332 (Fluid Dynamics of the Solid Earth)

The upward arrows at the surface are reflected short-wave flux $\alpha F_r$, emitted long-wave flux $ϵ\sigma T_s^4$ from the Stefan–Boltzmann law, and, when positive, conductive heat arriving from the planetary ice shell. The downward arrow is the incident solar flux $F_r$. The planetary surface energy balance is thereforeThe relevant planetary albedo is the Bond albedo. In a steady one-dimensional shell, Fourier's law and $k(T)=a/T$ make the upward conductive fluxwhere $z$ points downward and $T(h)=T_m$. Substitution gives the requested implicit equation

Conservation of energy and Fourier's law giveAt steady state $k{T}_z=F_c$ is constant. Since $k=a/T$,The two temperature boundary conditions therefore giveAt the basal phase boundary, the Stefan condition isThe first term removes latent heat released by freezing, while the oceanic flux $F_w$ supplies heat to the interface. A steady shell has $\dot{h}=0$, hence $F_c=F_w$ andHere $T_s$ must simultaneously satisfy the surface balance from part a with $F_c=F_w$.

WriteThe quasi-steady model is the coupled pairFor a thin shell, $T_s$ is close to $T_m$. A first-order asymptotic expansion givesand thereforeTo leading order $F_c\sim F_{\ast }$, so an initially freezing shell grows linearly:

For a thick shell the conductive correction in the surface balance is small. The radiative equilibrium temperature and its first correction areThus, with $A=a\log (T_m/T_e)$,When $F_w$ is negligible over an intermediate range, this ordinary differential equation gives the square-root growth lawRetaining $F_w$, its implicit solution isConsequently a sketch of $h(t)$ starts approximately linearly, crosses to square-root growth, and approaches $h_{\infty }$ with exponential decay of $h_{\infty }-h$. A shell placed above $h_{\infty }$ instead thins because the basal oceanic heat flux exceeds the conductive loss.

Treat the warm ice shell as a very viscous layer and estimate its Rayleigh numberThe temperature-dependent ice viscosity should be evaluated carefully, because deformation is concentrated near the warm base; a useful first estimate uses a representative basal or depth-averaged viscosity. The onset time is obtained by inserting $h(t)$ and $\Delta T(t)=T_m-T_s(t)$ into this expression and finding when $\mathrm{Ra}$ first exceeds the critical value for the shell's mechanical boundary conditions, usually of order $10^3$.

Above onset, mantle convection within the ice transports heat more efficiently than thermal conduction, increases the basal heat loss for a given thickness, and generally limits further thickening. For otherwise similar bodies, a larger planetary radius normally gives larger gravitational acceleration and hence a larger Rayleigh number. Convection is therefore more likely and begins in a thinner or younger shell, although changes in pressure-dependent melting temperature, shell thickness, and viscosity can alter that comparison.

Take $z$ vertically downward from the cold surface. The ice occupies $0<z<h(t)$ and the water occupies $z>h(t)$. The temperature satisfies the heat equationin both phases, withThe remaining interfacial condition is the Stefan conditionA schematic should show the cold surface, the descending phase boundary, the ice temperature rising from $T_b$ to $T_m$, and the water temperature tending from $T_m$ to $T_{\infty }$.

The diffusion length suggests the similarity variableThe Neumann solution of the Stefan problem, written using the error function, isThese expressions satisfy all four thermal boundary conditions. At the interface,Since $\dot{h}=\lambda \sqrt{\kappa /t}$ and $k=\rho c_p\kappa$, the Stefan condition becomesThe second term represents heat supplied from the warm water and therefore reduces the freezing rate.

At times long compared with the encounter time at depth $H$, the fixed depth $H$ is negligible relative to the growing diffusion lengths. Salt obeys the diffusion equation with diffusivity $D$. Withthe self-similar concentration profile in the liquid isIt has $C(h,t)=C_i$ and tends to $C_0$ in the far field.

Because the solid contains no salt, conservation of solute at the moving boundary requires the rejected solute flux from Fick's first law to equal the rate at which the interface sweeps up salt:Substituting the similarity profile givesThis relation determines $C_i$ if $\lambda$ is already known. More generally it must be solved together with the thermal Stefan condition and the interfacial phase-equilibrium relation $T_i=-m{C}_i$ from the liquidus. For $ϵ\ll 1$, salt occupies a much thinner boundary layer than heat and $C_i$ can be much larger than $C_0$.

Before the ice reaches the salty layer, the state moves vertically in a temperature--concentration phase diagram along $C=0$: cooling crosses the pure-water melting point and produces pure ice. After the encounter, the liquid state starts at $(T_{\infty },C_0)$ and approaches the interface stateon the liquidus $T_L=-mC$. The solid remains on the $C=0$ axis.

When $ϵ\ll 1$, the solutal diffusion length $\sqrt{Dt}$ is much smaller than the thermal diffusion length $\sqrt{\kappa t}$. Immediately ahead of the interface, rejected salt makes $C$ fall steeply from $C_i$ to $C_0$, so the local liquidus $T_L=-mC$ rises steeply with depth while the actual temperature changes comparatively slowly. Constitutional supercooling occurs whereverIts local onset criterion at the interface isThe supercooled zone lies in the thin salty boundary layer immediately ahead of the planar front. A small forward protrusion then enters liquid that is already below its local freezing point, so the planar interface is susceptible to a morphological instability.

The shallow-ice approximation makes the pressure hydrostatic pressure,Horizontal Stokes flow balance and the stress-free boundary condition at $z=h$ givesoThe basal shear stress and the stated basal sliding law are thereforeIntegrating once more and imposing $u(0)=u_b$ yieldsFor the right half of an ice cap, $h_x<0$, so both basal sliding and internal deformation carry ice away from the ice divide.

The depth-integrated ice flux isLocal mass conservation with accumulation rate $a$ saysSubstituting the flux givesas required. Equivalently, after using $u_b=-\beta h_x$,

Use the dimensional analysisand define the dimensionless sliding parameterThe dimensionless flux and steady conservation law areOn the right half-cap, let the snowline be $X=X_s$ and the nose be $X=X_n$. The zero-flux condition at the ice divide gives $Q=X$ in the accumulation region. Continuity of $Q$ gives $Q=2X_s-X$ in the ablation region, and zero nose flux gives

Introduce the increasing functionSince ${\Phi }_X=-Q$, the ablation profile satisfying $H(X_n)=0$ isAt the snowline $H=1$, soIn the accumulation region,and the divide thickness $H_0$ is fixed byThese two implicit formulas give a continuous thickness and flux at the snowline. Reflection across $X=0$ gives the full two-dimensional ice cap.

The ratio of sliding flux to deformation flux at thickness $h$ isThus basal sliding dominates throughout the dynamically important thick part of the cap whenwith the divide thickness giving a more conservative choice if it is known. In this limit $s\gg 1$, and the divide-to-nose extent from part c becomesThe full width of the symmetric cap is $2x_n$.

When sliding dominates and the whole cap is below the snowline, its thickness obeysWith $\alpha =0$, this is a one-dimensional porous medium equation. Its mass-preserving Barenblatt solution isWithin its support,Now set $h=h_{ss}+ct$. Since $h_{ss}$ solves the unforced equation,Matching the uniform ablation term $-\alpha$ givesThe similarity solution has a virtual origin. Writing $s=t+t_0$ gives the exact decaying familyIf its initial centre thickness is $h^{\ast }$, then $C=(h^{\ast }+3\alpha t_0/4)t_0^{1/3}$ and the exact extinction time isIts value depends on the initial cap width through the virtual age $t_0$, but its dimensional scale is unambiguouslyThe additive melting correction alone has the corresponding time $4h^{\ast }/(3\alpha )$.

For $\mathbf{u}=\nabla \times (\psi \hat{\mathbf{y}})$,Hence $u_x+w_z=-{\psi }_{zx}+{\psi }_{xz}=0$, so mass conservation holds identically. Removing the hydrostatic part by writingand taking the curl of the Stokes flow equation gives the biharmonic equation

For a Fourier mode proportional to $e^{ikx}$ with $k>0$, decay as $z\to -\infty$ selectsThe velocity and pressure fields are thereforeDirect substitution verifies the momentum equation.

For the surface $F(x,z,t)=z-\zeta (x,t)=0$, an outward normal vector and positively oriented tangent vector areWith $\zeta =ϵe^{ikx}$ and $|kϵ|\ll 1$, linearization gives

To first order, vanishing tangential velocity on the elastic plate is simply $u(0)=0$. In the mode from part a this givesConsequentlyso $u(0)=0$, $w(0)=ikA{e}^{ikx}$, and $w_z(0)=0$. The dynamic pressure at the surface is

Evaluating the hydrostatic pressure at $z=\zeta$ contributes the restoring normal stress $\rho g\zeta$. Since ${\nabla }^4\zeta =k^4\zeta$ for this mode, the linearized normal-stress balance isIt follows thatThe two restoring effects are buoyancy and the plate's bending stiffness.

The kinematic boundary condition is ${\zeta }_t=w$ at $z=0$. Part c therefore givesand henceIn terms of wavelength $\lambda =2\pi /k$,Without elasticity, $\tau \sim 2\mu k/(\rho g)$: shorter wavelengths penetrate less deeply but relax more slowly because their viscous gradients are larger. At short wavelength, elastic bending dominates andso the strong $B{k}^4$ restoring stress makes very short waves relax rapidly. The graph therefore tends to zero at both wavelength extremes and has one maximum. Differentiating with respect to $k$ places that maximum atThis exponential relaxation is the simplest half-space model of postglacial rebound and glacial isostatic adjustment.

For finite mantle depth, retain all four vertical solutions of the biharmonic equation,rather than discarding the terms that grow as $z\to -\infty$. Apply the plate conditions at $z=0$ and matching conditions at $z=-H$. At the mantle--core interface these include continuity of normal velocity and normal stress, an appropriate tangential-stress condition for the liquid core, and a kinematic boundary condition for interface displacement. The normal-stress balance gains the stable buoyancy termbecause ${\rho }_c>\rho$.

The resulting homogeneous linear system for the mode amplitudes has a nontrivial solution only when its determinant vanishes; that condition replaces the half-space decay rate. Modes with $kH\gg 1$ decay before sensing the core and recover the half-space result. Modes with $kH\lesssim 1$ involve the full mantle depth, feel the basal density contrast and core mobility, and acquire a modified relaxation time. This is the planar finite-depth extension of flexural isostasy.