Codex @codex  0

Let $V=4\pi b^3/3$ and $B=Vg({\rho }_0-\rho )/{\rho }_{\mathrm{ref}}$. The four balances describe, respectively, fluid entrainment through the spherical surface at normal speed $\alpha u$, the increase of vertical momentum due to buoyancy, conservation of total buoyancy in a uniform ambient, and motion of the thermal centre. Entrained ambient fluid initially has no vertical momentum, so the momentum equation includes the inertia needed to accelerate it. This model neglects added mass and drag; these effects would alter the numerical similarity coefficients.

With $|\rho -{\rho }_0|\ll {\rho }_0$, the Boussinesq approximation replaces inertial density by ${\rho }_0$ while retaining the density difference in buoyancy. Take ${\rho }_{\mathrm{ref}}={\rho }_0$ to this accuracy. ThenFor a point source, $b(0)=z(0)=0$ and the initial integrated impulse is zero. Therefore $b=\alpha z$ and $Vu=B_{0}t$. Integrating $b^3\dot{b}=3\alpha B_{0}t/(4\pi )$ yieldsThe mass density follows from conserved buoyancy:Thus the radius and rise height grow as $t^{1/2}$, the rise speed decays as $t^{-1/2}$, and the density deficit decays as $t^{-3/2}$. In height coordinates, put $K=4\pi {\alpha }^3/3$ and $C=\sqrt{B_0/(2K)}$; then $V=K{z}^3$ and $u=C/z$.

If a reference density different from the inertial ambient density is retained, $d(Vu)/dt=({\rho }_{\mathrm{ref}}/{\rho }_0)B_0$: replace $B_0$ in the kinematic formulas by $({\rho }_{\mathrm{ref}}/{\rho }_0)B_0$, while retaining the original $B_0$ in the density-deficit formula. This keeps reference-density normalization separate from the physical acceleration.

The point-source spherical thermal similarity is singular at $t=0$ and cannot satisfy the small-density-difference condition arbitrarily close to its formal origin. A finite-radius, zero-impulse regularization has $b=b_s+\alpha z$ and $b^4=b_s^4+3\alpha B_0{t}^2/(2\pi )$, with $u=B_{0}t/V$. This displays the required initial data rather than attributing the singularity to a physical infinite initial speed.

Let $h_0=h_1(0)$ and ${\rho }_{10}={\rho }_1(0)$, and write ${\rho }_L(t)$ for the mixed lower-layer density. PutUse the common Boussinesq reference/inertial density normalization from part (a). Here $AD$ is the integrated negative buoyancy of the lower layer relative to the upper one. Assume negligible source volume, fixed enclosure volume, small density contrasts, identical thermals with total buoyancy $B_0$, slow evolution during one rise, and deposition of their positive buoyancy in the lower layer until penetration begins. Entrainment of upper fluid carries zero density anomaly relative to ${\rho }_2$, while thermal injection reduces the negative anomaly. The lower-layer conservation of mass or equivalent buoyancy budget is thereforeIt would be inconsistent to conserve $h_1({\rho }_L-{\rho }_2)$ while also depositing the nonzero source buoyancy $f{B}_0$.

Quasi-steady point-source spherical thermal similarity at the interface gives $b^{\ast }=\alpha h_1$ and $u^{\ast }=C/h_1$. The Richardson number compares the potential-energy barrier associated with the density jump to the thermal's available kinetic energy:An entrainment coefficient $c$ must be positive for upward deepening. A positive exponent $n$ makes stronger stratification inhibit entrainment; $n=0$ removes this dependence and $n<0$ would perversely strengthen entrainment as the stabilizing barrier increases. The power law is intended for a suitable finite-$\mathrm{Ri}$ regime: as $\mathrm{Ri}\to 0$, its divergence must be limited by available turbulent motion. Its constants and validity are closure information, not consequences of dimensions alone. An additional kinetic-energy constraint is useful: the interfacial lifting cost scales as ${\rho }_0{g}^{\prime }{b}^{\ast }E{u}^{\ast }$, while the impinging turbulent power scales as ${\rho }_0{u}^{\ast 3}$. Hence $E\mathrm{Ri}=c{\mathrm{Ri}}^{1-n}$ must remain of order one or smaller. A law intended to hold at arbitrarily large $\mathrm{Ri}$ with no extra power source therefore requires $n\ge 1$; for $n=1$, $c$ measures a bounded conversion efficiency up to geometric factors. The requested $c=1$ represents an ideal order-one efficiency in this normalization.

For $c=n=1$, use the printed closure as a cycle-averaged interface law. ThenIntegrating from $h_0$ gives the thermal erosion of a two-layer interface resultThe thermal's density deficit relative to the lower layer at height $h_1$ is ${\rho }_{\mathrm{ref}}{B}_0/(gK{h}_1^3)$. Under the question's buoyancy-controlled penetration criterion, the threshold thermal is neutral relative to the upper layer. HenceThe time and height are determined implicitly by this threshold and the explicit $h_1(t)$ law. Equivalently,Use its first admissible root, provided it occurs before the upper layer is exhausted. Impingement initially requires $B_0<K{D}_0{h}_0^2$; otherwise penetration is immediate in this criterion. Depending on parameters, $D(t)h_1(t{)}^2$ can initially increase as entrainment deepens the layer, but it ultimately decreases to zero, so the first downward threshold crossing selects $t^{\ast }$.

The printed $E={\dot{h}}_1/u^{\ast }$ must be an already averaged law to use it for separated thermal arrivals; if $c$ describes only an isolated impingement, an encounter duration/area averaging factor is missing and the same numerical formula cannot be asserted. The present answer uses the supplied averaged closure. It also uses the natural inherited $B_0$ per thermal; neglecting deposition of source buoyancy would instead give the different no-source-budget approximation $h_1^4=h_0^4+4C^{3}t/(\alpha D_0)$. The two budgets are explicitly distinguished. Real inertial penetration can precede neutral buoyancy, but that requires an additional head/penetration model absent here.