= Solution
Let $d=\psi_1-\psi_2$ and $C=f_0^2/g'=H_1F_1=H_2F_2$. Multiply each <two-layer quasi-geostrophic potential vorticity> equation by $H_i\psi_i$ and sum. The time-derivative terms are
$$
\sum_iH_i\psi_iq_{it}
=\nabla_h\cdot\sum_iH_i\psi_i\nabla_h\psi_{it}
-\partial_t\left[\frac12\sum_iH_i|\nabla_h\psi_i|^2+\frac C2d^2\right].
$$
The interface terms combine to $-Cd\,d_t$, since the two layers share the same depth-weighted coupling $H_iF_i$. The planetary term $\beta y$ has no time derivative. For the nonlinear terms, $\nabla_h\cdot\mathbf u_i=0$ and $\mathbf u_i\cdot\nabla_h\psi_i=0$ imply
$$
H_i\psi_iJ(\psi_i,q_i)=\nabla_h\cdot(H_i\psi_iq_i\mathbf u_i).
$$
Consequently the local <two-layer quasi-geostrophic energy conservation> law is
$$
\boxed{\partial_tE+\nabla_h\cdot\mathbf F_E=0,\qquad
E=\frac12\left[H_1|\nabla_h\psi_1|^2+H_2|\nabla_h\psi_2|^2+C(\psi_1-\psi_2)^2\right],}
$$
$$
\boxed{\mathbf F_E=-\sum_{i=1}^2H_i\left(\psi_i\nabla_h\psi_{it}+\psi_iq_i\mathbf u_i\right).}
$$
The flux expression is one convenient form obtained directly from the requested multiplication; divergence-free modifications would represent the same local balance. $E$ is physical energy per horizontal area divided by the common reference density. Multiplying both $E$ and $\mathbf F_E$ by that density restores the dimensional physical-energy convention.
Back to article page