Past exam of the mathematics course of the University of Cambridge / 2019 / iii / Paper 331 (Hydrodynamic Stability)

Let $V=U-c$ and take $k>0$ without loss of generality. The Boussinesq approximation replaces density by ${\rho }_0$ in inertia and pressure acceleration, while retaining the density perturbation in buoyancy. Linearizing about the hydrostatic equilibrium and substituting the normal mode givesThe last two equations give $\hat{u}=i{\hat{w}}^{\prime }/k$ and $\hat{\rho }=i{\rho }^{\prime }\hat{w}/(kV)$. The horizontal momentum equation givesSubstituting in the vertical momentum equation yields the Taylor–Goldstein equation:The buoyancy frequency is real for stable density stratification. The impermeable walls impose $\hat{w}(\pm L)=0$.

For a possibly unstable normal mode, $c_i\ne 0$ ensures $V=U-c$ never vanishes, so choose one continuous branch of $V^a$. Substitute $\hat{w}=V^{a}q$ into the Taylor–Goldstein equation and collect terms:Multiplication by $V^{2a}$ puts its first two terms in divergence form:Multiply by the complex conjugate $\overline{q}$ and apply integration by parts. Since $q(\pm L)=0$, the boundary term vanishes and the power-transformed Taylor–Goldstein energy identity isThe weight is $V^{2a}$, not $|V{|}^{2a}$: preserving its complex phase is essential for the subsequent stability proofs.

Choose $a=1/2$ in the power-transformed Taylor–Goldstein energy identity and let $Q=N^2-(U^{\prime }{)}^2/4$. The result isBecause $\mathrm{Im}V=-c_i$ and $\mathrm{Im}(1/V)=c_i/|V{|}^2$, its imaginary part givesIf $4N^2>(U^{\prime }{)}^2$ everywhere, the integral is strictly positive for every nonzero eigenfunction, so a mode with $c_i\ne 0$ is impossible. Thus the flow is neutrally stable to these inviscid normal modes. This is the Miles–Howard theorem, expressed as a lower bound $\mathrm{Ri}>1/4$ on the gradient Richardson number. It excludes exponential modal growth; it does not by itself exclude transient growth.

Choose $a=1$. The power-transformed Taylor–Goldstein energy identity reduces toThe right side is real. For an unstable mode, $c_i>0$, so taking the imaginary part givesConsequently the phase velocity is a weighted average:The weights are nonnegative and not identically zero. For a nonconstant shear profile, equality at an extremum would force $q$ to vanish on an interval where $U$ differs from that extremum; uniqueness for the regular ordinary differential equation would then give the zero mode. This is the phase-speed bound for unstable stratified shear modes.

Let $\hat{\eta }$ be the displacement amplitude of an interface. Its kinematic boundary condition is $ik(U-c)\hat{\eta }=\hat{w}$, on either side, so the first of the jump conditions for stratified inviscid shear flow isFor zero surface tension, pressure is continuous at the displaced interface. Since the background hydrostatic pressure has derivative $-g\rho$, linearization gives $[\hat{p}-g\rho \hat{\eta }{]}_{-}^{+}=0$. Using $\hat{p}={\rho }_0[(U-c){\hat{w}}^{\prime }-U^{\prime }\hat{w}]/(ik)$ and the kinematic condition givesThis derivation also works when $U$ jumps: the displacement is shared, while the vertical velocity need not be. It avoids ambiguous products of distributions at a discontinuous velocity profile.

Write $V_0=\Delta U/2$, $z_{\pm }=\pm h/2$, $B=g\Delta \rho /(2{\rho }_0)$ and $d=B/(2k)$. The velocity is linear through all three layers, so $U(z_{\pm })=\pm V_0$ and $U^{\prime }$ has no jump. Each interface has density drop $\Delta \rho /2$. The jump conditions for stratified inviscid shear flow becomeBetween interfaces, the Taylor–Goldstein equation is ${\hat{w}}^{\prime \prime }-k^2\hat{w}=0$. A decaying solution can therefore be writtenIts derivative jumps are $-2k{b}_{\pm }$, while its values at the interfaces are $b_{\pm }+e^{-kh}{b}_{\mp }$. Hence the coefficients satisfyThe determinant must vanish. With $\tilde{c}=c/V_0$, $\alpha =kh/2$ and $J=g\Delta \rho h/({\rho }_0\Delta U^2)$, one has $d/V_0^2=J/(2\alpha )$. Expanding gives the dispersion relation for two density interfaces in uniform shear:

Put $y={\tilde{c}}^2$ and $s=J/(2\alpha )>0$. The dispersion relation for two density interfaces in uniform shear isIts discriminant is $16s+4s^2{e}^{-4\alpha }>0$, so both roots $y_{\pm }$ are real. Their sum is positive; therefore one is negative precisely when their product is negative. That condition isSolving for $J$ givesIn this interval $y_{-}<0$, so $\tilde{c}=i\sqrt{-y_{-}}$ is the growing normal mode. At either endpoint the corresponding root is zero and the mode is neutral.

For large $\alpha$, the coupling $e^{-kh}=e^{-2\alpha }$ is exponentially weak. If each interface were isolated, its interfacial gravity waves would have dimensional phase velocities $U(z_{\pm })\pm \sqrt{d}$, orThe upper-interface wave travelling against its positive local current has speed $1-\sqrt{J/(2\alpha )}$. The lower-interface wave travelling against its negative current has speed $-1+\sqrt{J/(2\alpha )}$. Their laboratory speeds match at zero when $J=2\alpha$. Weak interaction then produces counterpropagating wave instability rather than two independent travelling waves.

The band has the large-wavenumber form $J=2\alpha \pm 2\alpha e^{-2\alpha }+O(\alpha e^{-4\alpha })$. Exactly at resonance, $\mathrm{Im}\tilde{c}\sim e^{-2\alpha }/2$, soThus the instability window and growth rate both become exponentially small as the interfaces become weakly coupled.

For distinct eigenvalues, the formula holds for $n=0,1$. If the lower-left entry of $L^n$ is $b_n=({\lambda }_1^n-{\lambda }_2^n)/({\lambda }_1-{\lambda }_2)$, multiplication by $L$ gives $b_{n+1}={\lambda }_1{b}_n+{\lambda }_2^n=({\lambda }_1^{n+1}-{\lambda }_2^{n+1})/({\lambda }_1-{\lambda }_2)$. This proves by mathematical induction thatSumming the absolutely convergent matrix exponential series givesFor the repeated eigenvalue ${\lambda }_1={\lambda }_2=\lambda$, the continuous limit is $A(t)=e^{\lambda t}\left(\begin{array}{cc}1& 0\\ t& 1\end{array}\right)$. This is also the exponential of a two-dimensional Jordan block.

Termwise differentiation of the matrix exponential gives $A^{\prime }=LA$ and $A(0)=I$. Thus $x(t)=A(t)x(0)$ satisfies $x^{\prime }=Lx$ and the prescribed initial condition. Uniqueness for a linear differential equation makes it the solution.

For real ${\lambda }_1,{\lambda }_2$, computeThus $L$ is a non-normal matrix. For distinct eigenvalues, takeThese satisfy $L{v}_j={\lambda }_j{v}_j$ but have inner product $v_1^T{v}_2=1$. The eigenvectors are therefore nonorthogonal. When the eigenvalues coincide, only the direction $v_2$ remains as an eigenvector; $L$ is a defective Jordan block.

Differentiate the quadratic energy:For $x_1\ne 0$, put $r=x_2/x_1$. Immediate energy growth occurs precisely when ${\lambda }_2{r}^2+r+{\lambda }_1>0$. With both eigenvalues negative, such directions exist iff $1-4{\lambda }_1{\lambda }_2>0$. Writing $\delta =\sqrt{1-4{\lambda }_1{\lambda }_2}$, the growing initial conditions are the two opposite open wedgesIf $x_1=0$, energy decreases. No energy growth at any time, for any initial condition, is possible exactly whenIndeed, this makes the symmetric part of a matrix $H$ a negative semidefinite matrix, so $\dot{E}\le 0$ along every trajectory. If the product is smaller, the wedges already supply counterexamples at $t=0$. This is the instantaneous energy-growth criterion for a linear system.

At each fixed time, the maximum over nonzero initial conditions is the optimal energy amplification of a linear system:Write $a=e^{{\lambda }_{1}t}$, $d=e^{{\lambda }_{2}t}$ and $b=(a-d)/({\lambda }_1-{\lambda }_2)$. The symmetric matrix $A^{T}A$ has trace $\tau =a^2+b^2+d^2$ and determinant $a^2{d}^2$. Its largest eigenvalue isAn optimal initial condition is a right singular vector of $A(t)$ associated with its largest singular value.

If ${\lambda }_2>{\lambda }_1$, let $\alpha =1/({\lambda }_2-{\lambda }_1)$. Then $A(t)/e^{{\lambda }_{2}t}\to \left(\begin{array}{cc}0& 0\\ \alpha & 1\end{array}\right)$, so the precise long-time asymptotic statement isThe Cauchy-Schwarz inequality shows that the unit initial conditions achieving this leading factor areThis is an adjoint eigenvector, since $L^T(\alpha ,1{)}^T={\lambda }_2(\alpha ,1{)}^T$. It differs from the right eigenvector $v_2=(0,1{)}^T$: the initial condition maximizes the projection onto the slow mode, whereas the eventual state aligns with that right eigenvector. For ${\lambda }_2<0$, even the optimal energy ultimately decays to zero; the prefactor describes enhanced excitation of the slow mode.

The zero streamfunction gives zero velocity. The temperature ${\theta }_0=-z$ obeys both prescribed temperature values, has no horizontal gradient, and has zero Laplacian. Thus it satisfies both equations with ${\psi }_0=0$:This is the conductive state of Darcy-Bénard convection.

Linearizing about the conductive state givesThe perturbations vanish at $z=0,1$. For a normal mode with vertical dependence $\sin (n\pi z)$, put $q=k^2+n^2{\pi }^2$. The two modal equations are $-q\hat{\psi }=-ikR\hat{\theta }$ and $s\hat{\theta }+ik\hat{\psi }=-q\hat{\theta }$. Eliminating $\hat{\psi }$ givesNeutrality requires $R=(k^2+n^2{\pi }^2{)}^2/k^2$. The arithmetic-geometric mean inequality bounds this below by $4n^2{\pi }^2$, with equality at $k=n\pi$. The first unstable vertical mode is therefore $n=1$, givingThe wavelength statement applies at the first onset. At a general larger $R$, maximizing the growth rate within a fixed vertical mode gives $k^2=n\pi \sqrt{R}-n^2{\pi }^2$, so its preferred wavelength need not remain exactly two. This is the onset of convection in a horizontal Darcy layer.

The full perturbation equations arewhere the streamfunction advection bracket is $J(f,g)=f_z{g}_x-f_x{g}_z$. Put $R=R_c+{ϵ}^2$ and introduce the slow time $T={ϵ}^{2}t$, so ${\partial }_t={ϵ}^2{\partial }_T$ on the amplitude-dependent fields. ThenThe weakly nonlinear expansion balances the small linear growth against the cubic saturation on this slow time scale.

At first order, ${\nabla }^2{\psi }_1+R_c{\theta }_{1x}=0$ and ${\nabla }^2{\theta }_1-{\psi }_{1x}=0$. Since ${\nabla }^2{\psi }_1=-2{\pi }^2{\psi }_1$ and $R_c=4{\pi }^2$, the first equation fixes ${\theta }_{1x}={\psi }_1/2$. The homogeneous temperature boundary conditions and the second equation remove any horizontally uniform term, givingAt second order,The forcing is horizontally uniform. Taking no additional critical-mode component, a solution satisfying the homogeneous boundaries isThis mean-temperature correction in weakly nonlinear Darcy convection modifies the vertical temperature gradient and provides the feedback that saturates the convection.

Let ${\Psi }_3=A^3{\psi }_3$ and ${\Theta }_3=A^3{\theta }_3$ denote the full third-order coefficients. This notation keeps the detuning term, proportional to $A$, and the slow-time term, proportional to $A_T$, visible. DefineThe equations at order ${ϵ}^3$ areEquivalently, where $A\ne 0$, the question's coefficients satisfy ${\nabla }^2{\psi }_3+R_c{\theta }_{3x}=-A^{-2}{\theta }_{1x}$ and ${\nabla }^2{\theta }_3-{\psi }_{3x}=(A_T/A^3){\theta }_1+B$. Eliminate ${\Theta }_3$ from the full coefficients:Use the Fredholm solvability condition for a self-adjoint operator. With $⟨f,g⟩={\int }_0^1{\int }_0^{2}fgdxdz$, projection onto the critical eigenfunction ${\psi }_1$ makes the left side vanish. Therefore the Landau amplitude equation isThese integrals already answer the requested coefficient calculation. Evaluating them gives $⟨{\psi }_1,{\theta }_{1x}⟩=1/4$, $⟨{\psi }_1,{\partial }_x{\nabla }^2{\theta }_1⟩=-{\pi }^2/2$, and $⟨{\psi }_1,{\partial }_{x}B⟩={\pi }^2/32$. Hence, for the prescribed normalization,The cubic term opposes growth, so the nonzero steady amplitudes $A=\pm 2/\pi$ are stable equilibria within the chosen roll phase. The conductive state $A=0$ is unstable for this positive detuning. This is supercritical saturation of Darcy convection rolls.