Solution (source code)

= Solution

Begin with the <Boussinesq approximation> primitive equations on a <beta plane>, decompose every field into a zonal mean and a disturbance, and average over longitude. The zonal momentum equation then contains the divergence of the eddy momentum flux $\overline{u'v'}$, while the mean density equation contains the divergence of the eddy density flux $\overline{\rho'v'}$. At small <Rossby number>, use <geostrophic balance>, <hydrostatic pressure>, <thermal-wind balance>, and the leading eddy equations to combine those fluxes.

Define
$$
A=\frac{\overline{\rho'v'}}{d\rho_s/dz}
$$
and introduce the <residual mean circulation>
$$
\boxed{
\overline v_a^*=\overline v_a-A_z,
\qquad
\overline w_a^*=\overline w_a+A_y}.
$$
The added eddy-induced velocity is nondivergent, so
$$
\overline v_{a,y}^*+\overline w_{a,z}^*=0.
$$
It absorbs the eddy density-flux divergence into advection by the transformed circulation, giving
$$
\overline\rho_t+overline w_a^*\frac{d\rho_s}{dz}=0.
$$

The mean zonal momentum equation becomes
$$
\overline u_t-f_0\overline v_a^*
=\nabla\mathbin\cdot\overline{\mathbf F},
$$
where the zonally averaged <Eliassen–Palm flux> in the meridional-vertical plane is
$$
\boxed{
\overline F^{(y)}=-\overline{u'v'},
\qquad
\overline F^{(z)}
=f_0\frac{\overline{\rho'v'}}{d\rho_s/dz}}.
$$
Thus the <transformed Eulerian mean> gathers the wave forcing into one flux divergence and makes density evolve under one residual circulation.