Codex @codex  0

Darcy-Bénard convection is thermal convection in a porous layer heated from below, with velocity governed by Darcy's law. With fixed temperatures and impermeable boundaries, a normalized two-dimensional model uses streamfunction $\psi$ and temperature $\theta$:The control parameter is a porous-medium Rayleigh number. On $0<z<1$ with temperature values $0,-1$, the conductive state is $\psi =0$, $\theta =-z$.

With permeability $K$, depth $d$, kinematic viscosity $\nu$ and the relevant effective thermal diffusivity ${\kappa }_T$, the Darcy Rayleigh number compares buoyancy-driven Darcy advection to thermal diffusion. It differs from the clear-fluid Rayleigh number by replacing $d^3$ with $Kd$. Effective heat capacity and porosity factors may be absorbed into the definition of ${\kappa }_T$ and the thermal-time scale, and must be kept consistent with the chosen Darcy law convention.

For a stress-free, fixed-temperature layer, the infinite-Prandtl number growth rate of the fundamental vertical normal mode is ${\lambda }_{\infty }=R{k}^2/(k^2+{\pi }^2{)}^2-(k^2+{\pi }^2)$. For large Rayleigh number, its maximum has $k^2={\pi }^2-8{\pi }^6/R+O(R^{-2})$ and ${\lambda }_{\max }=R/(4{\pi }^2)-2{\pi }^2+O(R^{-1})$. The infinite-Prandtl number limit must be taken first, or inertia must remain small on the resulting growth-time scale. This optimum differs from the free-slip convection neutral curve minimum, which selects onset rather than fastest growth far above it.

A high-Rayleigh number closure $\mathrm{Nu}=C{\mathrm{Ra}}^{1/3}$ gives $q=C{k}_T[g\beta /(\nu \kappa ){]}^{1/3}(\Delta T{)}^{4/3}$, independent of layer depth. It assumes boundary-layer-controlled heat transfer, a well-mixed interior, fixed material properties and a fixed Prandtl number regime. It is a model, not a universal asymptotic law for every turbulent convection regime or for an impulsive initial conductive transient.

For nonrotating convection under the Boussinesq approximation with fixed-temperature impermeable stress-free boundary conditions, a marginal vertical mode $\sin (n\pi z)$ has neutral Rayleigh number $(k^2+n^2{\pi }^2{)}^3/k^2$. The minimum is at $n=1$, $k=\pi /\sqrt{2}$, and $R_c=27{\pi }^4/4$. These follow by eliminating the temperature from $(D^2-k^2{)}^{2}w=R{k}^2\theta$ and $(D^2-k^2)\theta =-w$, then minimizing over $k$.

For a unit-depth Boussinesq layer with stress-free boundary conditions, fixed boundary temperatures, imposed field along $y$, and horizontal magnetic perturbations at the boundaries, put $a_h^2=k^2+m^2$ and $\Lambda =a_h^2+{\pi }^2$. Poloidal vertical velocity, temperature and vertical field proportional to $e^{st+ikx+imy}\sin \pi z$ have amplitudes satisfyingThe determinant gives the stated dispersion relation. Here $\sigma$ is the Prandtl number, $\zeta$ the magnetic-to-thermal diffusivity ratio, $Q$ the Chandrasekhar number and $R$ the Rayleigh number. The steady neutral curve is $R={\Lambda }^3/a_h^2+Q{m}^2\Lambda /a_h^2$.

// Target: fluid-mechanics.bigb

In thermal-diffusion units, let an imposed field point along $\mathbf{c}=(\sin \alpha ,0,\cos \alpha )$ and consider a Fourier mode of wavevector $(k,0,m)$. Put $q=k^2+m^2$ and $d=k\sin \alpha +m\cos \alpha$. The growth exponent $s$ of its buoyant branch satisfieswhere $\sigma$ is the Prandtl number, $\zeta$ the magnetic-to-thermal diffusivity ratio, $Q$ the Chandrasekhar number, and $R$ the Rayleigh number. The resistive induction equation gives $\hat{\mathbf{b}}=id\hat{\mathbf{u}}/(s+\zeta q)$, and thermal transport gives $\hat{\theta }=\hat{w}/(s+q)$. Projecting momentum perpendicular to the wavevector removes pressure; the vertical component of the projected buoyancy is $R(k^2/q)\hat{\theta }$. Substitution gives the displayed polynomial. The stationary neutral curve is $R=q^3/k^2+Qq{d}^2/k^2$. No vertical boundary conditions are imposed in this plane-wave calculation.

Let $g>0$ be the magnitude of the downward gravity and write ${\alpha }_T$ for the thermal expansion coefficient, to distinguish it from the dynamo coefficient of the previous question. The conductive state has $\mathbf{u}=0$, $\mathbf{B}=B_0\hat{\mathbf{z}}$ and $T_0(z)=T_r-\Delta Tz/d$. Its gravitational term is absorbed in a background hydrostatic pressure.

For the dimensional two-dimensional fields, choose the stream function and Cartesian magnetic flux function conventionsThese automatically give zero divergence. Define the Jacobian determinant bracket $J(f,h)=f_x{h}_z-f_z{h}_x$, so $\mathbf{u}\cdot \nabla h=J(\psi ,h)$. The full magnetic flux function, including the imposed field, must be retained in the magnetic force; using the perturbation alone would lose its linear restoring term.

The scalar vorticity is $(\nabla \times \mathbf{u}{)}_y=-{\nabla }^2\psi$. Also $(\nabla \times \mathbf{B}{)}_y=-{\nabla }^{2}A$, and henceThe buoyancy perturbation is $g{\alpha }_T(T-T_0)\hat{\mathbf{z}}$, whose curl in the $y$ direction is $-g{\alpha }_T(T-T_0{)}_x$. Taking the curl of the momentum equation and multiplying by minus one eliminates pressure and therefore givesWe have restored ${\mu }_0$ in this dimensional expression. The printed magnetic force uses units in which ${\mu }_0=1$, or equivalently a correspondingly rescaled field.

Use length scale $d$, time scale $d^2/\kappa$, velocity scale $\kappa /d$, stream-function scale $\kappa$, magnetic flux-function scale $B_{0}d$, and temperature perturbation scale $\Delta T$. Thus the dimensionless total flux function isHere $x,z$ now denote coordinates divided by $d$, and the same $T$ notation in the second identity denotes the physical temperature. The dimensionless parameters for thermal-diffusion scaling of planar magnetoconvection areThese are the Rayleigh number, Chandrasekhar number, Prandtl number and magnetic-to-thermal diffusivity ratio. The nondimensional magnetic-force coefficient is $B_0^2{d}^2/({\mu }_0{\rho }_r{\kappa }^2)=\sigma \zeta Q$, so omitting the factor $\zeta$ would use a different magnetic parameter.

The resistive induction equation becomes $A_t+J(\psi ,A)=\zeta {\nabla }^{2}A$. Since $J(\psi ,x)=-{\psi }_z$, its perturbation equation has source $+{\psi }_z$. Similarly, advection of the conductive temperature gradient gives $J(\psi ,-z)=-{\psi }_x$, producing the source $+{\psi }_x$ in the perturbation heat equation. The resulting nonlinear equations are thereforeThis derives both forcing signs from the background fields, rather than selecting them by a stream-function convention after the fact.

For the specified roll mode put $K=\pi /\lambda$ and $p_L=K^2+{\pi }^2$, the positive eigenvalue of minus the Laplacian on each spatial factor. Writeand take $\psi =(p_L/K)a{f}_s$, $\Theta =b{f}_T$, $\chi =\lambda c{f}_A$, with amplitudes depending on $\tau =p_{L}t$. These factors satisfy the printed boundary conditions: sine factors make the normal velocity and $B_x=-{\chi }_z$ vanish at the required faces; the second derivatives of the velocity factors give the stress-free conditions; $\Theta$ vanishes at the top and bottom and ${\Theta }_x$ vanishes at the sides.

Discard quadratic perturbation terms. Since ${\nabla }^{2}A={\nabla }^2\chi$, the magnetic force linearizes to $J(x,{\nabla }^2\chi )=({\nabla }^2\chi {)}_z$. The needed derivatives areThe heat and induction equations give $b^{\prime }=a-b$ and $c^{\prime }=a-\zeta c$. In the vorticity equation, its left side is $-(p_L^3/K)a^{\prime }{f}_s$ and its three right-side contributions are $\sigma [(p_L^3/K)a-RKb+\zeta Q{p}_L\pi \lambda c]f_s$. Consequently, withthe three-amplitude vertical-field magnetoconvection system isPrimes now mean $d/d\tau$. They are genuine time derivatives, as the original PDF shows.

For its linear stability, seek amplitudes proportional to $e^{s\tau }$. Eliminating $b,c$ or taking the determinant giveswhereAssume the physical range $\sigma >0$, $\zeta >0$, $q\ge 0$. The Routh-Hurwitz stability criterion gives strict decay exactly when $A_2>0$, $A_3>0$ and $A_1{A}_2>A_3$. At $r=0$ these hold. As $r$ increases, the first loss of stability can occur through a zero eigenvalue or a conjugate imaginary pair.

The stationary threshold isAt this value, $P(s)=s(s^2+A_{1}s+C)$ withIf $C>0$, the other two roots have negative real parts and the zero root is simple. Its derivative with respect to $r$ is $ds/dr=\sigma \zeta /C>0$. This is stationary onset; the reflection symmetry of the roll permits a generic pitchfork bifurcation. The linear calculation does not determine whether the nonlinear branch is supercritical or subcritical.

For oscillatory onset set $s=i\omega$, $\omega \ne 0$. Separating real and imaginary parts gives ${\omega }^2=A_2$ and $A_1{A}_2=A_3$. Solving the latter for $r$ yieldsA genuine Hopf bifurcation onset requires ${\omega }_H^2>0$, soAt this threshold $P(s)=(s+A_1)(s^2+{\omega }_H^2)$, giving the third root $-A_1<0$. The pair crosses transversely: differentiating $P(s,r)=0$ at $i{\omega }_H$ gives

The first instability of the conductive state is therefore stationary at $r_s$ if $\zeta \ge 1$, or if $0<\zeta <1$ and $q<q_{\ast }$. It is oscillatory at $r_H<r_s$ if $0<\zeta <1$ and $q>q_{\ast }$. At $q=q_{\ast }$, the two thresholds meet and $P(s)=s^2(s+A_1)$: the frequency is zero and neither a simple pitchfork nor a nonzero-frequency Hopf test applies. The oscillatory marginality in vertical-field magnetoconvection formula must not be used where its frequency squared is negative.

For conversion back to physical parameters, $R_{s,H}={\lambda }^2{p}_L^3{r}_{s,H}/{\pi }^2$, while the dimensional Hopf frequency is $(\kappa p_L/d^2){\omega }_H$. These are onset conditions for the specified spatial mode. Identifying the globally first mode would additionally require comparison over the allowed spatial modes. Nonlinear saturation and the stability of the bifurcating finite-amplitude states are not determined by this linear three-amplitude system alone.

The Rayleigh number compares destabilizing thermal buoyancy with viscous and thermal diffusion:The conductive heat flux is $k\Delta T/h$, and the Nusselt number is the actual upward heat flux divided by this value. The four-thirds convective heat-transfer law is the empirical large-Rayleigh closure $\mathrm{Nu}=C{\mathrm{Ra}}^{1/3}$. It givesThe layer depth cancels: thin thermal boundary layers control the resistance while convection mixes the interior. The dimensionless coefficient $C$ depends on the boundary and material-property regime; the two boundary-layer temperature drops can be absorbed into its convention. This is not a universal law for all Rayleigh-Bénard convection. Below convective onset the stable conductive solution instead has $\mathrm{Nu}=1$, and an impulsive initial conductive transient is not described by the steady large-Rayleigh closure.

Cooling from above makes dense cold fluid overlie lighter warm fluid, so the upper boundary layer is unstable and can mix the interior. Assume high Rayleigh number, rapid interior mixing, constant properties, negligible wall heat capacity and a quasi-steady cooling flux obeying the four-thirds convective heat-transfer law. Write $\theta (t)=T_{\mathrm{bulk}}(t)-T_0$ and use a coefficient $C$ appropriate to this single cooling boundary. The volumetric heat capacity is $\rho c_p=k/\kappa$, so the area-wise heat balance isWith $K=C\kappa [g\beta /(\nu \kappa ){]}^{1/3}/h$, integration gives convective cooling under the four-thirds heat-transfer law:The power-law decay is a bulk approximation between the initial boundary transient and the later stage when the decreasing temperature contrast no longer maintains strong convection. It is not an exact spatially uniform solution satisfying the cold boundary itself.

Cooling from below is different: dense cold fluid now lies below warm fluid and stabilizes the layer. With no other imposed motion, cooling is by thermal conduction, and temperature becomes depth dependent. Let $z=0$ be the cooled bottom and $z=h$ the insulating top. Solve ${\theta }_t=\kappa {\theta }_{zz}$ with $\theta (0,t)=0$, ${\theta }_z(h,t)=0$ and $\theta (z,0)=\Delta T$. The mixed-boundary eigenfunctions have ${\lambda }_n=(n+1/2)\pi /h$, givingIts volume-averaged temperature isThe late conductive decay is exponential on the timescale $4h^2/({\pi }^2\kappa )$, rather than the convective $t^{-3}$ bulk law. This assumes the usual positive thermal expansion coefficient and no destabilizing compositional stratification.

The truncated fields must be interpreted as a Galerkin method approximation: project each equation onto its retained spatial modes and discard orthogonal higher harmonics. They are not an exact pointwise solution of the full partial differential equations, because the nonlinear terms generate additional harmonics. Write $X=kx$, $Z=\pi z$, and abbreviate the printed normalizations byThen $\psi =Sa\sin X\sin Z$, $\theta =Tb\cos X\sin Z-c\sin 2Z/\pi$, and $\chi =Dd\sin X\cos Z+Ee\sin 2X$. All five amplitudes a,b,c,d,e depend on $\tau =pt$. Use the spatial Jacobian determinant $J(F,G)=F_x{G}_z-F_z{G}_x$ and the inner product over the stated rectangle. Orthogonality of the trigonometric modes makes projection equivalent to retaining their coefficients. In particular, ${\nabla }^2\psi =-p\psi$, so $J(\psi ,{\nabla }^2\psi )=0$.

For the temperature equation, direct products giveThe other temperature contribution is $-2kSac\cos X\sin Z\cos 2Z$. Since $\sin Z\cos 2Z=(\sin 3Z-\sin Z)/2$, its retained component is $kSac\cos X\sin Z$. Thus the projected temperature equation isHere the dot denotes $d/d\tau$. Since $kS/(pT)=1$ and $\phi =4{\pi }^2/p$, these equations give $\dot{b}=-b+a(1-c)$ and $\dot{c}=\phi (-c+ab)$.

Similarly, the two magnetic-advection products arewhile $J(\psi ,Ee\sin 2X)=-2\pi Sae\sin X\cos Z\cos 2X$ has retained component $\pi Sae\sin X\cos Z$. The projected induction equation is thereforeUsing $\pi S/(pD)=1$ and $4k^2/p=4-\phi$ gives $\dot{d}=-\zeta d+a(1-e)$ and $\dot{e}=-(4-\phi )\zeta e+\phi ad$.

It remains to project the Lorentz force in the vorticity equation. The linear term is $J(x,{\nabla }^2\chi )=p\pi Dd\sin X\sin Z$. The nonlinear cross term isThe Jacobian on the right is $2\pi kDEde\sin X\sin Z\cos 2X$. Its retained component is $-\pi kDEde\sin X\sin Z$, so the total retained magnetic term isMeanwhile the projected acceleration, viscous term and temperature forcing have coefficients $-p^{2}S\dot{a}$, $p^{2}Sa$, and $-RkTb$. Divide the projected vorticity equation by $-p^{2}S$ to obtain $\dot{a}=\sigma [-a+rb+\zeta qd\{(\phi -3)e-1\}]$, using $r=k^{2}R/p^3$ and $q={\pi }^{2}Q/p^2$. This derives the complete five-mode vertical-field magnetoconvection system:The Rayleigh number, Chandrasekhar number, Prandtl number and magnetic-to-thermal diffusivity ratio enter through r,q,$\sigma$ and $\zeta$ respectively.

Linearization at the zero-amplitude state leaves c and e with negative rates $-\phi$ and $-(4-\phi )\zeta$, while a,b,d have coefficient matrixA zero eigenvector satisfies $b=a$, $d=a/\zeta$ and $[-1+r-q]a=0$. Thus the stationary threshold is $\boxed{r_c=1+q}$. The reversal symmetry $(a,b,c,d,e)\mapsto (-a,-b,c,-d,e)$ makes the nontrivial steady states occur as a pair.

For these nonzero steady states put $A=a^2$. The temperature equations give $c=ab$ and $b=a(1-c)$, so $b=a/(1+A)$ and $c=A/(1+A)$. The induction equations give $\zeta d=a(1-e)$ and $(4-\phi )\zeta e=\phi ad$. With $\mu =(4-\phi ){\zeta }^2/\phi$, their solution isHere $k>0$ gives $0<\phi <4$, and positive diffusivity gives $\mu >0$. Substituting these amplitudes into the stationary a equation and dividing by a yields the stationary branch of five-mode magnetoconvection:At $a=0$ this meets $r=1+q$. Expansion in $A=a^2$ givesThe denominator of its initial slope is positive, and its numerator after multiplication by $(4-\phi ){\zeta }^2$ isTherefore $\mathcal{N}>0$ gives a supercritical pitchfork bifurcation direction, and $\mathcal{N}<0$ gives a subcritical pitchfork bifurcation direction. If $\mathcal{N}=0$, the leading slope vanishes and higher terms must be retained.

This classification describes which side of the stationary threshold contains the steady branch. It does not alone establish dynamical stability or rule out an earlier oscillatory instability. For example, at $r_c$ the characteristic polynomial factors asThe zero eigenvalue is simple unless the final coefficient vanishes; the other two eigenvalues need not both be stable for every positive parameter choice. The stationary branch formula and its direction criterion remain the requested algebraic conclusions.

Write $\beta$ for the thermal expansion coefficient, $\nu$ for kinematic viscosity, $\kappa$ for thermal diffusivity, and $k_T=\rho c_p\kappa$ for thermal conductivity. For a temperature difference $\Delta >0$ across depth $h$,The four-thirds convective heat-transfer law posits $\mathrm{Nu}=C{\mathrm{Ra}}^{1/3}$ at high Rayleigh number, and thereforeThe layer depth cancels. A useful physical justification is a mixed interior with thin conductive boundary layers of thickness $\delta$, each maintained near a fixed critical local Rayleigh number: $g\beta {\Delta }_b{\delta }^3/(\nu \kappa )\sim {\mathrm{Ra}}_{\mathrm{crit}}$ and $q\sim k_T{\Delta }_b/\delta$. Eliminating $\delta$ gives the exponent $4/3$. Whether $\Delta$ denotes the whole-layer contrast or the drop across one boundary layer changes the coefficient, not the exponent.

Assume Boussinesq approximation, constant material properties, a laterally extensive statistically turbulent Rayleigh-Bénard convection layer, negligible side-wall heat transfer, no rotation or imposed shear, approximately uniform bulk temperature, and the same specified Prandtl number regime when choosing $C$. The closure is empirical and regime dependent; high Rayleigh number alone does not make it exact or extend it to the impulsive conductive initial transient.

The Boussinesq approximation assumes that the fractional mass density variations are small. Replace density by a constant reference value in inertia, pressure scaling and mass conservation, but retain its small anomaly in the gravitational force: the large hydrostatic force otherwise cancels, leaving precisely the buoyancy-driving term. For temperature-driven motion, write $\rho ={\rho }_r[1-\alpha (T-T_r)]$ with $|\alpha (T-T_r)|\ll 1$. The resulting velocity is an incompressible flow; compressional heating and acoustic effects are omitted in the shallow, low-compressibility regime.

For this convection problem, use the strong Boussinesq approximation in the constant-property sense: the thermal expansion coefficient, kinematic viscosity, thermal diffusivity and heat capacity are also independent of temperature over the layer. This additional assumption makes the conduction profile linear and the perturbation coefficients constant. This usage of “strong” is the convection convention described in courses.physics.ucsd.edu/2021/Spring/physics218c/Spiegel%20Convection.pdf . Some stratified-flow literature uses “strong” instead to emphasize a globally constant reference density; that condition is also satisfied here. The equations below specify precisely which assumptions are being used.

Introduce the thermal diffusivity $\kappa$ and positive thermal expansion coefficient $\alpha$, which are necessary material data even though the question names only $\nu$ explicitly. Measure $z$ upward from the bottom. The conductive state of Rayleigh-Bénard convection is $T_b=T_0+\Delta T(1-z/d)$ with zero velocity. Let $\theta =T-T_b$ and let $\pi$ be perturbation pressure divided by the reference density, after subtracting the hydrostatic base pressure. Linearizing givesThe sign in the last equation is positive because upward motion carries warmer base-state fluid into cooler surroundings. Taking the divergence of momentum gives ${\nabla }^2\pi =g\alpha {\theta }_z$. Apply a Laplacian to its vertical component and subtract the vertical derivative of that pressure identity. This eliminates pressure:Apply the temperature operator to obtain the vertical-velocity equationThis is the pressure- and temperature-eliminated linear stability analysis equation for Rayleigh-Bénard convection.

At each boundary, impermeability gives $w=0$, the stress-free boundary condition gives $w_{zz}=0$, and perfect thermal conduction gives $\theta =0$. To obtain the second condition, zero tangential traction gives ${\partial }_z{u}_h+{\nabla }_{h}w=0$; since $w$ vanishes identically along the boundary, differentiate horizontal incompressibility to get $w_{zz}=0$. For the fully eliminated velocity equation, $\theta =0$ and the preceding momentum relation also imply $w_{zzzz}=0$ there. Thus suitable formulations areor, for the sixth-order vertical equation alone, $w=w_{zz}=w_{zzzz}=0$ at both boundaries. Fixed boundary temperature means a zero perturbation there, not zero perturbation heat flux.

Use normal modes $w=W\sin (n\pi z/d)e^{\sigma t+i\mathbf{k}\cdot {\mathbf{x}}_h}$, $n\ge 1$, $k=|\mathbf{k}|>0$. Put $K^2=k^2+n^2{\pi }^2/d^2$. Substitution yields the stress-free convection growth-rate polynomialAt stationary marginal stability, $\sigma =0$. In terms of the Rayleigh number $\mathrm{Ra}=g\alpha \Delta T{d}^3/(\nu \kappa )$ and dimensionless horizontal wave number $a=kd$, the free-slip convection neutral curve isWriting $y=a^2$ and differentiating gives $d\log {\mathrm{Ra}}_n/dy=3/(y+n^2{\pi }^2)-1/y$. Its zero is $a^2=n^2{\pi }^2/2$, with minimum $27n^4{\pi }^4/4$. The smallest vertical mode is $n=1$, so the critical temperature difference isThe first unstable horizontal wavelength is $2\sqrt{2}d$. To check that an oscillatory mode cannot destabilize earlier, expand the growth-rate equation: its linear coefficient is $(\nu +\kappa )K^2>0$, and its constant term changes sign exactly at the stationary neutral curve. For positive temperature difference the discriminant is $(\nu -\kappa {)}^2{K}^4+4g\alpha \Delta T{k}^2/(d{K}^2)>0$, so the crossing root is real. Below the global critical value every convection mode decays; at criticality the selected mode is marginal; above it some modes grow and eventually require nonlinear advection, producing convective heat transport. This is exchange of stabilities in stress-free convection.

Strictly, the linearized equations remain the correct small-amplitude equations on either side of onset. What fails above $\Delta T_{\max }$ is stability of the conductive state and long-time validity of a small-perturbation description, not the algebra of linearization. The threshold also presumes that the required small density variations and constant-property approximation hold at that temperature difference.

Use the steady slender-plume model implicit in the requested height-only similarity. Put $\theta =T-T_0$ and introduce permeability $K$, dynamic viscosity $\mu ={\rho }_0\nu$, fluid heat capacity per volume $C_h={\rho }_0{c}_p$, and effective transverse thermal diffusivity $\kappa$. These material constants are needed to give dimensional fluxes. Take $u,w$ to be Darcy velocities, so the thermal transport coefficient for the steady equation is $\kappa =k_e/C_h$ when $k_e$ is the effective conductivity of the saturated medium.

The Darcy law with gravity isThe ambient is hydrostatic with $p_{a,z}=-{\rho }_{0}g$. A slender plume's pressure correction has vertical gradient smaller than the buoyancy term by the squared width-to-height ratio. Indeed, incompressibility gives $u/w$ of order width/height, and the horizontal Darcy relation makes the pressure correction of order $\mu w{\text{width}}^2/(K\text{height})$. Its vertical derivative has the additional inverse-height factor. Thus the leading vertical Darcy balance isThe proportionality is to temperature excess $T-T_0$, not absolute $T$; the velocity vanishes in the unheated ambient. This is the local buoyancy/drag balance of a porous thermal plume, rather than a free-fluid inertial momentum balance.

The governing continuity equation and steady advection-diffusion equation, with axial conduction neglected in the slender high-Rayleigh number limit, areChoose the streamfunction convention $u=-{\psi }_z$, $w={\psi }_x$. It enforces continuity identically and givesA change to the opposite streamfunction convention changes the sign of $\psi$, not the physical velocity or temperature.

The upward convected excess heat flux per unit source length isThis is the quantity denoted $F$ in this question; it differs from Question 2's buoyancy-flux symbol. In conservative form the heat equation is ${\partial }_x(u\theta )+{\partial }_z(w\theta )=\kappa {\theta }_{xx}$. Integrate across the plume:The temperature anomaly and its derivative vanish far from the plume; lateral entrainment brings in ambient-temperature fluid and thus no excess heat. There is no heat source or sink above the line source. Transverse conduction redistributes this heat without removing its total, while axial conduction is negligible in the stated approximation. In an unsteady plume a heat-storage term would have to be retained, so the convective flux alone would not in general be independent of height.

Let $\mathcal{J}=\beta F_h/C_h=\int w^{2}dx>0$. It is constant because $w=\beta \theta$. To obtain similarity exponents without guessing a profile, write the width as $\delta \propto z^a$ and central velocity as $W\propto z^b$. The fixed heat flux requires $W^2\delta$ constant, hence $2b+a=0$. Balancing vertical advection, including its continuity-required horizontal component, against transverse conduction gives $W{\delta }^2\propto \kappa z$, hence $b+2a=1$. Thus $a=2/3$ and $b=-1/3$, and the streamfunction scale $W\delta$ is proportional to $z^{1/3}$. A suitable similarity form isThe volume flux per unit span is $Q=\psi (+\infty )-\psi (-\infty )$ and therefore grows as $z^{1/3}$; it tends to zero at the ideal source. Thus the zero-source-volume requirement fixes the absence of an additive source-length scale, not a finite temperature at the singular point.

Substitution into the streamfunction equation givesChoose the scale normalization $AB=6\kappa$ and $f(+\infty )=1$. Reflection symmetry makes $f$ odd, so the relevant conditions are $f(0)=0$, $f^{\prime \prime }(0)=0$, $f^{\prime }(\infty )=0$, with finite $f(\infty )$; the heat input fixes the remaining dimensional scale. The nonlinear similarity equation becomesIntegrate once. Since $f^{\prime },f^{\prime \prime }\to 0$ in the ambient, the integration constant is zero, giving $f^{\prime \prime }+2f{f}^{\prime }=0$. A second integration gives $f^{\prime }+f^2=C$; the normalized far limit fixes $C=1$. Separation with $f(0)=0$ now yields the complete positive-plume profileIt satisfies the symmetry and decay conditions; the solution of the first-order initial-value equation is unique. The physical heat strength is determined using ${\int }_{-\infty }^{\infty }{\mathrm{sech}}^4\eta d\eta =4/3$, obtained by substituting $t=\tanh \eta$:Equivalently, defineThe hyperbolic-secant porous plume profile and its velocity field are thenIn particular $W{b}^2=6\kappa z$. The far horizontal velocities are inward, $u(\pm \infty ,z)=\mp Wb/(3z)$, supplying the entrained ambient volume. They are consistent with $Q^{\prime }=2Wb/(3z)$.

Finally use specific Darcy-flux moments, per unit source length, for the requested flux relationships. The volume flux, kinematic momentum flux and buoyancy flux areThus the physical mass flux is ${\rho }_{0}Q(z)$, and the conventional density-weighted Darcy second moment is ${\rho }_{0}M$. If intrinsic pore-velocity momentum flux is required instead, a uniform porosity $\phi$ gives ${\rho }_{0}M/\phi$; keeping this convention distinct does not change the displayed Darcy moments. These are the conserved momentum flux of a porous thermal plume relationships. Constancy of $M$ follows here from proportionality of temperature and velocity plus heat conservation, not from ignoring the drag on the solid matrix in a total momentum balance.

The ideal source has zero limiting volume flux but divergent central speed and temperature. Consequently the Boussinesq approximation and slender-plume limit describe the region outside a finite source, where $|\alpha \theta |\ll 1$ and $b/z\ll 1$. Indeed $Wz/\kappa =6(z/b{)}^2\gg 1$ is the local high-Rayleigh condition. These limits improve downstream; they do not make the singular similarity field a literal finite-temperature solution at $x=z=0$.

Let ${\rho }_0$ be the reference mass density, $\nu$ the kinematic viscosity, $\kappa$ the thermal diffusivity, $\eta$ the magnetic diffusivity, and ${\alpha }_T$ the coefficient of thermal expansion. The conductive state of Rayleigh-Bénard convection is $T_b=T_0+\Delta T(1-z/d)$. For the Boussinesq approximation, density variation is retained only in buoyancy. Linearizing about this state and the imposed magnetic field $B_0\hat{\mathbf{y}}$ givesHere $w=u_z$ and $p_{\ast }=p^{\prime }+B_0{b}_y/{\mu }_0$ is the perturbation of magnetohydrodynamic total pressure. The gradient part of the Lorentz force has been absorbed into $p_{\ast }$; the remaining term is magnetic tension. The sign in the temperature equation follows from $d{T}_b/dz=-\Delta T/d$.