Solution (source code)

= Solution

Take inward <fluid entrainment> speed to be positive. The non-Boussinesq form of the <Batchelor entrainment hypothesis> is
$$
\boxed{v_e=\alpha w\sqrt{\rho/\rho_0},}
$$
where $\alpha$ is the <entrainment coefficient>. Define physical, cross-section-integrated fluxes for the <top-hat plume model> by
$$
\boxed{Q=\pi b^2\rho w,\qquad M=\pi b^2\rho w^2,\qquad
F=\pi b^2w g(\rho_0-\rho).}
$$
Thus $Q$ is <mass flux>, $M$ is <momentum flux>, and $F$ is a density-weighted <buoyancy flux>. The kinematic buoyancy flux in the <Boussinesq approximation> is $F/\rho_0$. These definitions keep the factor $\pi$; conventions suppressing it simply rescale all three fluxes consistently.

For a homogeneous ambient, the sectional volume, mass and momentum balances are
$$
\partial_t(\pi b^2)+\partial_z(\pi b^2w)=2\pi b v_e,
$$
$$
\partial_t(\pi b^2\rho)+\partial_zQ=2\pi b\rho_0v_e,\qquad
\partial_tQ+\partial_zM=\pi b^2g(\rho_0-\rho).
$$
The momentum source is the net upward <buoyancy> force; entrained ambient fluid initially supplies no vertical momentum. Subtracting the mass balance from $\rho_0$ times the volume balance gives conservation of the density deficit:
$$
\partial_t[\pi b^2g(\rho_0-\rho)]+\partial_zF=0.
$$
Now $w=M/Q$, $\pi b^2\rho=Q^2/M$, and $\pi b^2g(\rho_0-\rho)=QF/M$. Substitution produces the <non-Boussinesq top-hat plume equations>
$$
\boxed{\partial_t\left(\frac{Q^2}{M}\right)+\partial_zQ=2\alpha\sqrt{\pi\rho_0M},\qquad
\partial_tQ+\partial_zM=\frac{QF}{M},\qquad
\partial_t\left(\frac{QF}{M}\right)+\partial_zF=0.}
$$
For reconstructing the fields without a <Boussinesq approximation>,
$$
\rho=\frac{\rho_0gQ}{gQ+F},\qquad b^2=\frac{Q^2}{\pi\rho M},\qquad
\frac{g(\rho_0-\rho)}\rho=\frac FQ.
$$
The last expression uses the plume density in the acceleration denominator; it agrees with the ambient-density definition of <reduced gravity> to Boussinesq order.