Solution (source code)

= Solution

Keep the half-section convention of part (a). Define the <mass density>-weighted <mass flux>, <buoyancy flux> and <momentum flux> by
$$
\boxed{Q=\rho bw,\qquad F=(\rho_0-\rho)gbw,\qquad M=\rho bw^2.}
$$
The full-section fluxes are $2Q,2F,2M$. These factors must be changed consistently if full fluxes are used; mixing a full-width flux with a half-width <entrainment> source would introduce a spurious factor two.

For positive $Q,M$, $w=M/Q$, $\rho b=Q^2/M$, and $(\rho_0-\rho)gb=FQ/M$. The mass equation in part (a), the momentum equation, and the <buoyancy> equation therefore become
$$
\boxed{
\begin{aligned}
\partial_t(Q^2/M)+Q_z&=\rho_0\alpha M/Q,\\
Q_t+M_z&=FQ/M,\\
\partial_t(FQ/M)+F_z&=-\frac{\rho_0}{\rho}N^2Q.
\end{aligned}}
$$
These are direct algebraic rewritings with the <mass density> factors retained. In particular, $FQ/M$ is the <buoyancy> inventory per unit height as well as the integrated <buoyancy> force. If desired, $\rho$ can be reconstructed from $F/(gQ)=\rho_0/\rho-1$.

At leading <Boussinesq> order use a constant reference <mass density> $R$ in the inertial fluxes, $Q=Rbw$, $M=Rbw^2$, and write $g'=g(\rho_0-\rho)/R$ so $F=Rg'bw$. Set $\rho_0/\rho=1$ in the <stratification> source and use $N^2\simeq-g\rho_0'/R$. The <unsteady top-hat line-plume balances> are then
$$
\boxed{\partial_t(Q^2/M)+Q_z=R\alpha M/Q,\qquad
Q_t+M_z=FQ/M,\qquad
\partial_t(FQ/M)+F_z=-N^2Q.}
$$
Equivalently the kinematic fluxes $q=Q/R=bw$, $m=M/R=bw^2$ and $f=F/R=g'bw$ obey the same equations with $R$ removed from the first source. This explicit distinction between <mass density>-weighted and kinematic fluxes fixes the normalization for the power-law solution.