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

Use the Dupuit approximation: the saturated region is shallow enough that hydrostatic pressure is $p=\rho g(h-z)$ and flow is predominantly horizontal. Darcy's law then gives the horizontal Darcy velocity $u_x=-u_b(1+\beta z)h_x$, where $u_b=k_0\rho g/\mu$. Integrating through the saturated depth gives the volume flux per unit widthThe stored water volume per unit area is $\varphi h$, so mass conservation yields the unconfined aquifer with depth-dependent permeability equationDefine the positive river discharge by $Q(t)=-q(0,t)$. A useful check on all subsequent results is the integrated mass conservation law

In the steady state, mass conservation reduces to $q_x=R$. The zero-flux condition at the divide gives $q=R(x-L)$. Substituting the Darcy flux and integrating from the absorbing river boundary yieldsThus the nonnegative groundwater height is the unique nonnegative root ofThe height increases toward the divide, where $h_x=0$. The river receives all the rainfall, so its steady volume flux per unit width is $\boxed{Q=RL}$.

Before drainage reaches the divide, most of the aquifer fills locally to height $H=Rt/\varphi$. Only a growing boundary layer beside the river departs substantially from this height. Write $D=u_b\beta /2$.

At early times $H\ll 2/\beta$, the mobility is $u_{b}h$. The balance between storage and nonlinear diffusion givesSubstitution into the governing equation gives the forced filling similarity for power-law diffusion with exponent $m=1$:The positive outward volume flux per unit width isThe boundary slope is singular, but $F_e\sim (2c_e\eta {)}^{1/2}$ gives a finite flux.

At intermediate times $H\gg 2/\beta$, the large-depth mobility is $D{h}^2$. The corresponding similarity solution iswithConsequently,Here $F_i\sim (3c_i\eta {)}^{1/3}$ is the outer outlet profile. The height actually vanishes at the river, so an outlet layer in a deep unconfined aquifer restores the $u_{b}h$ mobility extremely close to the boundary while transmitting this same leading discharge.

The mobility changes when $H\sim 2/\beta$, whereas the divide first affects the deep filling solution when ${\ell }_i\sim L$. ThusA distinct intermediate regime requires $t_{\beta }\ll t_{\mathrm{fill}}$. If the divide is reached while the aquifer is still shallow, the early filling law instead crosses directly toward the steady state, at a time of order $\varphi L/\sqrt{u_{b}R}$. The faster intermediate discharge growth comes from increasing permeability of a porous medium as deeper parts of the aquifer become saturated; once the finite domain is felt, the discharge approaches $RL$.

After recharge stops, the finite interval sets the horizontal scale. Apply the separable draining profile for power-law diffusion to each limiting equationWith $\xi =x/L$ and virtual age $\tau =t+t_0$, the similarity solution isThe positive profile fixes $c_m={\lim }_{\xi \downarrow 0}{F}_m^m{F}_m^{\prime }={\int }_0^1{F}_{m}d\xi$. The discharge follows either from the boundary Darcy flux or integrated mass conservation:In the deep regime this givesIn the final shallow regime,The deep profile contains the thin outlet layer in a deep unconfined aquifer described above. The whole aquifer eventually becomes shallow, and its nonzero basal permeability controls the ultimate ${\tau }^{-2}$ discharge.

Estimating the mobility crossover by a typical height $2/\beta$ givesUsing the deep profile's height at the divide makes this ${\tau }_{\mathrm{cross}}\simeq \varphi \beta L^2{F}_2(1{)}^2/(4u_b)$. These are regime-dependent similarity approximations: the virtual origins are set by matching to the initial drainage and crossover, and a deep regime occurs only if the aquifer is initially sufficiently deep.

Let a representative fixed volume have mixture mass density $\bar{\rho }={\rho }_s\phi +{\rho }_l(1-\phi )$. Since the crystal framework is rigid and stationary, the liquid mass flux is ${\rho }_l\mathbf{u}$, where $\mathbf{u}$ is Darcy velocity. Mass conservation givesThe constant phase densities therefore imply the phase-change volume source in a rigid mushFor ice less dense than brine, increasing solid fraction creates an expansion flow; when the phase densities agree, this source vanishes.

In the heat-exchanger frame the crystal velocity is ${\mathbf{v}}_s=-V{\mathbf{e}}_z$, while the liquid velocity is ${\mathbf{v}}_l={\mathbf{v}}_s+\mathbf{u}/(1-\phi )$. The steady vertical mass flux is constant and equals $-{\rho }_{s}V$ in the fully solid eutectic material. HenceThis gives the density-change flow in a pulled mushy layerThe laboratory mixture volume flux is $-V[r+(1-r)\phi ]{\mathbf{e}}_z$. In the pure liquid $\mathbf{u}=(1-r)V{\mathbf{e}}_z$, while in the fully solid region $\mathbf{u}=0$ and the material moves at $-V{\mathbf{e}}_z$. The liquid velocity formula applies wherever liquid exists; at solid fraction one there is no liquid phase to assign a velocity.

For a steady field in the apparatus frame, the time derivative in the translating crystal frame is $-V{\partial }_z$. The given solute equation therefore becomesSubstituting the density-change flow in a pulled mushy layer reduces this to $(1-\phi )C_z=C{\phi }_z$. Thus the solute conservation in a steadily pulled mush integrates to $(1-\phi )C=C_0$, using the zero solid fraction at the top boundary.

Set $\mathcal{C}=m{C}_0/\Delta T$. The liquidus gives $C=C_0(1-\theta /\mathcal{C})$, soIn particular, $\phi =0$ at $\theta =0$, and immediately above the eutectic temperature, where $\theta =-1$, the mush has $\phi =1/(\mathcal{C}+1)$. It becomes fully solid across the eutectic front below.

Use $\zeta =Vz/\kappa$, $H=Vh/\kappa$, and $S=L/(c_p\Delta T)$, the latent-to-sensible heat ratio. Transforming the heat equation to the fixed apparatus frame and substituting Darcy flux gives the exact steady mush equationFor $\mathcal{C}\gg 1$, the solid fraction is $\phi =-\theta /\mathcal{C}+O({\mathcal{C}}^{-2})$. Since $S/\mathcal{C}=O(1)$, the latent term remains leading order, whereas the correction $(1-r)\phi {\theta }_{\zeta }$ is small. Define $\Omega =1+S/(r\mathcal{C})$. The large-concentration thermal profile of a pulled mush then obeysWith $\theta (H)=0$ and $\theta \to {\theta }_{\infty }$ in the liquid,The thermal conductivity agrees on both sides and the solid fraction vanishes at the mush–liquid interface, so there is no jump in latent production there: continuity of heat flux gives ${\theta }_m^{\prime }(H)={\theta }_l^{\prime }(H)=r{\theta }_{\infty }$. ThusThe boundary value ${\theta }_m(0)=-1$ determines $H$, as calculated next. The crystal fraction is obtained from the previous part's formula, with the temperature field understood to this leading asymptotic accuracy.

Impose the eutectic temperature at $\zeta =0$ in the large-concentration thermal profile of a pulled mush:ThereforeThe thickness is proportional to thermal diffusivity divided by pulling speed. Latent-heat release modifies the exponential decay rate through $\Omega$, while the far-liquid temperature fixes how far the profile can remain below the liquidus.

Use the common-density, equal-volume-heat-capacity model implicit in the stated speed law. Let the liquid fraction be $\varphi$ ahead of the melting front and $\varphi +\Delta \varphi$ behind it. Melting adds exactly the liquid needed to fill the newly created pores, so total mass conservation makes Darcy flux continuous. Consequently,Here $u$ denotes Darcy velocity and $v$ the pore-liquid velocity in the stationary rock frame. Treating the pore-space increase as storage without its simultaneous melting source would incorrectly change the Darcy flux.

Relative to the cold unmelted material, the bulk enthalpy increase behind the front is $\rho c_p\Delta T+\rho L\Delta \varphi$: all phases gain sensible heat and the melted ice consumes latent heat. The advective heat-flux difference is $\rho c_{p}U\Delta T$. Applying the Rankine-Hugoniot condition to this energy balance therefore givesso the advection-driven melting front in a porous matrix moves atThe advected latent contribution of the liquid is the same on both sides and cancels; the denominator measures sensible heating plus phase-change energy per unit bulk volume.

Put $k_a=k$ ahead of the front and $k_b=k+\Delta k$ behind. With viscosity $\mu$, Darcy's law gives the planar pressure gradients $p_{a0}^{\prime }=-\mu U/k_a$ and $p_{b0}^{\prime }=-\mu U/k_b$. Write the displacement as $\eta e^{\sigma t+i\alpha y}$ with $\alpha >0$, and take the unperturbed interface to be $x=0$ in its translating frame.

The liquid flux is divergence-free away from the interface, so the pressure perturbations are harmonic. The decaying normal modes areExpanding pressure continuity at the displaced interface and matching normal Darcy flux givesHenceThe common normal-flux perturbation is $\delta u_n=-k_b\alpha B/\mu =\alpha U\Delta k\eta /(k_a+k_b)$. The local enthalpy balance converts it into a front-speed perturbation, $\sigma \eta =\delta u_n/(1+S\Delta \varphi )$. Thus the melting-front instability due to permeability contrast hasFor a signed transverse wavenumber, replace $\alpha$ by $|\alpha |$.

Melting raises permeability of a porous medium. An advancing protrusion therefore offers a less resistive path to hot liquid, attracting larger Darcy flux and receiving more heat. The resulting faster melting reinforces the protrusion: this is a reactive infiltration instability. A permeability decrease would reverse this feedback, while zero contrast gives no growth in this ideal model.

For positive contrast, the growth rate is proportional to transverse wavenumber, so the model has no finite fastest-growing wavelength. The latent-to-sensible heat ratio slows the front but supplies no short-wave cutoff. Finite thermal transport and a resolved melting zone can introduce a length scale; pore geometry eventually limits the continuum description. If interface curvature changes the equilibrium melting temperature through the Gibbs--Thomson relation, that can also oppose short-wave corrugations. The modified dispersion relation must include such effects before a preferred wavelength can be calculated; a length scale alone does not guarantee that every added mechanism selects one.

Let $\ell =(H_m-H_0)/\sin \alpha$ be the along-slope distance to the snowline, so $H(x)=H_m-x\sin \alpha$ and $a(x)=A(1-x/\ell )$. Define $D=g\sin \alpha /(3\nu )$.

Treat ice as an incompressible Newtonian fluid, neglect inertia, and use lubrication theory with thickness measured normal to the slope. The bed has a no-slip boundary condition, the free surface has zero tangential stress, and normal pressure is hydrostatic. The small thickness slope allows its pressure-gradient contribution to be neglected against gravity along the mountain. With normal coordinate $y$, the tangential equation and boundary conditions areIntegration givesThe horizontal radius of a ring is $x\cos \alpha$, so mass conservation gives the gravity-driven ice flow on a conical slope equationNegative accumulation is ice ablation and applies only where ice exists; the ice-free region has $h=0$.

In a steady state, regularity and zero total flux at the apex require $xD{h}^3\to 0$ as $x\to 0$. Integrating givesRequiring a continuous zero-thickness steady terminus yields the steady conical ice cap with a linear accumulation gradient:The terminus lies below the snowline, allowing the ablation region to balance snowfall. The maximum thickness occurs at $x=3\ell /4$.

The volume follows from integrating the ring areas. With $s=x/\ell$,Substituting $D$ and $\ell =\Delta H/\sin \alpha$ givesThe ideal outer profile has steep slopes very close to the apex and terminus. Those small regions require local corrections to the assumed slope balance, while the bulk profile and leading volume follow from the stated approximation.

With both ice accumulation and ice ablation removed, the equation is $h_t+x^{-1}(xD{h}^3{)}_x=0$ and the conserved volume is $V_0=2\pi \cos \alpha M$, where $M={\int }_0^{x_N}xhdx$.

If the thickness and extent scales are $H(t)$ and $X(t)$, mass conservation gives $H{X}^2\sim M$, while the flow equation gives $H/t\sim D{H}^3/X$. Therefore $X\propto t^{1/5}$ and $H\propto t^{-2/5}$. Set $h={\tau }^{-2/5}f(\xi )$, $\xi =x{\tau }^{-1/5}$, with a possible virtual time origin $\tau =t+t_0$. The similarity solution satisfiesIntegration and regular zero total flux at the apex give $D\xi f^3={\xi }^{2}f/5$, hence the positive profile iswith dry bed beyond the front. Volume normalization givesso the volume-conserving conical ice-current similarity solution hasThe terminus has finite thickness and is a shock wave in the gravity-only kinematic wave equation. The Rankine-Hugoniot condition gives ${\dot{x}}_N=q_N/h_N=D{h}_N^2=x_N/(5\tau )$, exactly agreeing with the similarity extent. The characteristic speed behind the front is $3D{h}_N^2$, larger than its speed, while the dry-bed characteristic speed is zero, so the front is compressive and gives an entropy solution.

This solution describes the long-time spreading, rather than exactly matching the earlier steady profile at the instant snowfall stops. The initial transient can be described by the characteristic transformation for conical ice drainage: with starting point $s$, put $K(s)=s^{1/3}{h}_s(s)$; thenwhere characteristics remain smooth. Subsequent crossings are resolved by the same conservation and entropy conditions. Restoring the neglected local pressure gradient would smooth the idealized front.

Initially, away from the apex and snowline, transport is small and local ice accumulation dominates:The early ice volume is $V(t)\simeq \pi A{\ell }^2\cos \alpha t/3$. Balancing radial transport $D(At{)}^3/\ell$ with $A$ identifies the global transition scalesThere is also an earlier local adjustment near the apex: the factor $1/x$ makes the same flux significant at distance $x_a\sim D{A}^2{t}^3$. The apex filling layer of a conical ice cap is described byAt the apex, $F\sim (\eta /2{)}^{1/3}$ gives the local steady profile $h\sim (Ax/(2D){)}^{1/3}$. Thus gravity-driven redistribution reaches a fixed position $x$ on a time of order $[x/(D{A}^2){]}^{1/3}$, and its region grows out to the snowline when $t\sim t_{\ast }$.

At times comparable to $t_{\ast }$, ice transport becomes important across the mountain's accumulation region, carries ice below the snowline, and establishes an ablation zone. The cap then approaches the steady conical ice cap with a linear accumulation gradient, with $x_N\to 3\ell /2$. For a moving ice-covered disk, mass conservation and its front condition giveThe total input vanishes at $x_N=3\ell /2$, explaining the final extent. The narrow apex and front regions are governed by local adjustments to the bulk approximation; they do not create a further distinct large-scale time regime in this scaling description.