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

Take $z$ upwards and write $\zeta =z+Wt$, so that the unperturbed interface is $\zeta =0$, with carbon dioxide above and water below. Use the sharp-interface, incompressible flow version of Darcy's law, with constant porosity and permeability of a porous medium, no capillary pressure, and effectively unbounded layers. The distinction between Darcy flux and interface speed matters: the base flux is $q_z=-\varphi W$, since the printed $W$ is the speed of the interface. Thus the base pressure gradients satisfyThe mass density ordering is stabilizing, whereas injection of the less viscous fluid into the more viscous fluid produces the Saffman–Taylor instability.

For linear stability analysis, displace the interface by $\eta =\hat{\eta }{e}^{\sigma t+i\alpha x}$. Taking the divergence of Darcy's law shows that the perturbation pressure satisfies the Laplace equation. Decay away from the interface gives, for $a_{\alpha }=|\alpha |>0$,The kinematic boundary condition and continuity of normal Darcy flux giveMeanwhile pressure continuity, evaluated on the displaced interface, givesEliminating $A,B$ yields the buoyancy-modified Darcy fingering dispersion relation:Every nonzero wavenumber decays when $W<W_c$, grows when $W>W_c$, and is neutral at equality, whereThe zero wavenumber is a uniform displacement and is neutral. Without a short-scale regularization such as surface tension, this model has no fastest-growing finite wavenumber: on the unstable side its growth rate increases without bound as $|\alpha |$ increases. Finite layer depths would change the dispersion relation; none are specified here.

Taking $g=9.81\mathrm{m}{\mathrm{s}}^{-2}$ and a year of $365.25$ days, the data giveThe numerical paragraph supplies no porosity, so it cannot determine a unique numerical interface speed. If $W$ were instead intended to mean Darcy velocity, the numerical critical value would be $103.2\mathrm{m}\mathrm{y}{\mathrm{r}}^{-1}$, and the interface would move at $W/\varphi$. These two velocity conventions must not be conflated.

Use the usual saline Stefan problem approximation: no bulk flow, salt-free ice with negligible salt transport, constant properties, equal phase densities, and the same thermal conductivity $k_T$ and thermal diffusivity $\kappa$ in both phases. Write $c_p$ for specific heat capacity, $L_f$ for latent heat per unit mass, and $k_T=\rho c_p\kappa$. These thermal symmetries are needed for the arithmetic-mean interface temperature requested in the paper; unequal phase conductivities would give a weighted balance instead.

Set $\eta =x/(2\sqrt{\kappa t})$, $\xi =x/(2\sqrt{Dt})$, $ϵ=\sqrt{D/\kappa }$ and $a(t)=2\lambda \sqrt{Dt}$. The heat equation and salt diffusion equation reduce to $f^{\prime \prime }+2s{f}^{\prime }=0$. Their similarity solutions, in terms of the complementary error function, areThese have the required interface values and far-field limits. At each fixed $x\ne 0$ they recover the initial data as $t\downarrow 0$. The liquidus condition is $T_i=-m{C}_i$.

Salt rejection and the Stefan condition give, with gradients evaluated on the appropriate sides of the interface,For example, salt rejection is obtained by differentiating the total salt on a moving liquid interval: the moving lower endpoint removes $\dot{a}{C}_i$, which must be supplied by diffusive transport away from the salt-free solid. Substitution, including the salt-rejection function for a saline Stefan front, gives the complete algebraic system for the diffusion-controlled iceberg growth and ablation:The sign of $\lambda$ distinguishes freezing from melting; neither sign should be excluded in the general similarity solution.

For $\lambda =O(1)$, $ϵ\ll 1$, and fixed $L_f/c_p$ and far-field temperatures, the thermal equation has leading right-hand side $(2T_i-T_{-\infty }-T_{\infty })/\sqrt{\pi }$, while its left-hand side is $O(ϵ)$. HenceThis is a leading-order balance, not an exact cancellation of latent heat. The resulting $C_i=-T_i/m$ and the first algebraic equation determine the leading $\lambda$. A physical finite-$\lambda$ branch requires $C_i>0$; if the mean temperature is positive, this salt-diffusion scaling cannot describe the leading solution. Likewise, a latent-to-sensible heat ratio diverging as $1/ϵ$ changes the leading thermal balance.

Define ${\theta }_{-}=-m{C}_0-T_{-\infty }$ and ${\theta }_{+}=T_{\infty }+m{C}_0$. Put $T_0=-m{C}_0$ and $\Delta \theta ={\theta }_{-}-{\theta }_{+}$. Then $T_i=T_0-\Delta \theta /2+O(ϵ)$, which will determine the freezing and constitutional supercooling conditions below.

Let $E(\lambda )=\sqrt{\pi }\lambda e^{{\lambda }^2}\mathrm{erfc}(\lambda )$. The salt equation in the saline Stefan problem gives $C_i-C_0=E(\lambda )C_i$. Since the complementary error function is positive on the real axis and $C_i>0$, the signs of $C_i-C_0$ and $\lambda$ agree. Thus freezing is equivalent to $C_i>C_0$, or $T_i<T_0=-m{C}_0$. Using the leading thermal balance,at leading order away from a vanishing temperature difference. In fact the exact onset is also ${\theta }_{-}={\theta }_{+}$: at $\lambda =0$ the salt balance forces $C_i=C_0$, and the Stefan condition then forces $T_{-\infty }+T_{\infty }=2T_0$. The solution with $\lambda =0$ has a stationary interface and balancing conductive heat fluxes, despite nontrivial diffusion of the initial temperature discontinuity. For completeness, the exact freezing criterion follows without relying on that local crossing. By the salt-rejection function for a saline Stefan front, $T_i(\lambda )$ decreases continuously from $T_0$ to $-\infty$ as positive $\lambda$ increases. If $T_{-\infty }<T_0$, let ${\lambda }_s$ be where $T_i=T_{-\infty }$. On $[0,{\lambda }_s]$, the solid-side thermal term decreases, the liquid-side term increases, and the latent-heat term increases: the factors $e^{-{ϵ}^2{\lambda }^2}/\mathrm{erfc}(-ϵ\lambda )$ and $e^{-{ϵ}^2{\lambda }^2}/\mathrm{erfc}(ϵ\lambda )$ decrease and increase respectively. These monotonicities follow from the positive integral representation of $e^{v^2}\mathrm{erfc}v$. The thermal residual starts at $({\theta }_{-}-{\theta }_{+})/\sqrt{\pi }$ and is negative at ${\lambda }_s$, so a positive root exists uniquely precisely when ${\theta }_{-}>{\theta }_{+}$. Beyond ${\lambda }_s$ both thermal contributions to freezing are negative, so no further root is possible. If $T_{-\infty }\ge T_0$, the same sign argument excludes every positive root.

For freezing, the liquid concentration decreases away from the interface, so the liquidus $T_L=-mC$ increases there. Constitutional supercooling means $T_l<T_L$ somewhere ahead of the interface. Since these temperatures agree at the interface, its local onset isFor the similarity solution and salt rejection condition, this is exactlyIt is also the global criterion for this freezing solution when $D<\kappa$ and $T_{\infty }>T_i$: the ratio $T_L^{\prime }(x)/T_l^{\prime }(x)$ is a positive constant times $\exp [-x^2(1/D-1/\kappa )/(4t)]$. For $x\ge a\ge 0$ this ratio decreases with $x$. If the ratio is at most one at the interface, $T_l-T_L$ cannot become negative; if it is greater than one there, a supercooled interval appears immediately ahead. Equality marks onset, not a finite supercooled interval.

For fixed positive $\Delta \theta$ and $\lambda =O(1)>0$, the left side of the inequality is $O(ϵ)$ and the right side is positive of order one, so freezing is constitutionally supercooled for sufficiently small $ϵ$. To locate the boundary between the two regimes, resolve the smaller scale $\Delta \theta =O(ϵ)$, $\lambda =O(ϵ)$. The salt equation and liquidus then giveHere $ϵ\lambda =O({ϵ}^2)$, so the thermal equation improves to $T_i=(T_{\infty }+T_{-\infty })/2+O({ϵ}^2)$; this justifies the second expansion including latent heat. Also $T_{\infty }-T_i={\theta }_{+}+O(ϵ)$ and the left side of the exact criterion is $ϵ{\theta }_{+}/\sqrt{\pi }+O({ϵ}^2)$. Therefore the small-diffusivity constitutional-supercooling threshold isThis printed condition is an asymptotic onset condition, rather than an exact finite-$ϵ$ inequality. In its $O({ϵ}^2)$ transition window the complete algebraic system and exact gradient criterion must be used.

Introduce specific heat capacity $c_p$, latent heat per unit mass $L_f$, and $\Delta T=T_m-T_s>0$. Write the prescribed convective heat flux as $F=K(\overline{T}-T_m)$, whereThere is a dimensional convention to specify: with $Q$ measured in $\mathrm{W}\mathrm{k}{\mathrm{g}}^{-1}$ and dimensionless $\lambda$, this printed expression makes $K$ a heat transfer coefficient only if $\mu$ is kinematic viscosity. If $\mu$ denotes dynamic viscosity, the numerator needs ${\rho }_0^{2}g\alpha Q$ instead. The question says only “viscosity”; use the consistent kinematic interpretation, or equivalently the corrected dynamic formula. The following balances apply to either convention once $K$ has units $\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$.

For a thin lava lake crust with a linear temperature profile, Fourier's law gives the outward conductive flux $k\Delta T/a$. The Stefan condition balances it against upward thermal convection and freezing, while the well-mixed liquid loses the convective flux and gains radiogenic heating. To leading order in $a/H$ the radiogenically heated lava-lake crust model isThe fixed surface temperature already incorporates radiative cooling, so no extra radiative loss is added to the interior balance. The linear conductive profile is interpreted as a quasi-steady approximation, so sensible-heat storage within the crust is not evolved. Thinness alone does not establish this approximation: rapid thermal adjustment in the crust is also needed, for example a small Stefan number during the initial conduction-driven growth. Replacing the liquid depth by $H-a$ and retaining the moving liquid-volume enthalpy would introduce thin-crust corrections: the requested constant-temperature branch refers to the leading model above. A transient thermal profile in a thick or newly nucleating crust would require the full Stefan problem, rather than this linear-profile closure.

Using the thermal diffusion time $H^2/\kappa$, defineThe three independent dimensionless parameters are the convective-conductive coefficient $B$, the heating parameter $R$, and the Stefan number $S$. For positive $B,R,S$ the evolution equations becomewhere prime denotes $d/d\tau$. Temperature has the exact solutionConsequently the unique constant-temperature initial condition isAt this temperature, the convective flux equals the column's radiogenic heat production. Starting with no crust, the constant-temperature lava-crust growth law is obtained by separating $h^{\prime }=S(1/h-R)$:The crust grows towards $h_{\ast }=1/R$, or $a_{\ast }=k\Delta T/({\rho }_{0}QH)$, where conduction carries all the radiogenic input. The logarithm makes reaching this equilibrium take infinite time. For an initial thickness $h_0>0$, subtract $[-R{h}_0-\log (1-R{h}_0)]/(S{R}^2)$ on the same side of the equilibrium; the corresponding antiderivative uses $\log |1-Rh|$ on either branch. The thin-crust description through equilibrium requires $R\gg 1$.

For general initial temperature, the full early and intermediate history is specified by the exact temperature above and one scalar ordinary differential equation. A particularly convenient integral formulation for an initially bare surface isThis removes the infinite initial derivative of $h$. The right-hand side points into $y>0$ at zero, and is decreasing in $y$ when ${\theta }_0\ge 0$; two nonnegative solutions therefore cannot separate, which gives uniqueness. Its solution, together with the temperature formula, determines the entire model evolution and can be computed without a phase-boundary search. The early expansion follows by substitution:The dimensional leading crust thickness is $a\sim \sqrt{2k\Delta Tt/({\rho }_0{L}_f)}$. Thus hotter liquid reduces the first correction to conductive crust growth, whereas colder liquid permits more rapid thickening.

If ${\theta }_0>{\theta }_{\ast }$, the temperature decreases and the instantaneous stationary lava-crust thickness $h_e(\tau )=1/[B\theta (\tau )]$ increases towards $h_{\ast }$. Starting from zero, the crust always remains below this rising curve: at a first contact $h^{\prime }=0$ but $h_e^{\prime }>0$, so it cannot cross from below. Hence hot-start crust grows monotonically. It also stays thinner than the constant-temperature solution: the early expansion puts it below, and at a proposed first crossing their derivative difference is $-SB({\theta }_0-{\theta }_{\ast })e^{-B\tau }<0$, preventing crossing from below. If $0\le {\theta }_0<{\theta }_{\ast }$, the temperature increases and the instantaneous stationary lava-crust thickness decreases. The crust stays thicker than the constant-temperature solution by scalar comparison, but it may cross the decreasing $h_e$, reach a maximum and subsequently melt back towards $h_{\ast }$. Such overshoot is possible, not inevitable; the two relevant relaxation rates are $B$ and $S{R}^2$. The numerical plot uses $B=4$, $R=20$, $S=0.1$, so crust adjustment is faster than temperature relaxation and the cool-start example does overshoot:

The curves use ${\theta }_0/{\theta }_{\ast }=0.5,1,1.5$ and initially $h=0$. They solve the complete reduced time-dependent equations, rather than replacing the moving stationary thickness by the actual crust. At late times, set $h=h_{\ast }+v$, $\delta ={\theta }_0-{\theta }_{\ast }$; linearization gives $v^{\prime }=-S{R}^{2}v-SB\delta e^{-B\tau }$. For $S{R}^2\ne B$,At equal rates the forced term is $-SB\delta \tau e^{-B\tau }$. In particular, slow warming with $B<S{R}^2$ approaches equilibrium from a thicker crust, explaining the plotted remelting. These long-time expressions describe the linearized response; nonlinear corrections can have their own faster exponential rates. All predictions are restricted to $a/H\ll 1$, nonnegative liquid superheat, and the assumed vigorous-convection flux law. If crust growth violates those conditions, the reduced model must be replaced rather than extrapolated to total freezing.

Let $z$ increase upwards, so the bed is at $z=-b(x)$, the surface is at $z=h(x,t)$, and $H=h+b$. In the shallow-shelf approximation the incompressible flow has nearly depth-uniform horizontal speed $u$ and vertical strain $w_z=-u_x$. The Newtonian fluid stress tensor, together with the hydrostatic approximation, then givesThe factor four includes both horizontal extension and the pressure correction required by vertical compression. Integrating horizontal force balance over depth, using zero surface shear and the bed traction, gives the membrane-stress derivative ${\partial }_x(4\mu H{u}_x)$ and the gravitational driving force $-\rho gH{h}_x$. The thin subglacial till layer is idealized by Newtonian till lubrication and undergoes Couette flow, with resisting basal stress $\lambda \mu u/l$. Hence force balance and mass conservation yieldSince $b$ is fixed, $H_t=h_t$. No accumulation or ablation is included.

At the grounding line, the first boundary condition is ice-sheet flotation:For an unbuttressed ice shelf, the second is the extensional stress required to balance the difference between the integrated ice and seawater hydrostatic pressures. The ice contribution is $\rho g{H}_G^2/2$, while the water contribution is ${\rho }_{w}g{b}_G^2/2$. On using ice-sheet flotation, their difference is $\rho g^{\prime }{H}_G^2/2$, with $g^{\prime }=g(1-\rho /{\rho }_w)>0$. Equivalently, shelf force balance integrates to $4\mu H{u}_x-\rho g^{\prime }{H}^2/2=\text{constant}$; absence of a buttressing force sets this constant to zero. Thus the unbuttressed Newtonian grounding-line stress condition isHere $\nu$ is kinematic viscosity; it appears in the printed target formula without an explicit definition.

The ratio of longitudinal stress divergence to till drag is of order $4lH/(\lambda L^2)$. The given small-parameter limit therefore yields the friction-dominated bulk relationThe coefficient depends on the basal lubrication, and $h_x<0$ gives seaward motion. This is a bulk reduction: longitudinal stress has been neglected in the differential equation but its boundary traction remains specified. A complete uniformly valid approximation close to the grounding line could require a membrane-stress boundary layer. In the remainder, use the pointwise friction closure and retain the shelf stress condition as the boundary data, as requested in the paper.

Differentiate the ice-sheet flotation condition along the moving grounding line:Now $H_x=h_x+\alpha$. Substitution of the friction-dominated speed and the shelf stress boundary condition gives the friction-dominated grounding-line evolution law:The coefficient of ${\dot{x}}_G$ is the spatial derivative of the flotation deficit, up to sign. If it vanishes, the implicit equation remains the correct compatibility condition but division to determine a finite speed is invalid; tangential contact needs separate analysis.

In a steady state, volume flux per unit width is constant: $Hu=q_0$. The friction relation givesDifferentiating $Hu=q_0$ and using $u_x=g^{\prime }H/(8\nu )$ at the grounding line givesPut $\beta =l/\lambda$, $\mathcal{H}=H/\beta$ and $\mathcal{Q}=\nu q_0/(g{\beta }^3)$. Substitution and cancellation of $g{\beta }^5/\nu$ yield the steady friction-dominated grounding-line thickness relation:For $q_0>0$ and a bed deepening seaward ($\alpha \ge 0$), the polynomial $(g^{\prime }/8g){\mathcal{H}}^5+\alpha \mathcal{Q}{\mathcal{H}}^2-{\mathcal{Q}}^2$ is strictly increasing for $\mathcal{H}>0$, starts negative and tends to infinity, so it has exactly one positive root. At that root $\mathcal{Q}>\alpha {\mathcal{H}}^2$ and $H_x<0$, as required by steady outward extension. The relation belongs to the specified friction closure; it is not a flux law for arbitrary till rheology or arbitrary membrane-stress matching.