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

Take $x$ in the injection direction and $y$ along the flat interface. The specified $U$ must be a Darcy velocity, so the pore-fluid velocity and unperturbed interface speed are $U/\varphi$. Write the perturbed interface as $x=X(t)+\eta (y,t)$, with $\dot{X}=U/\varphi$, and choose the normal from fluid 1 into fluid 2. The background pressure gradients from Darcy law are ${\partial }_x{p}_i^{(0)}=-{\mu }_{i}U/k$.

By incompressibility and uniform permeability of a porous medium, pressure perturbations are harmonic functions. A transverse Fourier mode with positive wavenumber $\alpha$ has the decaying formsThe kinematic boundary condition and equal normal Darcy flux implyWithout capillarity the pressure is continuous at the displaced interface. Its linear perturbation is thereforeConsequently the planar viscous-fingering dispersion relation reduces toThis is the Saffman–Taylor instability: less viscous fluid advances preferentially through protrusions. Without a short-wave regularizing mechanism the idealized rate has no finite maximum.

For the prescribed apparent capillary pressure, the chosen normal is $\mathbf{n}=(1,-{\eta }_y)/\sqrt{1+{\eta }_y^2}$, so $\nabla \cdot \mathbf{n}=-{\eta }_{yy}+O({\eta }^2)={\alpha }^2\eta +O({\eta }^2)$. Subtract the flat-front pressure jump $\gamma (1+\beta X){\kappa }_0$. Linearization of the jump at $x=X+\eta$ givesThere is no first-order product of the surface-tension perturbation and the perturbed macroscopic curvature. ThusUsing the capillary number convention $\mathrm{Ca}={\mu }_{2}U/\gamma$ and ${\mu }_2/({\mu }_1+{\mu }_2)=(1+M)/2$, this isIf $x$ is measured from the unperturbed front and $\gamma$ denotes its local apparent tension, set $X=0$ to obtain the usual fixed-coefficient expression. For a physically fixed spatial gradient, $X=Ut/\varphi$ and the displayed rate is an instantaneous, generally time-dependent growth rate; the amplitude solves $\dot{a}=\sigma (\alpha ;X(t))a$, rather than one constant exponential law for the entire displacement. The parameter $\beta$ has dimensions of inverse length; the absolute tension gradient is $\gamma \beta$.

Assume $1+\beta X>0$, as required for positive apparent surface tension, and define $A=({\mu }_2-{\mu }_1)U-k\gamma \beta {\kappa }_0$. If $A>0$, unstable modes satisfy $0<{\alpha }^2<A/[k\gamma (1+\beta X)]$. Differentiating the cubic dispersion relation gives the fastest-growing modeWhen $A\le 0$, every nonzero mode decays, with a neutral zero-wavenumber translation in the infinite-domain limit. For $\beta {\kappa }_0>0$, surface-tension-gradient stabilization of viscous fingering suppresses the instability entirely at positive injection speeds satisfyingIf $\beta {\kappa }_0\le 0$, no positive injection speed eliminates every long-wave unstable mode in an unbounded interface. The prescribed wetting variation is treated as an effective normal pressure law for porous-media flow.

Let $w$ be downward Darcy flux, measured per unit total horizontal area. Neglect the resistance of the displaced phase, capillary pressure, and horizontal subsurface flow. Taking atmospheric pressure as zero, the upper end of the saturated column has pressure $\rho gh$ and its advancing lower end has pressure zero. Across a column of depth $l$, Darcy law therefore givesThe pressure head from the overlying layer supplies $h/l$, and gravity supplies one. Because the new liquid fills a pore volume $\varphi dl$ per unit area, the vertical imbibition under a draining liquid layer obeys $\varphi \dot{l}=w$. Volume conservation gives $h+\varphi l=h_0$, so, until the surface layer is exhausted,The initial condition is understood as a limiting solution of this singular equation. For $0<\varphi <1$, separation of variables givesAt small penetration, expansion of the logarithm yields $t\sim \varphi \nu l^2/(2gk{h}_0)$, henceA sufficient dimensional early-time condition is $t\ll \varphi \nu h_0/(gk)$, for which $l\ll h_0$. The printed $t\ll 1$ needs a chosen time unit or nondimensionalization. In the limiting case $\varphi =1$, $l^2=2gk{h}_{0}t/\nu$ holds exactly until $h=0$. The divergent initial Darcy flux is an idealization of sudden contact with a dry continuum; it cannot describe penetration below the pore scale.

For the spreading surface layer, measure $z$ upwards from the substrate. A thin, slow gravity current has hydrostatic pressure $p=\rho g(h-z)$ and a leading horizontal viscous balance $\nu u_{zz}=g{h}_x$. A no-slip boundary condition at $z=0$ and a stress-free boundary condition at $z=h$ giveThis lubrication gravity-current flux drives horizontal spreading down the surface-pressure gradient. Surface volume conservation is $h_t+q_x=-w$, while conservation in the pores is $\varphi l_t=w$. Combining these balances gives the deep-substrate drainage of a gravity current modelDrainage transfers liquid from the surface layer into the substrate. Indeed $(h+\varphi l{)}_t+q_x=0$.

For definiteness take a one-sided current on $0<x<x_N(t)$, supplied at $x=0$. The imposed volume is per unit span, so $V_0$ has dimensions of area. The boundary conditions and volume conservation constraint areThe source starts from zero injected volume. Initially dry points acquire $l>0$ only when the advancing current reaches them; the formula $h/l$ is not applied outside that region. At the advancing nose these are limiting conditions for a singular drainage equation. For a source feeding two identical sides, each side has half the volume and inlet flux; only the corresponding numerical prefactors change. The integral and inlet-flux conditions are equivalent once the local balance and nose conditions hold.

The physical approximations are lubrication theory, negligible inertia, slowly varying surface height, small pore-scale Reynolds number, and negligible capillary entry pressure. Treating the substrate as vertically draining columns also requires its lateral flow to be negligible and its depth to exceed the penetration depth. The model is valid while surface liquid remains available above each draining column.

For a similarity solution, let $h$ and $l$ scale as $t^a$, and let the current length scale as $t^b$. Since $h/l$ then has no time dependence, drainage requires $a-1=0$. Balancing $h_t$ with the horizontal spreading term gives $a-1=4a-2b$, so $a=1$ and $b=2$. These exponents also give injected volume $t^{a+b}=t^3$, explaining the particular supply law.

Define scaleswhich satisfy $CD=V_0/{\tau }^3$ and $D^2=(g/\nu )C^3$. The cubic-input similarity for a draining gravity current isSubstitution into the two conservation equations giveswhereThus only $K$ and $\varphi$ remain in the dimensionless problem. Its constraints areIntegrating the sum of the two ordinary differential equations verifies that the inlet flux and integral normalization agree. The dimensionless endpoint ${\xi }_N={\xi }_N(K,\varphi )$ is the undetermined multiplicative constant in the length law; solving for it is not required.

For an axisymmetric thin elastic plate, define ${\nabla }_r^2=r^{-1}{\partial }_r(r{\partial }_r)$. Uniform excess pressure and bending stiffness $B$ give $B{\nabla }_r^{4}h=p$. The roof is attached to undeformed rock at the front, giving clamped boundary conditions $h(R)=h_r(R)=0$, and regularity excludes singular displacement or curvature at $r=0$. A regular particular solution is $p{r}^4/(64B)$; adding a constant and a multiple of $r^2$ and enforcing the two edge conditions gives the axisymmetric clamped-plate deflectionIntegrating the magma intrusion volume,therefore yieldsHere the axisymmetric phrase requires the radial plate biharmonic operator; a Cartesian one-dimensional beam equation would not give the same pressure coefficient. Hydrostatic contributions are excluded as instructed.

For the early front, let $U_p=\dot{R}$ be the leading propagation speed. The fluid and fracture radii differ only within a local region of length $l_p\ll R$. In lubrication theory, the elastic pressure gradient near the front is of order $B{h}_f/l_p^5$, and the volume flux through the gap is of order $B{h}_f^4/(\mu l_p^5)$. In a frame following the fluid front that flux is of order $U_p{h}_f$. The viscous peeling of an elastic plate law is thereforeMatch the peeling region to the interior curvature: $h_f/l_p^2\sim {\kappa }_i$. The stated vapour-gap relation $h_f\simeq p_T{l}_p^4/(24B)$, with $p_T>0$ interpreted as the relevant pressure magnitude, then givesSubstitution gives the vapour-tip peeling of a magma intrusion speedWith constant supplied volume flux $Q$ and negligible initial intrusion volume, $V=Qt$. Integrating $R^{14}\dot{R}\sim B^{3/2}{Q}^{7/2}{t}^{7/2}/(\mu p_T^{1/2})$ gives $R^{15}\sim B^{3/2}{Q}^{7/2}{t}^{9/2}/(\mu p_T^{1/2})$, orOrder-one coefficients require the full local peeling profile. This early regime is a front-controlled asymptotic regime, not an assertion that a continuum peeling zone remains thin at arbitrarily small $t$.

When fracture toughness instead imposes fixed ${\kappa }_f>0$, equating it to the outer edge curvature gives the toughness-controlled magma intrusion lawIts pressure can be written in either useful formThe pressure decreases as the intrusion grows on this advancing branch.

If $V_{r0}$ is the reservoir volume before any magma was transferred to the initially negligible intrusion, volume conservation and negligible conduit storage give $V_r=V_{r0}-V$. With the supplied pressure-volume law and conduit conductance, the reservoir-fed magma intrusion obeysIn this law $E$ is an effective reservoir pressure-volume stiffness, with dimensions of pressure per volume. If the stated initial reservoir volume refers instead to the start of the late regime, when the intrusion already has volume $V_i$, replace $V_{r0}$ by the total $V_T=V_{r0}+V_i$ throughout. An independent value of $V_i$ would then be needed.

For an advancing solution, inflow ceases when the reservoir and intrusion pressures equalize. Its final volume is the larger positive root in $0<V<V_{r0}$ ofEquivalently, $y=\sqrt{V_{\infty }}$ solves $y^3-V_{r0}y+\mathcal{A}/E=0$. The larger root is essential: $\sqrt{V}(V_{r0}-V)$ rises from zero to a maximum $2V_{r0}^{3/2}/(3\sqrt{3})$ at $V=V_{r0}/3$ and then decreases. A nonempty interval of positive inflow exists only ifThere are then two positive equilibria $V_{-}<V_{r0}/3<V_{+}$. At equilibrium the derivative of the volume ordinary differential equation is $\beta E(V_{r0}-3V)/(2V)$, so $V_{-}$ is unstable and $V_{+}$ has asymptotic stability. An initial advancing intrusion with $V_{-}<V_i<V_{+}$ approaches $V_{\infty }=V_{+}$. The equality case has a double equilibrium and no positive-inflow interval; above the threshold there is no advancing equilibrium branch of this model.

The late formula $p\propto V^{-1/2}$ is singular at $V=0$ and cannot nucleate the intrusion: the earlier propagation regime must supply its finite seed. Nor does this fixed-edge-curvature propagation law prescribe how an irreversible fracture closes under negative inflow. The final-size answer is conditional on reaching the advancing late branch. In particular the reservoir is generally not emptied, and $\beta$ controls the approach rate rather than the equilibrium size.

Use a quasistatic surface energy balance with negligible heat storage in the surface or ice, a linear ice temperature profile, and a fixed freezing temperature at the ice-ocean interface. These are the assumptions behind the single-category sea-ice thermodynamic model used below. Write $\theta =t/\tau$, $F_O={\lambda }_O(T_O-T_m)>0$ andHere ${\lambda }_A$ and ${\lambda }_O$ are heat transfer coefficients, $k$ is thermal conductivity, $\kappa =k/(\rho c_p)$ is thermal diffusivity, and $\mathcal{S}$ is the inverse Stefan number. Here $L$ is latent heat per unit mass. The thermal $k$ denotes conductivity; the previous questions use $k$ for permeability.

While there is open water, the upward ocean heat supply balances the loss to air: ${\lambda }_O(T_O-T_s)={\lambda }_A(T_s-T_A)$. ThusFreezing begins when this surface reaches $T_m$ on the cooling part of the annual cycle, which requires $r<1$ and givesFor $r\ge 1$ the surface never cools strictly below freezing and no positive ice thickness develops in this model. The limit $r\to 0$ gives $t_1\to \pi \tau /2$.

With ice present and a cold surface, thermal conduction and Newton's law of cooling give $k(T_m-T_s)/h={\lambda }_A(T_s-T_A)$. The surface temperature isDuring surface melting the temperature is capped at $T_s=T_m$; the ice is not permitted to follow the above formula to a temperature higher than its melting point. The upward conductive heat flux during the cold phase isThe Stefan condition at the ocean interface gives $\rho L\dot{h}=F_C-F_O$. Under the permitted approximation of negligible ocean heat during growth,Starting from $h(t_1)=0$, integration yields the implicit growth lawThe leading maximum occurs when the air returns to freezing from below:The equality for $h_m$ is within the ocean-heat-neglected growth model, with the freezing time set by the open-water balance. If $aH\gg 1$, conductive resistance through the ice dominates the atmospheric boundary resistance at maximum thickness. ConsequentlyThis reproduces the displayed large-conductance limit. A sufficient condition for the quasistatic temperature assumption is $h_m^2/\kappa \ll \tau$; large $\mathcal{S}$ ensures this in the thick-ice scaling.

The approximation concerning ocean heat needs care. Its ratio to atmospheric heat removal during the cold phase isIt is uniformly small away from the beginnings and ends of winter providedFor the thick-ice limit this requires both $r\ll 1$ and $raH\ll 1$. This also makes the integrated ocean heat small relative to the latent heat of the maximum ice column. It is not possible to neglect $F_O$ pointwise throughout the closed growth interval: at the actual freezing onset it balances the atmospheric removal, and at the true maximum the growth rate is zero because the two fluxes again balance. The ocean-heat correction to seasonal ice growth therefore gives narrow endpoint corrections. In the full model the maximum is slightly earlier than $3\pi \tau /2$, withThus $t_2=3\pi \tau /2-O(\tau r[1+a{h}_m])$ in the stated approximation, rather than an exact maximum time for finite ocean heat.

For $3\pi \tau /2\le t\le 2\pi \tau$, while ice survives, both the warm air and ocean melt it. The surface is at $T_m$, and the total Stefan condition isIntegrating from the leading maximum givesuntil its first zero. The physical thickness is then zero for the rest of this warm interval and the open-water temperature applies. Therefore at the end of the first year,Positive ice at that time is equivalent, in the specified approximation, toThis establishes the algebraic criterion printed in the paper.

There is, however, an actual qualification to the paper's word “perennial.” This inequality tests peak-temperature sea-ice survival, since $t=2\pi \tau$ is the next atmospheric temperature maximum. Surface melting continues until at least $t=5\pi \tau /2$, so positive thickness at the maximum does not guarantee survival throughout the summer. In the same approximation, survival through that full warm half-cycle would require the stronger inequalityIn fact, starting ice-free under the specified zero-mean sinusoidal forcing, even the zero-ocean-heat upper bound on winter growth cannot meet it. In that limit $t_1=\pi \tau /2$ andsince $h_m^2/2>0$. The expression on the right is already the atmospheric melt over the entire subsequent warm half-cycle. Positive ocean heat can only reduce the winter ice and increase summer melt. Thus the printed condition is a first-year peak-temperature survival test, not a sufficient criterion for literal year-round ice in this model. For example, with $aH=1/2$ and very small positive $r$, it gives $h_m/H\simeq 2(\sqrt{2}-1)\simeq 0.828$, while first-year melt is approximately $0.5H$ but full warm-season melt is approximately $H$. Ice is present at the year's endpoint and disappears later that same summer. No reinterpretation of the printed inequality as guaranteed perennial cover is justified without an additional seasonal assumption.