HTML
--> --> -->The EP flux was first derived by (Eliassen and Palm, 1961). Then, Andrews and McIntyre (1976, 1978) proposed the generalized EP flux in the transformed Eulerian mean (TEM) equations by bringing in the residual mean flow. The generalized EP flux transfers the disturbance heat into the momentum equation to describe the dynamic and thermodynamic effects. The TEM equation is \begin{equation} \label{eq1} \frac{\partial\bar {u}}{\partial t}-f\bar {v}^\ast=\frac{1}{\rho_0}\nabla\cdot{F} , \ \ (1)\end{equation} where, \begin{equation} \label{eq2} {F}=-\rho_0(\overline{u'v'}){j}+\frac{f\rho_0}{N^2}(\overline{v'\phi'_z}){k} \ \ (2)\end{equation} is the 2-D EP flux; $\bar v^\ast$ is the residual mean flow; $\bar u$ is the zonal mean flow; u' and v' are disturbance wind speed; φz is vertical change of potential height; j and k represent the horizontal and vertical component; f is the Coriolis parameter; ρ0 is the basic density; N2 is the buoyancy frequency; and F is the EP flux. The residual mean flow and EP flux are powerful tools to diagnose the meridional circulation and wave-mean interaction in the atmosphere. The residual mean circulation is the sum of the Eulerian-mean flow and the Stokes drift in small-amplitude theory. The Stokes drift is the difference between the Eulerian-mean flow and the Lagrangian-mean flow under the small-amplitude assumption. The EP flux is expressed as the product of the wave activity density and the group velocity under the Wentzel-Kramers-Brillouin (WKB) approximation, representing the wave propagation (Edmon et al., 1980). When the divergence of the EP flux is positive, the zonal mean flow is accelerated. For linear, steady and conservative waves in a purely zonal basic flow, the divergence of the EP flux vanishes; thus, the waves neither drive nor accelerate the zonal mean flow. This is called the non-acceleration theorem (Eliassen and Palm, 1961). Subsequent studies (e.g., Zeng, 1982; Andrews, 1983; Wu and Chen, 1989; Pfeffer, 1992; Ding and Shen, 1998; Gao et al., 2004) have made many efforts to extend the EP flux to different situations. However, the two-dimensional EP flux is limited when analyzing the 3D wave propagation and wave activity; there are many studies generalizing the TEM equation to three dimensions. (Hoskins et al., 1983) deduced the three-dimensional (3D) EP flux based on the horizontal velocity correlation tensor, but their wave activity flux failed to parallel to the group velocity. (Plumb, 1986) obtained the 3D EP flux and the corresponding residual mean flow of Rossby waves under the quasi-geostrophic approximation using the potential vorticity conservation theorem. His wave activity flux is equal to the product of the wave activity density and the group velocity. However, the residual mean flow is unlike the TEM equation. (Takaya and Nakamura, 1997) derived phase-independent 3D wave-activity flux using the pseudomomentum wave energy conservation laws, which are applicable to quasi-stationary Rossby waves in zonally varying basic flow. (Miyahara, 2006) formulated the 3D EP flux and the residual mean flow applicable to gravity waves in a log-pressure coordinate system in a Boussinesq fluid. His 3D EP flux is equal to the product of the wave activity density and the group velocity, and the residual mean circulation is equal to the Stokes drift under the WKB approximation. (Noda, 2010) derived the 3D wave activity flux for a plane wave, but his formulas can only be used for a purely monochromatic wave. Recently, Kinoshita et al. (Kinoshita et al., 2010; Kinoshita and Sato, 2013a, Kinoshita and Sato, 2013b) formulated the 3D EP flux and corresponding residual mean circulation for both inertia-gravity waves and Rossby waves under the WKB approximation and proved that the 3D EP flux is equal to the product of the wave activity density and the group velocity, and the residual mean flow is the sum of the Eulerian-mean flow and the Stokes drift. But, his three dimensional wave activity flux is assumed under the hydrostatic approximation, and only applicable to large-scale systems.
However, inertia-gravity waves are the most ubiquitous in mesoscale systems, such as strong convection events or torrential rain events. In addition, the wave activity flux is an effective tool to analyze the influence of inertia-gravity waves on mesoscale systems. Thus, the purpose of this study is to formulate a 3D TEM equation and 3D EP flux that is applicable to non-hydrostatically balanced and non-geostrophically balanced systems, to examine the feedback to mesoscale systems from the inertia-gravity wave activity.
The paper is arranged as follows. In section 2, the 3D activity flux is derived from the time-mean horizontal momentum equation. Section 3 proves that the 3D EP flux is parallel to the group velocity and the residual mean flow satisfies the continuity equation. Then, section 4 takes a heavy precipitation event as an example to explore the application of the 3D EP flux. Finally, a summary and conclusions are provided in section 5.
Consider an adiabatic, inviscid, uncompressible and moist atmosphere. The governing equations in f-plane Cartesian coordinates are given by \begin{eqnarray} \label{eq3} &&\frac{du}{dt}-fv=-\frac{1}{\rho}\frac{\partial p}{\partial x} ,\ \ (3)\\ \label{eq4} &&\frac{dv}{dt}+fu=-\frac{1}{\rho }\frac{\partial p}{\partial y} ,\ \ (4)\\ \label{eq5} &&\frac{dw}{dt}=-\frac{1}{\rho }\frac{\partial p}{\partial z}-g ,\ \ (5)\\ \label{eq6} &&\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0 ,\ \ (6)\\ \label{eq7} &&\frac{d\theta}{dt}=0 , \ \ (7)\end{eqnarray} where x,y,z represent horizontal and vertical component of physical quantity; t is time; u, v and w are the zonal, meridional and vertical velocities, respectively; T is the temperature; θ is the potential temperature; p is the pressure; ρ is the air density; g is gravitational acceleration; and f is the Coriolis parameter.
Each variable is separated into the time mean (-) and the deviation from the time mean ('). The governing equations [Eqs. (4)-(8)] for the time-mean field are given by \begin{eqnarray} \label{eq8} &\overline{u_t}\!+\!\bar {u}\overline{u_x}\!+\!(\overline{u_y}-f)\bar {v}+\bar {w}\overline{u_z}\!=\!-\frac{1}{\bar {\rho}}\bar {p}_x -(\overline{u'^2})_x-(\overline{u'v'})_y\!-\!(\overline{u'w'})_z ,\ \ (8)&\nonumber\\ \\ \label{eq9} &\overline{v_t}\!+\!\bar {u}(\overline{v_x}+f)\!+\!\bar {v}\overline{v_y}+\bar {w}\overline{v_z}\!=\!-\frac{1}{\bar {\rho}}\bar {p}_y -(\overline{u'v'})_x-(\overline{v'^2})_y\!-\!(\overline{v'w'})_z ,\ \ (9)&\nonumber\\ \\ \label{eq10} &\overline{w_t}+\bar {u}\overline{w_x}\!+\!\bar {v}\overline{w_y}\!+\!\bar {w}\overline{w_z}\!=\!-\frac{1}{\bar {\rho}}\bar {p}_z \!-\!g\!-\!(\overline{u'w'})_x\!-\!(\overline{v'w'})_y-(\overline{w'w'})_z ,\ \ (10)&\nonumber\\ \\ \label{eq11} &\overline{u_x}+\overline{v_y}+\overline{w_z}\!=\!0 ,\ \ (11)&\\ \label{eq12} &\overline{\theta_t}+\bar {u}\overline{\theta_x}+\bar {v}\overline{\theta_y}+\bar {w}\overline{\theta_z} \!=\!-(\overline{u'\theta'})_x-(\overline{\theta'v'})_y-(\overline{\theta'w'})_z ,\ \ (12)& \end{eqnarray} where the subscripts represent the partial derivatives of time and space. The last three terms on the right-hand side of Eqs. (9)-(11) are stress terms representing the feedback from the disturbance field to the mean field.
By modifying the equations of Kinoshita and Sato (2013a,b), the residual mean circulation (u*,v*,w*) satisfying the continuity equation here is defined as \begin{eqnarray} \label{eq13} \overline{u^\ast}&=&\bar {u}+\frac{(\overline{S})_y}{f}-\left(\frac{\overline{u'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_z ,\ \ (13)\\ \label{eq14} \overline{v^\ast}&=&\bar {v}-\frac{(\overline{S})_x}{f}-\left(\frac{\overline{v'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_z ,\ \ (14)\\ \label{eq15} \overline{w^\ast}&=&\bar {w}+\left(\frac{\overline{u'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_x +\left(\frac{\overline{v'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_y , \ \ (15)\end{eqnarray} where $\overline{S}=\{\overline{u'^2}+\overline{v'^2}+\overline{w'^2}-[(\overline{\theta'^2}g^2)/(N^2\bar{\theta}^2)]\}/2$ is the difference between the kinetic and potential energy and $N^2=(\partial\bar{\theta}/\partial z)(g/\bar{\theta})$ is the Brunt-Vaisala frequency. The Coriolis force (f) here is constant. The residual mean circulation presented above contains not only the horizontal gradient of the disturbance energy difference, but also the 3D gradient of the disturbance heat flux. This 3D residual mean circulation is non-hydrostatic and applicable to mesoscale systems, which is different from those of Kinoshita and Sato (2013a,b).
Substituting Eqs. (14) and (15) into Eqs. (9) and (10) yields the non-hydrostatic TEM equations \begin{eqnarray} \label{eq16} \overline{u_t}+\bar {u}\overline{u_x}+\bar {v}\overline{u_y}+\bar {w}\overline{u_z}-f\overline{v^\ast}&=& -\frac{1}{\bar {\rho}}\bar {p}_x-\nabla\cdot {F}_1 ,\ \ (16)\\ \label{eq17} \overline{v_t}+\bar {u}\overline{v_x}+\bar {v}\overline{v_y}+\bar {w}\overline{v_z}+f\overline{u^\ast}&=& -\frac{1}{\bar {\rho}}\bar {p}_y-\nabla\cdot {F}_2 . \ \ (17)\end{eqnarray} And, from Eq. (11) we obtain \begin{equation} \label{eq18} \overline{w_t}+\bar {u}\overline{w_x}+\bar {v}\overline{w_y}+\bar {w}\overline{w_z}= -\frac{1}{\bar {\rho}}\bar {p}_z-g-\nabla\cdot {F}_3 , \ \ (18)\end{equation} where F1=(F1,1,F1,2,F1,3), F2=(F2,1,F2,2,F2,3) and F3=(F3,1,F3,2,F3,3) are the zonal, meridional and vertical components of the 3D wave activity flux. Their expressions are given by \begin{eqnarray} \label{eq19} F_{1,1}&=&\overline{u'^2}-\overline{S} ,\quad F_{1,2}=\overline{u'v'} ,\quad F_{1,3}=\overline{u'w'}-f\frac{\overline{v'\theta'}}{N^2}\frac{g}{\bar {\theta}} ,\ \ (19)\\ \label{eq20} F_{2,1}&=&\overline{u'v'} ,\quad F_{2,2}=\overline{v'^2}-\overline{S} ,\quad F_{2,3}=\overline{v'w'}+f\frac{\overline{u'\theta'}}{N^2}\frac{g}{\bar {\theta}} ,\ \ (20)\qquad\\ \label{eq21} F_{3,1}&=&\overline{u'w'} ,\quad F_{3,2}=\overline{v'w'} ,\quad F_{3,3}=\overline{w'w'} . \ \ (21)\end{eqnarray}
Here, the 3D wave activity flux is extended, including the zonal, meridional and vertical components. The primary difference with the wave activity flux from (Miyahara, 2006) and Kinoshita and Sato (2013a,b) is the vertical momentum equations [Eq. (20)], which are not considered in their formulation.
In contrast to the TEM equations in large-scale systems, this residual mean circulation only appears in the Coriolis force terms of the horizontal TEM equation. In addition, although the disturbance energy difference and the disturbance heat flux are introduced into the equations, they appear on both sides of the equations simultaneously. The TEM equations are based on non-geostrophic and on non-hydrostatic equilibrium; thus, they are sufficient to describe the active interaction between the mesoscale systems and the background circulation.
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.
Figure1. (a) Model domains. The shading denotes the topography (units: km), and the red line denotes the position of the cross-section discussed below. (b, c) 6-h accumulated precipitation (units: mm) according to (b) observation and (c) simulation data, during 1200-1800 UTC 17 August 2014.2
4.1. Numerical model description
The numerical experiment is performed using the Weather Research and Forecasting model to simulate the rainstorm. The initial and lateral boundary conditions of the model are based on the Global Forecast System (GFS) analysis dataset from the National Centers for Environmental Prediction, with horizontal grid spacing of 0.5°× 0.5°. The numerical experiment is one-way nested with two domains (Fig. 1a). The two domains both have 303× 303 horizontal grids, with horizontal grid spacings of 3 km and 1 km. Also, the domains have 51 vertical layers, with the model top set at 50 hPa. The model is integrated for 36 h, starting at 0000 UTC 17 August 2014, with output every 10 min. The physical parameterization schemes include the New Thompson microphysics scheme (Thompson et al., 2008), the RRTM longwave radiation transfer scheme (Mlawer et al., 1997), the Goddard shortwave radiation transfer scheme, and the Mellor-Yamada-Janjic scheme of the planetary boundary layer (Janji?, 1994).Figures 1b and c show the 6-h accumulated precipitation from the observations and the simulation in domain 1. The observed precipitation data are an hourly precipitation gridded dataset, which is from China automatic station and CMORPH precipitation data (http://data.cma.cn/data/detail/dataCode/SEVP_CLI_CHN_MERGE_CMP_PRE_HOUR_GRID_0.10/keywords/cmorph.html ). The heavy rainfall occurred predominantly in Sichuan, along the border of the Sichuan Basin and the Tibetan Plateau. The rainband moved from the northwest to the southeast. Compared to the observations, the simulated precipitation is slightly stronger, possibly because the observational data are from meteorological stations, but the terrain is complex in northern Sichuan and in the Tibetan Plateau; thus, there are fewer meteorological stations and the observational data have a lower resolution. However, the simulation can accommodate changes in the rainfall location, variation and duration. Next, the simulation data are used to study the 3D wave activity flux during the precipitation event.
2
4.2. Inertia-gravity waves
Lu et al. (2005a,b) and (Gao et al., 2015) discussed that the polarization——meaning that the directions of wave propagation and wave oscillation are perpendicular——is a representative characteristic of inertia-gravity waves. The phase difference between the perturbation divergence and vorticity of π/2 is an important characteristic of inertia-gravity waves. This polarization feature is crucial for extracting the inertia-gravity waves simply from simulations and observational data.Figures 2a and b show the vertical distribution of vertical vorticity and horizontal divergence at 1030 UTC and 1400 UTC. There are clearly positive and negative changes in the vertical vorticity and horizontal divergence at the height of 12 km in the rainstorm system. Also, it can be seen that the divergence region is one-quarter of a wavelength ahead of the vorticity, which is a typical characteristic of inertia-gravity waves. A detailed analysis of the characteristics of inertia-gravity waves during the heavy precipitation event can be found in (Liu et al., 2018). Notably, when inertia-gravity waves gradually enhance, the corresponding surface precipitation is also increasing. At 1400 UTC (Fig. 2b), the strength of the inertia-gravity waves and the accumulated precipitation peak almost simultaneously. It can be seen that the precipitation has a close relationship with the inertia-gravity waves. Therefore, in the next section, the 3D EP flux is introduced to analyze the influence of inertia-gravity waves on the development of precipitation.
Figure2. Vertical profile of vertical vorticity (shaded; units: 10-3 s-2) and horizontal divergence (contours; units: 10-3 s-2) along 31.1°N at (a) 1030 UTC and (b) 1400 UTC. The grey region denotes the topography. The red line represents the 30-min accumulated precipitation (units: mm).2
4.3. 3D EP flux analysis
In this section, to demonstrate the application of the divergence of the 3D EP flux, a simple case of a heavy precipitation event that occurred in the Sichuan area is analyzed. The time-mean field is obtained by a moving average with a period of 100 min, which is approximately the period of the inertia-gravity waves during the precipitation (Liu et al., 2018). The disturbance field is defined as the deviation from the time-mean field.Figure 3 is a vertical profile of area-average $\bar u_t$ (the averaged area is the domain 2 region). The red line represents the $\bar u_t$ resulting from the time partial derivative of the simulated u-wind, and the green line represents the $\bar u_t$ resulting from the left-hand side of the TEM equation [Eq. (15)]. Although the red line is slightly stronger than the green line, the two lines have the same variation tendency with height: positive in the lower and higher level, indicating that the u-wind tends to increase, and negative in the middle level, meaning that the u-wind tends to decrease. Thus, the TEM equation is valid and can be used in the following.
Figure3. Vertical profile of area-average $\bar{u}_t$ at 0900 UTC (the red line represents the time partial derivative from simulated data; the green line is calculated from the TEM equation; units: m s-2).The EP flux is proportional to the group velocity deduced above. Figure 4 shows the longitude-pressure cross section of the EP flux F1 and the EP flux divergence ?· F1 during the mature stages of the rainstorm. The contour is the 30-min accumulated precipitation, denoting the rainstorm location. The maximum of the EP flux divergence is primarily concentrated at the height of 10-14 km, which is consistent with the above-mentioned location of the inertia-gravity waves. In terms of the TEM equation, the EP flux divergence corresponds to the wave forcing to the mean flow. The positive EP flux divergence could weaken the mean flow, and the negative EP flux divergence could enhance the mean flow. There were two parts of the EP flux divergence in the rainstorm. One part is in the region of 103°-104°E, with negative EP flux divergence, indicating that the EP flux divergence accelerates the westerly wind at the height of 12 km. The other is in the region of 104°-105°E, where the EP flux divergence alternates between positive and negative. This means that there is almost no acceleration or deceleration. The EP flux (vector) is proportional to the group velocity, representing the propagation of the inertia-gravity wave energy. The EP flux vector in Fig. 4, with a maximum value in the region of 103°-104°E, points forward (eastward) and upward in the rainstorm, indicating that the wave energy propagates forward (eastward) and upward. The reason for the forward and upward propagation of the wave energy can be explained insofar as the strong westerly airflow runs over the mountain, causing powerful ascending velocity due to the mountain lift. The powerful upward motion transports the wave energy upward and the strong westerly airflow transports the wave energy forward. Whereas, in the decaying stage, the EP flux and the EP flux divergence almost vanish (figure omitted).
Figure4. Cross sections of EP flux (vectors) and EP flux divergence (shaded; units: 10-3 m s-2) along 31.1°N at 1030 UTC. The red line denotes the 30-min accumulated precipitation (units: mm).To examine the contribution of each term of the TEM equation to the mean flow, the longitude-pressure cross sections of $\overline{V}\cdot\nabla\bar{u}$, $\overline{v^{\ast}}f$ and ?· F1 are shown in Fig. 5. The EP flux divergence ?· F1 is strongest (Fig. 5c) and plays a key role in the environmental flow variation. The EP flux divergence is negative in the region of 103°-104°E, indicating an acceleration of westerly wind. There is a high-level jet at the height of 12-14 km during the rainstorm, and therefore the negative EP flux divergence is maintained and enhances the high-level jet, which is conducive to enhancing the vertical wind shear. In (Liu et al., 2018), it was shown that intensified high-level vertical wind shear favors inertia-gravity wave development and enhances convection, thus promoting rainstorm development. Many researchers (e.g., Uccelini, 1975; Li, 1978; Koch and Siedlarz, 1999; Lane and Zhang, 2011) have examined the interaction between inertia-gravity waves and deep convection and proposed that mesoscale gravity waves are capable of causing convective instability and initiating thunderstorm development; in turn, the latent heat release feeds back to the background waves, triggering or amplifying the gravity waves. In addition, there are positive and negative changes in the EP flux divergence in the region of 104°-105°E, indicatingthe westerly wind accelerates and decelerates alternately. The acceleration and deceleration cancel each other out, meaning there is almost no acceleration of the westerly wind. Therefore, strong convergence of EP flux occurs near the region of 103.5°-104°E, providing beneficial conditions for heavy precipitation development. The wind advection term $\overline{V}\cdot\nabla\bar{u}$ is the second major contributor to the mean flow (Fig. 5a). The wind advection term is positive at the height of 12-14 km and negative at the height of 14-16 km, indicating deceleration of the westerly wind at relatively low levels and acceleration at relatively high levels. The residual mean circulation term $\overline{v^{\ast}}f$ (Fig. 5b) is much weaker, with a distribution at height of 12-14 km.
Figure5. Cross sections of terms of the 3D TEM equation at 1030 UTC along 31.1°N: (a) $\bar{v}\cdot\nabla\bar{u}$; (b) $\overline{v^{\ast}}f$; (c) ?· F1 (units: 10-3 m s-2). The black line denotes the 30-min accumulated precipitation (units: mm).Figure 6 shows the vertical profile of the regional average of $\overline{V}\cdot\nabla\bar{u}$, $\overline{v^{\ast}}f$, ?· F1 and $\bar{u}_t$. The selected region is where the 6-h accumulated precipitation is greater than 10 mm. During the developing stage of the rainstorm (Fig. 6a), the EP flux divergence is negative at the height of 12-14 km, which is advantageous to enhancing the westerly wind and the high-level jet. Also, the EP flux divergence is positive at approximately the height of 10 km and 15 km, which is conducive to weakening the westerly wind above and below the high-level jet, thus enhancing the vertical wind shear of the environment. The non-geostrophic balance and vertical wind shear near the high-level jet is conducive to inertia-gravity wave enhancement; in turn, the developing inertia-gravity waves enhance the perturbation of the atmosphere, and thus the EP flux markedly strengthens. The residual mean circulation term $\overline{v^{\ast}}f$ is weaker and almost opposite to the distribution of the wind advection term $\overline{V}\cdot\nabla\bar{u}$, but their effects are the same with the deceleration of westerly wind. In the declining stage (Fig. 6b), each term is much weaker than before. The distribution of the EP flux is unchanged, but its intensity is obviously reduced. The wind advection terms contributes the most to the wind tendency $\bar{u}_t$, causing the reduction of the westerly wind.
Figure6. Vertical profile of area-average terms from the 3D TEM equation (a) at 1030 UTC and (b) at 1500 UTC. The red line denotes the EP flux divergence ? · F1; the green line denotes the advection terms $\bar{v}\cdot\nabla\bar{u}$; the blue line denotes the residual mean circulation $\overline{v^{\ast}}f$; and the yellow line denotes the local variation $\bar{u}_t$ (units: 10-3 m s-2).To study the causes of the EP flux divergence variation, we disassemble the 3D EP flux divergence into three parts. Figure 7 shows the longitude-pressure cross sections of ? F1,1/? x, ? F1,2/? y and ? F1,3/? z, where F1,1 denotes the zonal transport of the zonal momentum, F1,2 denotes meridional transport of zonal momentum, and F1,3 denotes the vertical transport of zonal momentum and the meridional transport of the heat flux. It can be seen that ? F1,1/? x and ? F1,3/? z are much stronger than ? F1,2/? y. Comparison between ? F1,1/? x, ? F1,3/? z and ? · F1 shows that ? F1,1/? x makes the greatest contribution to ? · F1, and their distributions are similar. F1,1 consists of the zonal transport of horizontal momentum $\overline{u'u'}$ and the difference between the perturbation kinetic energy and perturbation potential energy S, which are shown in Figs. 8a and b. The two terms with the same magnitude are both important to F1,1. $\overline{u'u'}$ is primarily located over the peak of the precipitation, possibly due to the high-level jet. S is positive to the east of the maximal precipitation and negative to the west, meaning the perturbation kinetic energy is stronger to the east of the maximal precipitation and the perturbation potential energy is stronger at the west. F1,3 consists of the vertical transport of horizontal momentum $\overline{u'w'}$ and the meridional transport of the heat flux $f(\overline{v'\theta'}/N^2)(g/\bar{\theta})$, which are shown in Figs. 8c and d. $\overline{u'w'}$ is much stronger than $f(\overline{v'\theta'}/N^2)(g/\bar{\theta})$, so the vertical transport of horizontal momentum $\overline{u'w'}$ is the primary contributor to the F1,3 variation. This is because there is strong upward motion over the peak of the precipitation due to the mountain uplift. Previous studies have shown that the vertical transportation of horizontal momentum flux is very important to convective activity (e.g., Houze, 2004; Chong et al., 1987; Chong, 2010). This study provides another explanation for the importance of $\overline{u'w'}$.
Figure7. Cross sections of (a) ? F1,1/? x, (b) ? F1,2/? y and (c) ? F1,3/? z (units: m-2 s-2) at 1030 UTC along 31.1°N. The red line denotes the 30-min accumulated precipitation (units: mm).
Figure8. Cross sections of (a) $\overline{u'u'}$, (b) S, (c) $\overline{u'w'}$ and (d) $f(\overline{v'\theta'}/N^2)(g/\bar{\theta})$ (units: m-2 s-2) at 1030 UTC along 31.1°N. The red line denotes the 30-min accumulated precipitation (units: mm).A case of heavy precipitation that occurred in the Sichuan area during 0600-1800 UTC 17 August 2014 is taken as an example to study the formulation of the 3D TEM equations and the 3D EP flux deduced here. Examining the contribution of each term of the TEM equation to the mean flow, the EP flux divergence ?· F1 is identified as the primary contributor to the mean flow tendency. During the mature stage of the precipitation, the maximum EP flux divergence is primarily concentrated at the height of 10-14 km, where the energy of inertia-gravity waves propagates forward (eastward) and upward. The EP flux divergence is negative at the height of 10 km and 15 km, decelerating the u-wind above and below the high-level jet, and positive at the height of 12-14 km, accelerating the high-level jet. This distribution of the EP flux divergence maintains and enhances the high-level jet, thus strengthening the vertical wind shear of the environment and promoting rainstorm development.
Finally, the factors influencing EP flux divergence are examined, revealing that the zonal component ? F1,1/? x and the vertical component ? F1,3/? z make the greatest contribution to ?· F1. The zonal transport of horizontal momentum $\overline{u'u'}$ and the difference value between the perturbation kinetic energy and perturbation potential energy S are both crucial to F1,1, and the vertical transport of horizontal momentum $\overline{u'w'}$ plays the key role in the variation of F1,3. From the above, the EP flux divergence has a significant influence on environmental flow, thus promoting the development of the heavy precipitation.
However, the diabatic effect in mesoscale systems is also very important to rainfall development, but the formulation of the 3D TEM equation deduced here does not yet consider its influence. This will be addressed in future work.
3
APPENDIX A
From the governing equations, it is assumed that the basic flow is still, and the linear disturbance equations for 3D inertia-gravity waves are \begin{eqnarray} \frac{\partial{u}'}{\partial t}-f{v}'+\frac{1}{\bar {\rho}}\frac{\partial{p}'}{\partial x}=0 ,\ \ (A1)\\ \frac{\partial{v}'}{\partial t}+f{u}'+\frac{1}{\bar {\rho}}\frac{\partial{p}'}{\partial y}=0 , \ \ (A2)\end{eqnarray} \begin{eqnarray} \frac{\partial{w}'}{\partial t}+\frac{1}{\bar {\rho}}\frac{\partial{p}'}{\partial z} -g\frac{{\theta}'}{\bar {\theta}}=0 , \ \ (A3)\\ \frac{\partial{u}'}{\partial x}+\frac{\partial{v}'}{\partial y}+\frac{\partial{w}'}{\partial z}=0 , \ \ (A4)\\ \frac{\partial{\theta}'}{\partial t}+N^2\frac{\bar {\theta}}{g}{w}'=0 . \ \ (A5)\end{eqnarray} According to multi-scale methods of waves, the slow spatial and temporal scale is defined as X=ε x, Y=ε y, Z=ε z, T=ε t, where x, y, and t are the physical quantities of the quick spatial and temporal scale, and parameter |ε|? 1.It is assumed that a form of a plane wave is considered, \begin{equation} {A}'=A_0(X,Y,Z,T){\rm e}^{i\phi} , \ \ (A6)\end{equation} where φ=kx+ly+mz-ω t is the wave phase; k, l and m are the wave numbers in the x, y, and z directions; $\tilde\omega=\omega-k\bar{u}$ is the intrinsic frequency; and ω is the local frequency. Because the basic flow is still $(\bar{u}\approx0)$,the intrinsic frequency is thus equal to the local frequency. The phase of φ satisfies the relationships \begin{equation} \frac{\partial\phi}{\partial t}=-\omega ,\quad \frac{\partial\phi}{\partial x}=k ,\quad \frac{\partial\phi}{\partial y}=l ,\quad \frac{\partial\phi}{\partial z}=n . \ \ (A7)\end{equation}
Then, it can be proved that the wave numbers satisfy the following relationships: \begin{eqnarray} \frac{\partial\omega}{\partial X}&=&-\frac{\partial k}{\partial T} ,\quad \frac{\partial\omega}{\partial Y}=-\frac{\partial l}{\partial T} ,\quad \frac{\partial\omega}{\partial Z}=-\frac{\partial m}{\partial T} ,\nonumber\\ \frac{\partial k}{\partial Y}&=&\frac{\partial l}{\partial X} ,\quad \frac{\partial k}{\partial Z}=\frac{\partial m}{\partial X} ,\quad \frac{\partial l}{\partial Z}=\frac{\partial m}{\partial Y} . \ \ (A8)\end{eqnarray}
The amplitude of waves is expressed by the method of small parameters: \begin{equation} A_0(X,Y,Z,T)=\sum_{n=0}^\infty\varepsilon_nA_{0,n}(X,Y,Z,T) . \ \ (A9)\end{equation}
Substitution of Eq. (A9) into Eqs. (A1)-(A5) yields the ε0 approximate equations: \begin{eqnarray} -i\omega u_{0,0}-fv_{0,0}+ik\frac{p_{0,0}}{\bar {\rho}}=0 ,\ \ (A10)\\ -i\omega v_{0,0}+fu_{0,0}+il\frac{p_{0,0}}{\bar{\rho}}=0 ,\ \ (A11)\\ -i\omega w_{0,0}+in\frac{p_{0,0}}{\bar {\rho}}-g\frac{\theta_{0,0}}{\bar{\theta}}=0 ,\ \ (A12)\\ ku_{0,0}+lv_{0,0}+mw_{0,0}=0 ,\ \ (A13)\\ -i\omega\theta_{0,0}+N^2\frac{\bar {\theta}}{g}w_{0,0}=0 , \ \ (A14)\end{eqnarray} and the ε1 approximate equations: \begin{eqnarray} -i\omega u_{1,0}-fv_{1,0}+ik\frac{p_{1,0}}{\bar {\rho}}=-\left(\frac{\partial u_{0,0}}{\partial T}+ \frac{1}{\bar {\rho}}\frac{\partial p_{0,0}}{\partial X}\right)&\!=\!&-A ,\ \ (A15)\qquad\ \ \\ -i\omega v_{1,0}+fu_{1,0}+il\frac{p_{1,0}}{\bar {\rho}}=-\left(\frac{\partial v_{1,0}}{\partial T}+ \frac{1}{\bar {\rho}}\frac{\partial p_{0,0}}{\partial Y}\right)&\!=\!&-B ,\ \ (A16)\\ -i\omega w_{1,0}+in\frac{p_{1,0}}{\bar {\rho}}-g\frac{\theta_{1,0}}{\bar {\theta}} =-\left(\frac{\partial w_{0,0}}{\partial T}+\frac{1}{\bar {\rho}}\frac{\partial p_{0,0}}{\partial Z}\right)&\!=\!&-C , \ \ (A17)\end{eqnarray} \begin{eqnarray} ku_{1,0}+lv_{1,0}+mw_{1,0}=i\frac{\partial u_{0,0}}{\partial X}+i\frac{\partial v_{0,0}}{\partial Y}+i\frac{\partial w_{0,0}}{\partial Z}&\!=\!&-D ,\ \ (A18)\qquad\quad\\ -i\omega\theta_{1,0}+N^2\frac{\bar {\theta}}{g}w_{1,0}=-\frac{\partial\theta_{0,0}}{\partial T}&\!=\!&-E .\ \ (A19) \end{eqnarray}
Eliminating the zero-order disturbance in the ε1 approximate equations produces \begin{eqnarray} &&(i\omega k+fl)mA+(i\omega l-fk)mB-i\omega(k^2+l^2)C-\qquad\nonumber\\ &&(\omega^2-f^2)mD+(k^2+l^2)\frac{g}{\bar {\theta}}E=0 . \ \ (A20)\end{eqnarray}
Substituting Eqs. (A15)-(A19) of A, B, C, D and E into Eq. (A20), the 3D wave activity equations for inertia-gravity waves are obtained: \begin{equation} \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 , \ \ (A21)\end{equation} where H=(k2+l2+m2)ω02 is the wave activity density and Cg=(Cg,x,Cg,y,Cg,z) are the group velocities of the 3D inertia-gravity waves: \begin{eqnarray} C_{{g},x}&=&\frac{(\omega^2-f^2)n^2k}{\omega(k^2+l^2+n^2)(k^2+l^2)} ,\ \ (A22)\\ C_{{g},y}&=&\frac{(\omega^2-f^2)n^2l}{\omega(k^2+l^2+n^2)(k^2+l^2)} ,\ \ (A23)\\ C_{{g},z}&=&\frac{(\omega^2-f^2)n}{\omega(k^2+l^2+n^2)} , \ \ (A24)\end{eqnarray} H satisfies the wave activity density conservation equation as follows: \begin{equation} \frac{\partial A}{\partial t}+\nabla\cdot({C}_{g}A)=0 . \ \ (A25)\end{equation} Based on the dispersion relation for inertia-gravity waves [Eq. (32)] and the local wave number [Eq. (A8)], the wave parameter tendency equations are given in the following: \begin{eqnarray} \frac{\partial\omega}{\partial T}+C_{{g},x}\frac{\partial\omega}{\partial X}+ C_{{g},y}\frac{\partial\omega}{\partial Y}+C_{{g},z}\frac{\partial\omega}{\partial Z}&=&0 ,\ \ (A26)\\ \frac{\partial k}{\partial T}+C_{{g},x}\frac{\partial k}{\partial X}+C_{{g},y}\frac{\partial k}{\partial Y}+C_{{g},z}\frac{\partial k}{\partial Z}&=&0 ,\ \ (A27)\\ \frac{\partial l}{\partial T}+C_{{g},x}\frac{\partial l}{\partial X}+C_{{g},y}\frac{\partial l}{\partial Y}+C_{{g},z}\frac{\partial l}{\partial Z}&=&0 ,\ \ (A28)\\ \frac{\partial m}{\partial T}+C_{{g},x}\frac{\partial m}{\partial X}+C_{{g},y}\frac{\partial m}{\partial Y}+C_{{g},z}\frac{\partial m}{\partial Z}&=&0 . \ \ (A29)\end{eqnarray}
3
APPENDIX B
Andrews and McIntyre (1976, 1978) noted that the residual mean circulation has to satisfy the relationship between the Eulerian-mean flow and the Stokes drift. Kinoshita and Sato (2013a,b studied and demonstrated that their 3D residual mean circulation is the sum of the Eulerian-mean flow and the Stokes drift under different assumptions at large scales. In this section, considering disturbances with small amplitudes in the basic state, we demonstrate that the residual mean circulation satisfies the relationship between the Eulerian-mean flow and the Stokes drift under non-hydrostatic equilibrium in mesoscale systems.The Stokes drift describes the difference between the Lagrangian and Eulerian averages of flow fields: \begin{equation} \bar {u}_{\rm S}=\bar {u}_{\rm L}-\bar {u} , \ \ (B1)\end{equation} where $\bar{u}_L$ is the Lagrangian average and $\bar{u}$ is the Eulerian average. The generic expression for Stokes drift by Taylor-expansion is \begin{equation} \bar {u}_{\rm S}=\overline{\xi'u'_x}+\overline{\eta'v'_y}+\overline{\varsigma'u'_z}+o(a^2) , \ \ (B2)\end{equation} where $\xi'(x,y,z,t)$, η'(x,y,z,t) and ζ'(x,y,z,t) are the zonal, meridional and vertical displacement of the air parcel, respectively. Ignoring the third and higher order, the 3D Stokes drift can be written as \begin{eqnarray} \bar {u}_{\rm S}&=&(\overline{\xi'u'})_x+(\overline{\eta'u'})_y+(\overline{\zeta'u'})_z ,\ \ (B3)\\ \bar{v}_{\rm S}&=&(\overline{\xi'v'})_x+(\overline{\eta'v'})_y+(\overline{\zeta'v'})_z ,\ \ (B4)\\ \bar {w}_{\rm S}&=&(\overline{\xi'w'})_x+(\overline{\eta'w'})_y+(\overline{\zeta'w'})_z , \ \ (B5)\end{eqnarray} where $\xi'$, η' and ζ' satisfy the Lagrangian constraint relationships \begin{eqnarray} &&{\xi}'_x+{\eta}'_y+{\zeta}'_z=0 ,\ \ (B6)\\ &&\overline{{\xi}'}=\overline{{\eta}'}=\overline{{\zeta}'}=0 ,\ \ (B7)\\ &&\overline{D}{\xi}'={u}' ,\quad \overline{D}{\eta}'={v}' ,\quad \overline{D}{\zeta}'={w}' , \ \ (B8)\end{eqnarray} and \begin{equation} \overline{D}=\frac{\partial}{\partial t}+\bar {u}\frac{\partial}{\partial x}+ \bar {v}\frac{\partial}{\partial y}+\bar {w}\frac{\partial}{\partial z} \ \ (B9)\end{equation}
Ignoring basic flow, Eq. (B8) is expressed as follows \begin{equation} \frac{\partial{\xi}'}{\partial t}={u}' ,\quad \frac{\partial{\eta}'}{\partial t}={v}' ,\quad \frac{\partial{\zeta}'}{\partial t}={w}' . \ \ (B10)\end{equation}
It is assumed that the particle disturbance displacement has the form of a plane sinusoidal wave $(\xi',\eta',\zeta')=(\xi_0,\eta_0,\zeta_0)\rm e^i\phi$, which is substituted into Eq. (B10). The particle displacement amplitude is obtained as \begin{eqnarray} \xi_0&=&\frac{i}{\omega}u_0=i\frac{k\omega p_{{\rm 0r}}-lfp_{{\rm 0i}}}{\omega(\omega^2-f^2)\rho_0}-\frac{k\omega p_{{\rm 0i}}+lfp_{{\rm 0r}}}{\omega(\omega^2-f^2)\rho_0} ,\ \ (B11)\\ \eta_0&=&\frac{iv_0}{\omega}=i\frac{kfp_{{\rm 0i}}-l\omega p_{{\rm 0r}}}{\omega(\omega^2-f^2)\rho_0}-\frac{l\omega p_{{\rm 0i}}-kfp_{{\rm 0r}}}{\omega(\omega^2-f^2)\rho_0} ,\ \ (B12)\\ \zeta_0&=&\frac{iw_0}{\omega}=-i\frac{(k^2+l^2)\omega p_{0{\rm r}}}{\omega n(\omega^2-f^2)\rho_0}+\frac{(k^2+l^2)\omega p_{{\rm 0i}}}{\omega n(\omega^2-f^2)\rho_0} .\ \ (B13)\quad \end{eqnarray} By using the particle displacement amplitudes, Eqs. (B11)-(B13), and the disturbance amplitudes, Eqs. (25)-(27), along with the Reynolds average relationship, it can be proved that the products of $\overline{\xi'u'}$, $\overline{\eta'v'}$ and $\overline{\zeta'w'}$ are zero because two quantities of the respective products are out of phase by 90°: $\overline{\xi'u'}= (1/2)\rm R(\xi_0^\otimes u_0)=0$, $\overline{\eta'v'}=(1/2)\rm R(\eta_0^\otimes v_0)=0$, $\overline{\zeta'w'}=(1/2)\rm R(\zeta_0^\otimes w_0)=0$.
From Eqs. (B3)-(B5), the 3D Stokes drift is obtained: \begin{eqnarray} \bar {u}_{\rm S}&=&(\overline{\eta'u'})_y+(\overline{\zeta'u'})_z ,\ \ (B14)\\ \bar {v}_{\rm S}&=&(\overline{\xi'v'})_x+(\overline{\zeta'v'})_z ,\ \ (B15)\\ \bar {w}_{\rm S}&=&(\overline{\xi'w'})_x+(\overline{\eta'w'})_y . \ \ (B16)\end{eqnarray}
By using the particle displacement amplitudes, Eqs. (B11)-(B13) and the Reynolds average relationship, the following relation is obtained: \begin{eqnarray} \overline{\eta'u'}&=&\frac{f(k^2+l^2)}{(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B17)\\ \overline{\xi'v'}&=&-\frac{f(k^2+l^2)}{(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B18)\\ \overline{\varsigma'u'}&=&-\frac{lf(k^2+l^2)}{n(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B19)\\ \overline{\xi'w'}&=&\frac{lf(k^2+l^2)}{n(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B20)\\ \overline{\eta'w'}&=&-\frac{kf(k^2+l^2)}{n(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B21)\\ \overline{\varsigma'v'}&=&\frac{kf(k^2+l^2)}{n(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B22)\\ \overline{u'\frac{\theta'}{N^2}\frac{g}{\bar {\theta}}}&=&\frac{fl(k^2+l^2)}{n(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B23)\\ \overline{v'\frac{\theta'}{N^2}\frac{g}{\bar {\theta}}}&=&-\frac{kf(k^2+l^2)}{n(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} ,\ \ (B24)\\ \overline{S}&=&\frac{f^2(k^2+l^2)}{(f^2-\omega^2)^2}\frac{\overline{p'^2}}{\bar {\rho}^2} . \ \ (B25)\end{eqnarray}
Substituting Eqs. (B17)-(B22) into the expression of Stokes drift, Eqs. (b14)-(b16), yields \begin{eqnarray} \bar {u}_{\rm S}&=&\frac{(\overline{S})_y}{f}-\left(\frac{\overline{u'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_z ,\ \ (B26)\\ \bar {v}_{\rm S}&=&-\frac{(\overline{S})_x}{f}-\left(\frac{\overline{v'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_z ,\ \ (B27)\\ \bar {w}_{\rm S}&=&\left(\frac{\overline{u'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_x+\left(\frac{\overline{v'\theta'}}{N^2}\frac{g}{\bar {\theta}}\right)_y . \ \ (B28)\end{eqnarray}
Comparing the expressions of Stokes drift, Eqs. (B26)-(B28), with the expressions of the residual mean circulation, Eqs. (12)-(14), the following relation is obtained: \begin{eqnarray} \overline{u^\ast}&=&\bar {u}+\bar {u}_{\rm S} ,\ \ (B29)\\ \overline{v^\ast}&=&\bar {v}+\bar {v}_{\rm S} ,\ \ (B30)\\ \overline{w^\ast}&=&\bar {w}+\bar {w}_{\rm S} . \ \ (B31)\end{eqnarray}
Thus, 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 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.
