Previous studies have shown that the large-scale zonal mean EP flux is parallel to the group velocity of Rossby waves. Thus, this section examines the relationship between the group velocity and the EP flux, demonstrating that the zonal, meridional and vertical components of the 3D EP flux in Eqs. (20)-(22) is parallel to the group velocity of inertia-gravity waves, representing the wave propagation.
The sinusoidal wave form for a' can be considered: a'=A0 eiφ, where φ=kx+ly+nz-ω t is the phase of the waves; ω is the local frequency; k, l and n are the zonal, meridional and vertical wave numbers; i is an imaginary number; and A0 is the amplitude of the disturbance. By substituting the sinusoidal wave form into the linear equations [see Eqs. (A1)-(A5) in Appendix A], the disturbance amplitudes of the inertia-gravity waves are obtained as \begin{eqnarray} \label{eq22} u_0&=&\frac{k\omega+ilf}{(\omega^2-f^2)}\frac{p_0}{\bar {\rho}} ,\ \ (22)\\ \label{eq23} v_0&=&\frac{l\omega-ikf}{(\omega^2-f^2)}\frac{p_0}{\bar {\rho}} ,\ \ (23)\\ \label{eq24} w_0&=&\frac{\omega(l^2+k^2)}{(f^2-\omega^2)n}\frac{p_0}{\bar {\rho}} =\frac{n\omega}{(\omega^2-N^2)}\frac{p_0}{\bar {\rho}} ,\ \ (24)\\ \label{eq25} \theta_0&=&\frac{inN^2\bar {\theta}}{(N^2-\omega^2)g}\frac{p_0}{\bar {\rho}} . \ \ (25)\end{eqnarray}
It is assumed that the disturbance amplitudes p0 can be separated into two parts, p0=p 0r+ip 0i (p 0r is the real part; p 0i is the imaginary part); thus, the disturbance amplitudes in Eqs. (23)-(26) can be written as \begin{eqnarray} \label{eq26} u_0&=&\frac{k\omega p_{\rm 0r}-lfp_{\rm 0i}}{(\omega^2-f^2)\rho_0}+i\frac{k\omega p_{\rm 0i}+lfp_{\rm 0r}}{(\omega^2-f^2)\rho_0} ,\ \ (26)\\ \label{eq27} v_0&=&\frac{kfp_{\rm 0i}-l\omega p_{\rm 0r}}{(\omega^2-f^2)\rho_0}+i\frac{l\omega p_{\rm 0i}-kfp_{\rm 0r}}{(\omega^2-f^2)\rho_0} ,\ \ (27)\\ \label{eq28} w_0&=&-\frac{(k^2+l^2)\omega p_{\rm 0r}}{n(\omega^2-f^2)\rho_0}-i\frac{(k^2+l^2)\omega p_{\rm 0i}}{n(\omega^2-f^2)\rho_0} ,\ \ (28)\\ \label{eq29} \theta_0&=&-\frac{nN^2\bar {\theta}p_{\rm 0i}}{(N^2-\omega^2)g\rho_0}+i\frac{nN^2\bar {\theta}p_{\rm 0r}}{(N^2-\omega^2)g\rho_0} . \ \ (29)\end{eqnarray}
For arbitrary waves, a'=A eiφ and b'=B eiφ, and the Reynolds averaging of their product is $$ \overline{{a}'{b}'}=\frac{1}{2}{\rm R}(A^\otimes B)=\frac{1}{4}(A^\otimes B+B^\otimes A), $$ where A and B are the amplitudes, "$\otimes$" denotes complex conjugation, and R denotes taking the real part. According to the formula above, the disturbance energy E can be written as \begin{eqnarray} \label{eq30} \overline{E}&=&\frac{1}{2}\left(\overline{{u}'^2}+\overline{{v}'^2}+\overline{{w}'^2}+ \frac{\overline{{\theta}'^2}}{N^2}\frac{g^2}{\bar {\theta}^2}\right)\nonumber=\frac{1}{4}{\rm R}\left(u_0^\otimes u_0+v_0^\otimes v_0+w_0^\otimes w_0+ \frac{\theta_0^\otimes\theta_0}{N^2}\frac{g^2}{\bar {\theta}^2}\right) . \ \ (30)\end{eqnarray}
By substituting the disturbance amplitude expressions [Eqs. (27)-(30)] and the dispersion relation ω2(k2+l2+n2)=N2(k2+l2)+f2n2 into Eq. (31), the disturbance energy is obtained: \begin{equation} \label{eq31} \bar {E}=\frac{\omega^2(k^2+l^2)(k^2+l^2+n^2)}{(f^2-\omega^2)^2n^2}\frac{\overline{{p}'^2}}{\bar {\rho}^2} . \ \ (31)\end{equation}
The wave activity density W, which is defined as the ratio of the perturbation energy to the intrinsic phase velocity (Kinoshita and Sato, 2013a), can be written as \begin{eqnarray} \label{eq32} \overline{W_1}&=&\frac{\overline{E}}{C_x}=\frac{\omega k(k^2+l^2)(k^2+l^2+n^2)}{(f^2-\omega^2)^2n^2} \frac{\overline{{p}'^2}}{\bar {\rho}^2}=\frac{k(k^2+l^2+n^2)}{\omega(k^2+l^2)}w_0^2 ,\ \ (32)\nonumber\\\\ \label{eq33} \overline{W_2}&=&\frac{\overline{E}}{C_y}=\frac{\omega l(l^2+k^2)(n^2+l^2+k^2)}{(\omega^2-f^2)^2n^2} \frac{\overline{{p}'^2}}{\bar {\rho}^2}=\frac{l(k^2+l^2+n^2)}{\omega(l^2+k^2)}w_0^2 .\ \ (33)\nonumber\\ \end{eqnarray} where $\overline{W_1}$ and $\overline{W_2}$ are two directions of horizontal wave activity density; Cx=ω/k and Cy=ω/l are the zonal and meridional components of the phase velocity, respectively. Previous work has proven that the wave activity density H satisfies the 3D wave activity conservation equation on the f-plane (for the detailed derivation, see Appendix A). The wave activity conservation equation is \begin{equation} \label{eq34} \frac{\partial H}{\partial t}+\frac{\partial}{\partial x}(C_{{g},x}H)+ \frac{\partial}{\partial y}(C_{{g},y} H)+\frac{\partial}{\partial z}(C_{{g},z}H)=0 , \ \ (34)\end{equation} where H=(k2+l2+n2)ω02 is the wave activity density and Cg=(Cg,x,Cg,y,Cg,z) are the group velocities of the 3D inertia-gravity waves [see Eqs. (A22)-(A24) in Appendix A].
According to the above expression, the wave activity density of Eqs. (33) and (34) can be written as $$ \overline{W_1}=\frac{kH}{\omega(k^2+l^2)},\quad \overline{W_2}=\frac{lH}{\omega(k^2+l^2)}. $$
Then, the tendency equation of the wave activity density is obtained as \begin{eqnarray} \label{eq35} &&\frac{\partial\overline{W_1}}{\partial t}+\frac{\partial}{\partial x}(C_{{g},x}\overline{W_1}) +\frac{\partial}{\partial y}(C_{{g},y}\overline{W_1})+\frac{\partial}{\partial z}(C_{{g},z}\overline{W_1})\nonumber\\ &&=\frac{k}{\omega(k^2+l^2)}\left[\frac{\partial}{\partial t}H+\nabla\cdot({C}_{g}H)\right]+\nonumber\\ &&\quad H\left[\frac{\partial}{\partial t}\left(\frac{k}{\omega(k^2+l^2)}\right)+{C}_{g}\cdot \nabla\left(\frac{k}{\omega(k^2+l^2)}\right)\right] ,\ \ (35)\\ \label{eq36} &&\frac{\partial\overline{W_2}}{\partial t}+\frac{\partial}{\partial x}(C_{{g},x}\overline{W_2}) +\frac{\partial}{\partial y}(C_{{g},y}\overline{W_2})+\frac{\partial}{\partial z}(C_{{g},z}\overline{W_2})\nonumber\\ &&=\frac{l}{\omega(k^2+l^2)}\left[\frac{\partial}{\partial t}H+\nabla\cdot({C}_{g}H)\right]+\nonumber\\ &&\quad H\left[\frac{\partial}{\partial t}\left(\frac{l}{\omega(k^2+l^2)}\right)+{C}_{g}\cdot \nabla\left(\frac{l}{\omega(k^2+l^2)}\right)\right] . \ \ (36)\end{eqnarray}
By using the local frequency, the horizontal wave number and wave activity conservation equation [Eq. (35)] of inertia-gravity waves, Eqs. (36) and (37), can be written as \begin{eqnarray} \label{eq37} \frac{\partial\overline{W_1}}{\partial t}+\frac{\partial}{\partial x}(C_{{g},x}\overline{W_1})+ \frac{\partial}{\partial y}(C_{{g},y}\overline{W_1})+\frac{\partial}{\partial z}(C_{{g},z}\overline{W_1})=0 ,\ \ (37)\\ \label{eq38} \frac{\partial\overline{W_2}}{\partial t}+\frac{\partial}{\partial x}(C_{{g},x}\overline{W_2})+ \frac{\partial}{\partial y}(C_{{g},y}\overline{W_2})+\frac{\partial}{\partial z}(C_{{g},z}\overline{W_2})=0 . \ \ (38)\end{eqnarray}
The above equations demonstrate that the wave activity densities $\overline{W_1}$ and $\overline{W_2}$ satisfy the wave activity conservation equation, which means that $\overline{W_1}$ and $\overline{W_2}$ are conserved during wave propagation.
By using the Reynolds averaging and the disturbance amplitudes, Eqs. (27)-(30), the 3D EP fluxes, Eqs. (20)-(22), can be rewritten as \begin{eqnarray} \label{eq39} F_{1,1}&=&\frac{k^2}{(\omega^2-f^2)}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\quad F_{1,2}=\frac{lk}{(\omega^2-f^2)}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\quad F_{1,3}=\frac{k(l^2+k^2)}{n(f^2-\omega^2)}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\ \ (39)\\ \label{eq40} F_{2,1}&=&\frac{lk}{(\omega^2-f^2)}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\quad F_{2,2}=\frac{l^2}{(\omega^2-f^2)}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\quad F_{2,3}=\frac{l(l^2+k^2)}{n(f^2-\omega^2)}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\ \ (40)\\ \label{eq41} F_{3,1}&=&-\frac{\omega^2k(l^2+k^2)}{n(f^2-\omega^2)^2}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\quad F_{3,2}=-\frac{\omega^2l(l^2+k^2)}{n(f^2-\omega^2)^2}\frac{\overline{{p}'^2}}{\bar {\rho}^2} ,\quad F_{3,3}=\frac{\omega^2(l^2+k^2)^2}{n^2(f^2-\omega^2)^2}\frac{\overline{{p}'^2}}{\bar {\rho}^2} . \ \ (41)\end{eqnarray}
By using the group velocities and the wave activity densities [Eqs. (33) and (34)], it can be proven that \begin{eqnarray} \label{eq42} C_{{g},x}\overline{W_1}&=&F_{1,1} ,\quad C_{{g},y}\overline{W_1}=F_{1,2} ,\quad C_{{g},z}\overline{W_1}=F_{1,3} ,\ \ (42)\\ \label{eq43} C_{{g},x}\overline{W_2}&=&F_{2,1} ,\quad C_{{g},y}\overline{W_2}=F_{2,2} ,\quad C_{{g},z}\overline{W_1}=F_{2,3} ,\ \ (43)\\ \label{eq44} C_{{g},x}\Re&=&F_{3,1} ,\quad C_{{g},y}\Re=F_{3,2} ,\quad C_{{g},z}\Re=F_{3,3} ,\ \ (44)\qquad \end{eqnarray} where $$\Re=\dfrac{\omega^3(l^2+k^2)^2(l^2+k^2+n^2)}{n^3(f^2-\omega^2)^3}\left(\dfrac{\overline{{p}'^2}}{\bar {\rho}^2}\right)$$ is another form of the wave activity density. This wave activity is essentially disturbance energy, which represents the combined effect of thermodynamic disturbance and dynamic disturbance.
In summary, by defining the modified wave activity density as Eqs. (33) and (34) and the 3D EP flux as Eqs. (40)-(42), a relation for the 3D propagation of inertia-gravity waves is derived; that is, the 3D EP flux is equal to the product of the group velocity and the modified wave activity density in all directions: $F_1=C_g\overline{W_1}$, $F_2=C_g\overline{W_2}$ and $F_3=C_g\Re$. This is useful for examining inertia-gravity wave energy propagation.
In addition, it is verified that the obtained 3D residual mean circulation is equal to the sum of the Eulerian time-mean flow and the Stokes drift, which satisfies the relationship between the Stokes drift and the residual mean circulation (the detailed derivation process can be seen in Appendix B). The TEM equations derived above under non-hydrostatic equilibrium are useful in examining how the generation and dissipation of atmospheric waves drives the mean circulation, and in exploring the energy propagation of inertia-gravity waves in mesoscale systems.