Fast-Wave Averaging with Phase Changes: Asymptotics and Application to Moist Atmospheric Dynamics

Abstract

Many systems involve the coupled nonlinear evolution of slow and fast components, where, for example, the fast waves might be acoustic (sound) waves with a small Mach number or inertio-gravity waves with small Froude and Rossby numbers. In the past, for some such systems, an interesting property has been shown: the slow component actually evolves independently of the fast waves, in a singular limit of fast wave oscillations. Here, a fast-wave averaging framework is developed for a moist Boussinesq system with additional complexity beyond past cases, now including phase changes between water vapor and liquid water. The main question is: Do phase changes induce coupling between the slow component and fast waves? Or does the slow component evolve independently, according to moist quasi-geostrophic equations? Compared to the dry dynamics, a substantial challenge is that the method needs to be adapted to a piecewise operator with variable coefficients, due to phase changes. A formal asymptotic analysis is presented here. For purely saturated flow without phase changes, it is shown that precipitation does not induce coupling, and the slow modes evolve independently. With phase changes present, the limiting equations show that phase boundaries could possibly induce coupling between the slow modes and fast waves.

Appendices

Appendix A: Non-dimensional Equations and Distinguished Limit

For the moist Boussinesq equations with phase changes, the dimensional form is shown in (1a)–(1d) of Smith and Stechmann (2017), and a non-dimensional version is described in the appendix of Smith and Stechmann (2017) in terms of buoyancy variables $b_\mathrm{u}$ and $b_\mathrm{s}$. Here, a different, but equivalent, non-dimensional version is described, using $\theta _e$ and $q_t$ as the moist thermodynamic variables:

$$\begin{aligned}&\dfrac{D_h\vec {u}_h}{Dt}+w\dfrac{\partial u_h}{\partial z}+R_0^{-1}u_h^{\bot }+E_\mathrm{u}\nabla _h\phi =0 \end{aligned}$$

(A.1)

$$\begin{aligned}&A^2\left( \dfrac{D_hw}{Dt}+w\dfrac{\partial w}{\partial z}\right) +E_\mathrm{u}\dfrac{\partial \phi }{\partial z}-\Gamma A^2 b=0 \end{aligned}$$

(A.2)

$$\begin{aligned}&\nabla _h\cdot u_h+\dfrac{\partial w}{\partial z}=0 \end{aligned}$$

(A.3)

$$\begin{aligned}&\dfrac{D_h\theta _e}{Dt}+w\dfrac{\partial \theta _e}{\partial z}+{\mathrm{Fr}_1}^{-2}\left( \Gamma A^2\right) ^{-1}w=0 \end{aligned}$$

(A.4)

$$\begin{aligned}&\dfrac{D_h q_t}{Dt}+w\dfrac{\partial q_t}{\partial z}-{\mathrm{Fr}_2}^{-2}\left( \Gamma A^2\right) ^{-1}w-V_rC_{cl}\dfrac{\partial q_r}{\partial z}=0 \end{aligned}$$

(A.5)

along with the relationships

$$\begin{aligned} b = b_\mathrm{u}H_\mathrm{u} + b_\mathrm{s}H_\mathrm{s}, \quad b_\mathrm{u} = \theta _e + \left( \dfrac{c_p \theta _0}{L_v}-1\right) q_t, \quad b_\mathrm{s} = \theta _e - \dfrac{c_p \theta _0}{L_v}q_t, \end{aligned}$$

(A.6)

where $(\mathrm{Ro}, \mathrm{Eu}, A, \Gamma , V_r)$ are the Rossby number, Euler number, aspect ratio, buoyancy parameter, and rain fall speed, respectively. Note that there are two moist thermodynamic variables ($\theta _e$ and $q_t$) and two phases, as opposed to the dry case with one Froude number, one thermodynamic variable ($\theta$), and one phase. The two “Froude” numbers used here are

$$\begin{aligned} \mathrm{Fr}_1&=U(N_1H)^{-1}\quad L_{d_1}=\dfrac{N_1H}{f}, \end{aligned}$$

(A.7)

$$\begin{aligned} \mathrm{Fr}_2&=U(N_2H)^{-1}\quad L_{d_2}=\dfrac{N_2H}{f}, \end{aligned}$$

(A.8)

$$\begin{aligned} {N_1}^2&= \dfrac{g}{\theta _0} \dfrac{\hbox {d} \tilde{\theta }_e }{\hbox {d}z} = \dfrac{g}{\theta _0} \dfrac{\hbox {d}}{\hbox {d}z} \left( \tilde{\theta } + \dfrac{L_v}{c_p}\tilde{q}_{v}\right) = \dfrac{g}{\theta _0}\left( B + \dfrac{L_v}{c_p}B_{vs}\right) , \end{aligned}$$

(A.9)

$$\begin{aligned} {N_2}^2&= - \dfrac{g}{\theta _0} \dfrac{L_v}{c_p} \dfrac{\hbox {d}\tilde{q}_t }{\hbox {d}z} = -\frac{g}{\theta _0}\left( \frac{L_v}{c_p}B_{vs}\right) , \end{aligned}$$

(A.10)

where $L_{d1}$ and $L_{d2}$ are Rossby radii of deformation, and $N_1$ and $N_2$ are buoyancy frequencies. Note that the notation $\mathrm{Fr}_2$, $L_{d2}$, $N_2$ is used in analogy to Froude number, Rossby radius of deformation, and buoyancy frequency, respectively, although $\mathrm{Fr}_2$, $L_{d2}$, and $N_2$ are defined in terms of not buoyancy but total water. More detail information of reference scales and the non-dimensional quantities can be found in Smith and Stechmann (2017) (Table A1, Table A2).

To define the distinguished limit, we consider the asymptotic scalings of (A.1–A.5) with respect to small Froude and small Rossby number (a rapid rotating and strongly stably stratified flow), which gives

$$\begin{aligned} \mathrm{Ro}=\mathrm{Eu}^{-1}=\varepsilon , \quad \mathrm{Fr}_1=\mathrm{Ro}\dfrac{L}{L_{d_1}} = O{(\varepsilon )}, \quad {\mathrm{Fr}_2}=\mathrm{Ro}\dfrac{L}{L_{d_2}} = O{(\varepsilon )}, \quad \Gamma A^2={\mathrm{Fr}_1}^{-1}. \end{aligned}$$

(A.11)

Also, from Smith and Stechmann (2017) (equation (A7)), we have $\dfrac{c_p \theta _0}{L_v} = C_{cl}\mathrm{Ro}$. For simplicity, setting $C_{cl} = 1$, we have $\dfrac{c_p \theta _0}{L_v} = \varepsilon$.

With aforementioned asymptotic scaling and distinguished limit relationship, the non-dimensional model is displayed as:

$$\begin{aligned}&\dfrac{D_h\vec {u}_h}{Dt}+w\dfrac{\partial \vec {u}_h}{\partial z}+\varepsilon ^{-1}\vec {u}_{h}^{\bot }+\varepsilon ^{-1}{\nabla _h}\phi =0 \end{aligned}$$

(A.12)

$$\begin{aligned}&A^2\left( \dfrac{D_hw}{Dt}+w\dfrac{\partial w}{\partial z}\right) +\varepsilon ^{-1}\dfrac{\partial \phi }{\partial z}=\varepsilon ^{-1} \dfrac{L_{d_1}}{L} b \end{aligned}$$

(A.13)

$$\begin{aligned}&\nabla _h\vec {u}_h+\dfrac{\partial w}{\partial z}=0 \end{aligned}$$

(A.14)

$$\begin{aligned}&\dfrac{D_h\theta _e}{Dt}+w\dfrac{\partial \theta _e}{\partial z}+\varepsilon ^{-1}\dfrac{L_{d_1}}{L}w =0 \end{aligned}$$

(A.15)

$$\begin{aligned}&\dfrac{D_h q_t}{Dt}+ w\dfrac{\partial q_t}{\partial z}-\varepsilon ^{-1}\dfrac{L_{d_2}}{L}w- V_r \frac{\partial q_r}{\partial z} =0 \end{aligned}$$

(A.16)

Apart from the key non-dimensional parameter $\varepsilon ^{-1}$ shown above, $\varepsilon _1^{-1} = \varepsilon ^{-1}\dfrac{L_{d_1}}{L}, \varepsilon _2^{-1} = \varepsilon ^{-1}\dfrac{L_{d_2}}{L}$ will be defined, which are related to two Froude numbers. Furthermore, picking $L = L_{d_1} = L_{d_2}$ (implying $\varepsilon =\varepsilon _1=\varepsilon _2$) and $A = 1$ allows simple notation and gives:

$$\begin{aligned}&\dfrac{D_h\vec {u}_h}{Dt}+w\dfrac{\partial \vec {u}_h}{\partial z}+\varepsilon ^{-1}\vec {u}_{h}^{\bot }+\varepsilon ^{-1}{\nabla _h}\phi =0 \end{aligned}$$

(A.17)

$$\begin{aligned}&\dfrac{D_hw}{Dt}+w\dfrac{\partial w}{\partial z}+\varepsilon ^{-1}\dfrac{\partial \phi }{\partial z}=\varepsilon ^{-1} b \end{aligned}$$

(A.18)

$$\begin{aligned}&\nabla _h\vec {u}_h+\dfrac{\partial w}{\partial z}=0 \end{aligned}$$

(A.19)

$$\begin{aligned}&\dfrac{D_h\theta _e}{Dt}+w\dfrac{\partial \theta _e}{\partial z}+\varepsilon ^{-1}w =0 \end{aligned}$$

(A.20)

$$\begin{aligned}&\dfrac{D_h q_t}{Dt}+ w\dfrac{\partial q_t}{\partial z}-\varepsilon ^{-1}w- V_r \frac{\partial q_r}{\partial z} =0 \end{aligned}$$

(A.21)

Note that $V_r=0$, $V_r = 1$ or $V_r = \varepsilon ^{-1}$ is remained to be specified, since we consider different scenarios for rainfall (no rainfall, or normal speed $V_T = 0.1$ m/s or large speed $V_T = 1$ m/s). With the special choices above, where all O(1) constants were set equal to unity, we arrive at the advantageous situation where only one distinguished parameter $\varepsilon$ appears, to help simplify the notation.

Appendix B: Change of Variables in Different Environments

In this appendix, we will demonstrate a change of variables to a 4-dimensional state vector $\vec {v}^{\; \intercal } = (M,PV_e,W_1,W_2)$, which separates the zero-frequency variables $M,PV_e$ from the wave variables $W_1,W_2$, starting from the 5-d state vector $\vec {v}^{\; \intercal } = (u,v,w,\theta _e,q_t)$ (which is actually 4-dimensional due to the additional constraint of incompressibility, $u_x+v_y+w_z=0$ and the special horizontal mean flow case $u_m,v_m$ has been discussed in (2.30),(2.31)). Two cases will be considered: $V_r=0$ and $V_r\ne 0$.

1.1 $V_r = 0$ with Phase Changes

The starting point is the moist Boussinesq system with phase changes, which has a 5-d state vector $\vec {v}^{\; \intercal } = (u,v,w,\theta _e,q_t)$ with evolution equations

$$\begin{aligned}&\dfrac{D_h\vec {u}_h}{Dt}+w\dfrac{\partial \vec {u}_h}{\partial z}+\varepsilon ^{-1}\vec {u}_{h}^{\bot }+\varepsilon ^{-1}{\nabla _h}\phi =0 \end{aligned}$$

(B.1)

$$\begin{aligned}&\dfrac{D_hw}{Dt}+w\dfrac{\partial w}{\partial z}+\varepsilon ^{-1}\dfrac{\partial \phi }{\partial z}=\varepsilon _1^{-1}\left( b_\mathrm{u}H_\mathrm{u} + b_\mathrm{s}H_\mathrm{s} \right) \end{aligned}$$

(B.2)

$$\begin{aligned}&\nabla _h\cdot \vec {u}_h+\dfrac{\partial w}{\partial z}=0 \end{aligned}$$

(B.3)

$$\begin{aligned}&\dfrac{D_h\theta _e}{Dt}+w\dfrac{\partial \theta _e}{\partial z}+\varepsilon _1^{-1} w =0 \end{aligned}$$

(B.4)

$$\begin{aligned}&\dfrac{D_h q_t}{Dt}+ w\dfrac{\partial q_t}{\partial z}-\varepsilon _2^{-1}w =0 \end{aligned}$$

(B.5)

$$\begin{aligned}&\text {where } b_\mathrm{u} = \left[ \theta _e + \varepsilon q_t - q_t\right] , \quad b_\mathrm{s} = \left[ \theta _e - \varepsilon q_t\right] . \end{aligned}$$

(B.6)

Applying the curl operator ($\nabla _h\times$) on Eq. (B.1) leads to

$$\begin{aligned} \dfrac{\partial \xi }{\partial t}+\nabla _h\times \left( \vec {u}_h\cdot \nabla _h\vec {u}_h+w\dfrac{\partial \vec {u}_h}{\partial z}\right) +\varepsilon ^{-1}\delta =0, \end{aligned}$$

(B.7)

and applying the divergence operator ($\nabla _h\cdot$) on Eq. (B.1) leads to

$$\begin{aligned}&\Rightarrow \dfrac{\partial \delta }{\partial t}+\nabla _h\cdot \left( \vec {u}_h\cdot \nabla _h\vec {u}_h+w\dfrac{\partial \vec {u}_h}{\partial z}\right) -\varepsilon ^{-1}\xi +\varepsilon ^{-1}\nabla _h^2\phi =0, \end{aligned}$$

(B.8)

$$\begin{aligned}&\text {where }\delta =\nabla _h\times \vec {u}^\perp _h = u_x + v_y,\quad \xi =\nabla _h\times \vec {u}_h = v_x - u_y, \end{aligned}$$

(B.9)

$$\begin{aligned}&u_h^\perp =\left( \begin{array}{c} -v \\ u\end{array}\right) , \quad \nabla _h\cdot u_h^\perp =-v_x+u_y=-\xi . \end{aligned}$$

(B.10)

For simplicity, the usage of notation $NL_{\xi }$ denotes the nonlinear term in Eq. (B.7). Meanwhile, with the incompressibility condition given by Eq. (B.3), one may replace $\delta$ by $-w_z$, and thus (B.7) becomes

$$\begin{aligned} \dfrac{\partial \xi }{\partial t}+NL_{\xi }-\varepsilon ^{-1}w_z=0. \end{aligned}$$

(B.11)

By introducing a new variable M,

$$\begin{aligned} M=q_t+G_m\theta _e, \quad G_m = \dfrac{\varepsilon _2}{\varepsilon _1}, \end{aligned}$$

(B.12)

and adding (B.4) and (B.5) together, one finds

$$\begin{aligned} \dfrac{\partial M}{\partial t}+\vec {u}\cdot \nabla M=0, \quad \text{ or }\quad \dfrac{DM}{Dt}=0. \end{aligned}$$

(B.13)

By introducing a new variable $PV_e$,

$$\begin{aligned} PV_e=\xi +F\dfrac{\partial \theta _e}{\partial z}, \quad F = \dfrac{\varepsilon }{\varepsilon _1}, \end{aligned}$$

(B.14)

and applying the operator ($\partial _z$) on Eq. (B.4), one finds

$$\begin{aligned} \dfrac{\partial \left( \partial _z \theta _e \right) }{\partial t}+\partial _z\left( \vec {u}\cdot \nabla \theta _e\right) +\varepsilon _1^{-1}w_z=0. \end{aligned}$$

(B.15)

Adding (B.11) and (B.15) together leads to

$$\begin{aligned} \dfrac{\partial PV_e}{\partial t}+\partial _z\left( \vec {u}\cdot \nabla \theta _e\right) +NL_\xi =0. \end{aligned}$$

(B.16)

This completes the derivation of the $M, PV_e$ equations.

The next step is to present variables $W_1$ and $W_2$. Similarly one could substitute $-w_z$ for $\delta$ in Eq. (B.8) to arrive at

$$\begin{aligned}&\dfrac{\partial w_z}{\partial t}-\nabla _h\cdot \left( \vec {u}_h\cdot \nabla _h\vec {u}_h+w\dfrac{\partial \vec {u}_h}{\partial z}\right) +\varepsilon ^{-1}\xi =\varepsilon ^{-1}\nabla ^2_h \phi ,\nonumber \\&\quad \Rightarrow \dfrac{\partial w_{zz}}{\partial t}-\partial _z\nabla _h\cdot \left( \vec {u}_h\cdot \nabla _h\vec {u}_h+w\dfrac{\partial \vec {u}_h}{\partial z}\right) +\varepsilon ^{-1}\xi _z=\varepsilon ^{-1}\partial _z\nabla ^2_h \phi . \end{aligned}$$

(B.17)

By applying the operator ($\nabla _h^2$) on Eq. (B.2), one finds

$$\begin{aligned} \dfrac{\partial \nabla _h^2 w}{\partial t}+\nabla _h^2\left( \vec {u}\cdot \nabla w\right) + \varepsilon ^{-1} \nabla _h^2 \frac{\partial \phi }{\partial z} - \varepsilon _1^{-1}\nabla _h^2 \underbrace{\left( b_\mathrm{u}H_\mathrm{u}+b_\mathrm{s}H_\mathrm{s}\right) }_{b}=0. \end{aligned}$$

(B.18)

Combining (B.17) and (B.18) together will then cancel the pressure terms and yield:

$$\begin{aligned} \dfrac{\partial \nabla ^2 w}{\partial t}+\varepsilon ^{-1}\xi _{z}-\varepsilon _1^{-1}\nabla ^2_h b+\nabla ^2_h\left( \vec {u}\cdot \nabla w\right) -\partial _{z}\nabla _h\cdot \left( \vec {u}_h\cdot \nabla _h\vec {u}_h+w\dfrac{\partial \vec {u}_h}{\partial _z}\right) =0 . \end{aligned}$$

(B.19)

Based on the linear part of Eq. (B.19), we naturally generate two variables:

$$\begin{aligned} W_1=\nabla ^2 w, \end{aligned}$$

(B.20)

$$\begin{aligned} W_2= & {} \xi _{z}-F\nabla _h^2 b, \quad F = \dfrac{\varepsilon }{\varepsilon _1}. \end{aligned}$$

(B.21)

When $W_1,W_2$ are inserted into the linear part of (B.19); the result is

$$\begin{aligned} \dfrac{\partial W_1}{\partial t}+\varepsilon ^{-1}W_2+\nabla ^2_h\left( \vec {u}\cdot \nabla w\right) -\partial _{z}\nabla _h\cdot \left( \vec {u}_h\cdot \nabla _h\vec {u}_h+w\dfrac{\partial \vec {u}_h}{\partial _z}\right) =0 . \end{aligned}$$

(B.22)

In order to close the system, taking the time derivative of $W_2$ will lead to its evolution equation. Since the $W_2$ term contains b, we first focus attention on $\dfrac{\partial b}{\partial t}$ (note $b = H_\mathrm{u}b_\mathrm{u} + H_\mathrm{s}b_\mathrm{s}$). Recall the non-dimensional forms of $b_\mathrm{u}, b_\mathrm{s}$ in (B.6), which are just combinations of $\theta _e, q_t$. Hence $\dfrac{\partial b_\mathrm{u}}{\partial t}$, $\dfrac{\partial b_\mathrm{s}}{\partial t}$ easily yield following two equations for $b_\mathrm{u}$ and $b_\mathrm{s}$:

$$\begin{aligned}&\dfrac{\partial b_\mathrm{u}}{\partial t}+\vec {u}\cdot \nabla b_\mathrm{u}+\varepsilon ^{-1}_\mathrm{u}\cdot w=0, \end{aligned}$$

(B.23)

$$\begin{aligned}&\dfrac{\partial b_\mathrm{s}}{\partial t}+\vec {u}\cdot \nabla b_\mathrm{s}+\varepsilon ^{-1}_\mathrm{s}\cdot w=0, \end{aligned}$$

(B.24)

where $\varepsilon ^{-1}_\mathrm{u}$, $\varepsilon ^{-1}_\mathrm{s}$ are non-dimensional forms of the buoyancy frequencies and the corresponding dimensional forms are $N_\mathrm{u}^2$, $N_\mathrm{s}^2$ mentioned in (2.13). Thereby, together with (B.4) and (B.5), we can relate $\varepsilon _\mathrm{u}^{-1}$, $\varepsilon _\mathrm{s}^{-1}$ with $\varepsilon _1^{-1},$ $\varepsilon _2^{-1}$ as follows:

$$\begin{aligned} \varepsilon _\mathrm{u}^{-1} = \varepsilon _1^{-1} + \varepsilon _2^{-1} - \dfrac{\varepsilon }{\varepsilon _2}, \quad \varepsilon _\mathrm{s}^{-1} = \varepsilon _1^{-1} + \dfrac{\varepsilon }{\varepsilon _2}. \end{aligned}$$

(B.25)

Next, write down the time derivative for buoyancy,

$$\begin{aligned} \dfrac{\partial b}{\partial t}=\dfrac{\partial \left( b_\mathrm{u}H_\mathrm{u}+b_\mathrm{s}H_\mathrm{s}\right) }{\partial t}=\dfrac{\partial b_\mathrm{u}}{\partial t}H_\mathrm{u}+\dfrac{\partial b_\mathrm{s}}{\partial t}H_\mathrm{s}+\left( b_\mathrm{u}-b_\mathrm{s}\right) \partial _t H_\mathrm{u}. \end{aligned}$$

(B.26)

Note that $\left( b_\mathrm{u}-b_\mathrm{s}\right) \partial _t H_\mathrm{u}$ becomes zero because $\partial _t H_\mathrm{u}$ is a Dirac delta function at the phase interface, and $b_\mathrm{u}=b_\mathrm{s}$ at the phase interface. As a result, and using (B.23) and (B.24) described above, we find

$$\begin{aligned}&\dfrac{\partial b}{\partial t}=-H_\mathrm{u}\varepsilon ^{-1}_\mathrm{u} w-H_\mathrm{s}\varepsilon ^{-1}_\mathrm{s} w-H_\mathrm{u}\vec {u}\cdot \nabla b_\mathrm{u}-H_\mathrm{s}\vec {u}\cdot \nabla b_\mathrm{s},\nonumber \\&\quad \text{ or }\quad \dfrac{\partial b}{\partial t}+C_{\left( H\right) }w+H_\mathrm{u}\vec {u}\cdot \nabla b_\mathrm{u}+H_\mathrm{s}\vec {u}\cdot \nabla b_\mathrm{s}=0,\nonumber \\&\quad \text {where }\quad C_{\left( H\right) } = H_\mathrm{u}\varepsilon ^{-1}_\mathrm{u} + H_\mathrm{s}\varepsilon ^{-1}_\mathrm{s} = H_\mathrm{u}\left( \varepsilon _1^{-1} + \varepsilon _2^{-1} - \dfrac{\varepsilon }{\varepsilon _2} \right) + H_\mathrm{s}\left( \varepsilon _1^{-1} + \dfrac{\varepsilon }{\varepsilon _2} \right) . \end{aligned}$$

(B.27)

Note that $C_{(H)}$ as the coefficient of the linear part in (B.27) contains not only $O(\varepsilon ^{-1})$ terms but also O(1) terms. Pulling out the $\varepsilon ^{-1}$ part, one arrives at the following version of (B.27):

$$\begin{aligned} \dfrac{\partial b}{\partial t}+ \varepsilon ^{-1}\left[ H_\mathrm{u}\left( \dfrac{\varepsilon }{\varepsilon _1} + \dfrac{\varepsilon }{\varepsilon _2} - \dfrac{\varepsilon ^2}{\varepsilon _2} \right) + H_\mathrm{s}\left( \dfrac{\varepsilon }{\varepsilon _1} + \dfrac{\varepsilon ^2}{\varepsilon _2} \right) \right] w + H_\mathrm{u}\vec {u}\cdot \nabla b_\mathrm{u}+H_\mathrm{s}\vec {u}\cdot \nabla b_\mathrm{s}=0. \end{aligned}$$

(B.28)

Apply operator ($\nabla ^2_h$) on Eq. (B.27) leads to

$$\begin{aligned} \dfrac{\partial \nabla ^2_h b}{\partial t}+\nabla ^2_h \left( C_{\left( H\right) }w\right) +\nabla ^2_h\left( H_\mathrm{u}\vec {u}\cdot \nabla b_\mathrm{u}+H_\mathrm{s}\vec {u}\cdot \nabla b_\mathrm{s}\right) =0. \end{aligned}$$

(B.29)

With this information in hand, we can now return to $W_2$ itself. Taking the time derivative of variable $W_2 = \xi _{z}-F\nabla ^2_h b$ and combining the information from Eqs. (B.11) and (B.29), we find

$$\begin{aligned} \dfrac{\partial \left( \xi _{z}-F\nabla ^2_h b\right) }{\partial t}&=\varepsilon ^{-1}\partial _z^2\left( w\right) +F\nabla ^2_h\left( C_{\left( H\right) }w\right) \nonumber \\&\quad -\partial _z\left( NL_\xi \right) +F\nabla ^2_h\left( H_\mathrm{u}\vec {u}\cdot \nabla b_\mathrm{u}+H_\mathrm{s}\vec {u}\cdot \nabla b_\mathrm{s}\right) . \end{aligned}$$

(B.30)

With the replacement of $W_1 = \nabla ^2w$, $W_2 = \xi _{z} - F\nabla _h^2 b$ in linear part, one could update the previous equation as

$$\begin{aligned} \dfrac{\partial W_2}{\partial t}&=\varepsilon ^{-1}\partial _z^2\left( \nabla ^{-2}W_1\right) +F\nabla ^2_h\left( C_{\left( H\right) }\nabla ^{-2}W_1\right) \nonumber \\&\quad -\partial _z\left( NL_\xi \right) +F\nabla ^2_h\left( H_\mathrm{u}\vec {u}\cdot \nabla b_\mathrm{u}+H_\mathrm{s}\vec {u}\cdot \nabla b_\mathrm{s}\right) . \end{aligned}$$

(B.31)

This concludes the derivation for the case of $V_r=0$ with phase changes.

1.2 $V_r = 1$ within Purely Saturated Region

In now considering $V_r\ne 0$, in the following discussion, attention will be confined to purely saturated region, so that $H_\mathrm{u} = 0$ and $H_\mathrm{s} = 1$, without phase changes, but with the presence of rainfall in consideration. Consequently, the $q_t$ equation in (B.5) will have an extra $\dfrac{\partial q_t}{\partial z}$ term, as shown in

$$\begin{aligned} \dfrac{D_h q_t}{Dt}+ w\dfrac{\partial q_t}{\partial z}-\varepsilon _2^{-1}w = \dfrac{\partial q_t}{\partial z}. \end{aligned}$$

(B.32)

The above modification of the $q_t$ equation will go through in the derivations of the M equation and $W_2$ equation, which are constructed based on the variable $q_t$. By the definition in (B.12), one may rewrite (B.13) as

$$\begin{aligned} \dfrac{\partial M}{\partial t}+\vec {u}\cdot \nabla M=\dfrac{\partial q_t}{\partial z}, \quad \text{ or }\quad \dfrac{DM}{Dt}= \dfrac{\partial q_t}{\partial z}. \end{aligned}$$

(B.33)

Since in a purely saturated region we have $b_\mathrm{s} = \theta _e - \varepsilon q_t$, we observe that the impact of rainfall on the $W_2$ equation will emerge through (B.28). After inserting the rainfall term into the original (B.28), and restricting attention to the saturated, single-phase scenario, we find

$$\begin{aligned} \dfrac{\partial b_\mathrm{s}}{\partial t}+ \varepsilon ^{-1}\left( \dfrac{\varepsilon }{\varepsilon _1} + \dfrac{\varepsilon ^2}{\varepsilon _2} \right) w + \vec {u}\cdot \nabla b_\mathrm{s}=-\varepsilon \dfrac{\partial q_t}{\partial z}. \end{aligned}$$

(B.34)

Then we find the form of the $W_2$ equation in a purely saturated region, with rainfall impact:

$$\begin{aligned} \dfrac{\partial W_2}{\partial t}= & {} \varepsilon ^{-1}\partial _z^2\left( \nabla ^{-2}W_1\right) + \varepsilon ^{-1}F\nabla ^2_h\left( \left( \dfrac{\varepsilon }{\varepsilon _1} + \dfrac{\varepsilon ^2}{\varepsilon _2} \right) \nabla ^{-2}W_1\right) + \varepsilon F\nabla _h^2\left( \dfrac{\partial q_t}{\partial z}\right) \nonumber \\&-\partial _z\left( NL_\xi \right) +F\nabla ^2_h\left( \vec {u}\cdot \nabla b_\mathrm{s}\right) . \end{aligned}$$

(B.35)

Though the $\dfrac{\partial q_t}{\partial z}$ term has been introduced into this equation, it arises at order $O(\varepsilon )$, which will not explicitly show up in the leading orders of behavior of $W_2$ related to $\mathscr {L}_*,$ $\mathscr {L}_0$. Nevertheless, the rainfall term still impacts the M evolution at leading order, as shown in (B.33).

1.3 $V_r = O(\varepsilon ^{-1})$ Within Purely Saturated Region

A similar argument can be implemented here with $V_r = O(\varepsilon ^{-1})$. The corresponding adjusted M, $W_2$ equations are given by

$$\begin{aligned} \dfrac{\partial M}{\partial t}+\vec {u}\cdot \nabla M= & {} \varepsilon ^{-1}\dfrac{\partial q_t}{\partial z}, \end{aligned}$$

(B.36)

$$\begin{aligned} \dfrac{\partial W_2}{\partial t}= & {} \varepsilon ^{-1}\partial _z^2\left( \nabla ^{-2}W_1\right) + \varepsilon ^{-1}F\nabla ^2_h\left( \left( \dfrac{\varepsilon }{\varepsilon _1} + \dfrac{\varepsilon ^2}{\varepsilon _2} \right) \nabla ^{-2}W_1\right) + F\nabla _h^2\left( \dfrac{\partial q_t}{\partial z}\right) \nonumber \\&-\partial _z\left( NL_\xi \right) +F\nabla ^2_h\left( \vec {u}\cdot \nabla b_\mathrm{s}\right) . \end{aligned}$$

(B.37)

Appendix C: Inverse Change of Variables to Recover $(u,v,w,\theta _e,q_t)$

In this appendix, we show how to recover the variables $(u,v,w,\theta _e,q_t)$, given the variables $(PV_e,M,W_1,W_2,u_m,v_m)$. In a sense, this is a type of PV inversion, although also including M and waves $W_1,W_2,u_m,v_m$.

The definition of $b_\mathrm{u},$ and $W_2$ give

$$\begin{aligned} b_\mathrm{u}= & {} \left( \theta _e+\varepsilon q_t-q_t\right) , \quad b_\mathrm{s}=\left( \theta _e - \varepsilon q_t\right) , \end{aligned}$$

(C.1)

$$\begin{aligned} W_2= & {} \xi _{z}-F\nabla _h^2\left( H_\mathrm{u}b_\mathrm{u}+H_\mathrm{s}b_\mathrm{s} \right) , \end{aligned}$$

(C.2)

and when $b_\mathrm{u}$, $b_\mathrm{s}$ are inserted into (C.2), the $W_2$ equation in terms of $\theta _e, q_t$ yields

$$\begin{aligned} W_2=\xi _{z}-F\nabla _h^2\left( H_\mathrm{u}\left( \theta _e+\varepsilon q_t-q_t\right) +H_\mathrm{s}\left( \theta _e-\varepsilon q_t\right) \right) , \end{aligned}$$

(C.3)

or

$$\begin{aligned} \xi _{z}-W_2=F\nabla ^2_h\left( \theta _e-H_\mathrm{u}q_t+\varepsilon q_t \right) . \end{aligned}$$

(C.4)

Through neglecting $\varepsilon q$, we only put O(1) balanced terms into consideration, implying leading order inversion formula in the end. Replacing $q_t$ by $M- G_m \theta _e$ (for simplicity setting $G_m = 1$, $F = 1$) shows

$$\begin{aligned}&\xi _{z}-W_2=\nabla ^2_h\left( \theta _e-H_\mathrm{u}\left( M-\theta _e\right) \right) , \end{aligned}$$

(C.5)

$$\begin{aligned}&{\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) =\left( 1+H_\mathrm{u}\right) \theta _e-H_\mathrm{u}M , \end{aligned}$$

(C.6)

$$\begin{aligned}&{\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) +H_\mathrm{u}M=\left( 1+H_\mathrm{u}\right) \theta _e , \end{aligned}$$

(C.7)

$$\begin{aligned}&\theta _e = \dfrac{1}{2}H_\mathrm{u}\left[ {\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) + M\right] + H_\mathrm{s}\left[ {\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) \right] . \end{aligned}$$

(C.8)

The aforementioned straightforward work only depends on definitions of buoyancy $b_\mathrm{u}$, $b_\mathrm{s}$, $W_2$, and M, which simply express $\theta _e$ in terms of $M, \xi , W_2$. The next goal is to write down the inversion of $\xi$ with respect to M, $PV_e$ and $W_2$.

To find the inversion PDE, we first apply operator ($\partial _z$) to (C.8), and we see that $\partial _z (\theta _e)$ equals

$$\begin{aligned} \dfrac{\partial }{\partial z} \left\{ \dfrac{1}{2}H_\mathrm{u}[{\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) + M] + H_\mathrm{s}\left[ {\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) \right] \right\} . \end{aligned}$$

(C.9)

Now recall the definition of $PV_e = \xi + F\dfrac{\partial \theta _e}{\partial z}$ (for simplicity setting $F = 1$), and notice that $\dfrac{\partial \theta _e}{\partial z}$ could be replaced by (C.9) to yield

$$\begin{aligned} \xi + \dfrac{\partial }{\partial z} \left\{ \dfrac{1}{2}H_\mathrm{u}\left[ {\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) + M\right] + H_\mathrm{s}\left[ {\nabla ^{-2}_h}\left( \xi _{z}-W_2\right) \right] \right\} = PV_e. \end{aligned}$$

(C.10)

If a streamfunction $\psi = (\nabla _h^{-2})\xi$ is introduced, which also implies $\xi = (\nabla _h^2)\psi$, $(\nabla _h^{-2})\xi _{z} = \psi _{z}$, one can rewrite (C.10) as

$$\begin{aligned} \nabla _h^2\psi + \dfrac{\partial }{\partial z} \left\{ \dfrac{1}{2}H_\mathrm{u}\left[ \partial _z\psi -{\nabla ^{-2}_h}W_2 + M\right] + H_\mathrm{s}\left[ \partial _z\psi -{\nabla ^{-2}_h}W_2\right] \right\} = PV_e. \end{aligned}$$

(C.11)

This is an elliptic PDE for the streamfunction $\psi$, given $PV_e$, M, and $W_2$. It is an extension of PV-and-M inversion (Wetzel et al. 2019, 2020) and now includes the influence of waves via $W_2$.

An important point is that the PDE (C.11) illustrates how $\psi$ is influenced by fast waves in two ways. First, as mentioned above, the presence of $W_2$ is one clear influence of waves. Second, recall that the Heaviside functions $H_\mathrm{u},H_\mathrm{s}$ also introduce $t, \tau$ dependence. In fact, even if one considers the recovery of $\psi _{(M,PV_e)}$ (by considering a case of recovery from given $M,PV_e$ with setting $W_1 = 0, W_2 = 0$), the $\tau$-dependence of $H_\mathrm{u},H_\mathrm{s}$ will introduce a fast $\tau$-dependence to $\psi _{(M,PV_e)}$, even though M and $PV_e$ themselves have no $\tau$-dependence. It shows how waves can influence $\psi _{(M,PV_e)}$ via phase changes.

Solving the elliptic PDE in (C.11) provides $\psi$ in terms of $(M, PV_e, W_2)$. Accordingly, knowledge of $\psi$ helps us to derive the inversion formulas for the velocity field $\vec {u}^{\; \intercal } = (u,v,w)$, which could be determined from $\psi$, $W_1$ and finally be expressed as $(M, PV_e, W_1, W_2,u_m,v_m)$ only.

Similarly, the definition of $W_1 = {\nabla ^{2}} w$ demonstrates

$$\begin{aligned} w={\nabla ^{-2}} W_1. \end{aligned}$$

(C.12)

With the incompressibility condition

$$\begin{aligned} u_x+v_y=-w_z=-{\left( \partial _z\nabla ^{-2}\right) } W_1, \end{aligned}$$

(C.13)

and the definition of $\xi = v_x - u_y$, we arrive at

$$\begin{aligned} v_{xx}+v_{yy}= & {} \xi _x-{\left( \partial _y \partial _z \nabla ^{-2}\right) } W_1, \end{aligned}$$

(C.14)

$$\begin{aligned} u_{xx}+u_{yy}= & {} -\xi _y-{\left( \partial _x \partial _z \nabla ^{-2}\right) } W_1. \end{aligned}$$

(C.15)

The results of u, v are expressed as

$$\begin{aligned} v= & {} {\left( \nabla _h^{-2}\right) }\left( \xi _x-{\left( \partial _y \partial _z\nabla ^{-2}\right) } W_1\right) , \end{aligned}$$

(C.16)

$$\begin{aligned} u= & {} {\left( \nabla _h^{-2}\right) }\left( -\xi _y-{\left( \partial _x\partial _z\nabla ^{-2}\right) } W_1\right) . \end{aligned}$$

(C.17)

As a more physically revealing form, one can rewrite (C.16)–(C.17) as

$$\begin{aligned} v - v_m= & {} \partial _x\psi -\partial _y \partial _z\left( \nabla _h^{-2} \nabla ^{-2} W_1\right) , \end{aligned}$$

(C.18)

$$\begin{aligned} u - u_m= & {} -\partial _y\psi -\partial _x \partial _z\left( \nabla _h^{-2} \nabla ^{-2} W_1\right) , \end{aligned}$$

(C.19)

where $u_m$, $v_m$ are mean velocities and subscript m denotes the horizontal average. (C.18)–(C.19) displays the contributions from the streamfunction $\psi$, mean velocities and from the velocity potential $-\nabla _h^{-2} \nabla ^{-2} W_1$ that is due to waves. Since $\psi$ could be found from (C.11) and written in terms of $(M,PV_e,W_2)$, we see that the velocity field $\vec {u}^{\; \intercal } = (u,v,w)$ could be obtained through inverting state vector $\vec {v}^{\; \intercal } = (M, PV_e, W_1, W_2,u_m,v_m)$.

The following contents offer a special inversion formula for the single phase case (purely saturated region with $H_\mathrm{u} =0, H_\mathrm{s} = 1$), under no presence of wave ($W_1 = 0, W_2 = 0, u_m = 0, v_m = 0$), which supports conclusions demonstrated on Sect. 5. In a purely saturated region ($H_\mathrm{s} = 1, H_\mathrm{u} = 0$), (C.10) becomes

$$\begin{aligned} \xi + {\partial _z\left( \nabla ^{-2}_h\right) }\left( \xi _{z}-W_2\right) = PV_e. \end{aligned}$$

(C.20)

The remaining work is to introduce the streamfunction $\psi = (\nabla _h^{-2})\xi$, which implies $\xi = (\nabla _h^2)\psi$, $(\nabla _h^{-2})\xi _{zz} = \psi _{zz}$ in (C.20). Without considering the impact of waves, setting $W_2 = 0$ in (C.20) leads to

$$\begin{aligned} \nabla ^2 \psi = PV_e. \end{aligned}$$

(C.21)

Then $\vec {u}_{(M,PV_e)}$, as the slow part velocity field, coming from (C.12, C.18, C.19) with $W_1=0$, $u_m = v_m =0$, and $\xi = (\nabla _h^2)\psi$, is given by

$$\begin{aligned} u_{(M,PV_e)} = -\psi _y, \quad v_{(M,PV_e)} = \psi _x, \quad w_{(M,PV_e)} = 0. \end{aligned}$$

(C.22)

The slow thermaldynamic variable $\theta _{e(M,PV_e)}$, with contributions from $M, PV_e$ slow components only, is derived through (C.8), with $H_\mathrm{u} = 0, H_\mathrm{s} = 1, W_2 = 0, \xi = (\nabla _h^2)\psi$:

$$\begin{aligned} \theta _{e(M,PV_e)} = \psi _z. \end{aligned}$$

(C.23)

Finally, the definition of $M = \theta _e + q_t$ directly expresses slow variable $q_{t(M,PV_e)}$ as

$$\begin{aligned} q_{t\left( M,PV_e\right) } = M - \psi _z. \end{aligned}$$

(C.24)

Appendix D: Fourier decomposition of $\mathscr {L}_*$

Two different scenarios will be presented corresponding to the purely saturated region with two different rainfall speeds $V_r = 1$ and $V_r=\varepsilon ^{-1}$ (these two cases may be generalized to $V_r = O(1)$ and $V_r=O(\varepsilon ^{-1})$, respectively). The Fourier analysis in following “Appendix D”, “Appendix E” will answer the main question: Will the slow component $\bar{v}_\mathrm{slow}(t, \vec {x})$ evolve independently from the fast component, as in (3.7)–(3.8), even in the presence of precipitation $V_r$? Or will precipitation $V_r$ introduce an influence of the fast waves on the evolution of the slow component? Eventually, exactly analogous equations for suitably-defined potential vorticity variables displayed in “Appendix E” clarifies that independence between slow and fast components. In other words, there is no impact from rainfall on slow modes evolution.

Working through the Fourier decomposition of $\mathscr {L}_*$, we use dimensional variables in order to make explicit the appearance of the dimensional frequencies $N_1, N_2$ described in (2.12), Coriolis parameter f and dimensional rainfall speed $V_T$, helping to elucidate the dominant physics and to make contact with previous literature, e.g. Phillips (1968), Lelong and Riley (1991), Bartello (1995), Embid and Majda (1998), Majda and Embid (1998) and Remmel and Smith (2009). Based on the dimensional system (1a)–(1d) of Smith and Stechmann (2017) (see also 17(b) in Smith and Stechmann (2017) with $q_{vs}(z) = 0$), it is convenient to use rescaled variables

$$\begin{aligned} \theta _e' = \dfrac{g}{\theta _0}\frac{\theta _e}{N_1}, \text { and } q_t' = \dfrac{gL_v}{\theta _0c_p}\frac{q_t}{N_2}. \end{aligned}$$

(D.1)

Then the modified dynamic system in dimensional form will be given:

$$\begin{aligned}&\dfrac{D\vec {u}}{Dt}+f\hat{z}\times \vec {u}= -\nabla \dfrac{\phi }{\rho _0}+\hat{z} \left( N_1\theta _e' - \dfrac{\theta _0c_p}{L_v}N_2q_t'\right) \end{aligned}$$

(D.2)

$$\begin{aligned}&\nabla \cdot \vec {u}=0 \end{aligned}$$

(D.3)

$$\begin{aligned}&\dfrac{D\theta _e'}{Dt} + N_1w=0 \end{aligned}$$

(D.4)

$$\begin{aligned}&\dfrac{Dq_t'}{Dt} - N_2w - V_T\dfrac{\partial q_t'}{\partial z}=0 \end{aligned}$$

(D.5)

With the assumption of periodic boundary conditions in the spatial domain, we try to seek dispersion relation, writing special eigenfunction wave solution as

$$\begin{aligned} \vec {v} = e^{\left( i\vec {k}\cdot \vec {x} - i\sigma (\vec {k})t\right) }\vec {\phi }, \end{aligned}$$

(D.6)

where $\vec {k}$ is the wave number, $\sigma (\vec {k})$ is the eigenfrequencies, $\vec {\phi }$ is the eigenvector, and $\vec {v}$ should satisfy the incompressibility condition. Similarly, as described in Sect. 2.1, after non-dimensional process, one could fill the system (D.2 – D.5) in the abstract formulation (2.2) to construct concrete $\mathscr {L_*}$ and $\mathscr {L}_0$ as follows. (Note that the pressure term is rewritten using the expression $\Delta \phi =-\varepsilon \nabla \cdot (\vec {u}\cdot \nabla \vec {u})+\partial \theta _e/\partial z-\varepsilon \partial q_t/\partial z+\xi$.)

${\varvec{V}}_{\varvec{r}} = {\varvec{1:}}$

$$\begin{aligned} \mathscr {L}_*(\vec {v})= & {} \left( \begin{array}{ccccc} -\partial _x \Delta ^{-1}\partial _y &{} -1 + \partial x \Delta ^{-1}\partial _x &{} 0 &{}\partial _x \Delta ^{-1}\partial _z &{} 0\\ 1 - \partial _y \Delta ^{-1}\partial _y &{} \partial _y \Delta ^{-1}\partial _x &{} 0 &{}\partial _y \Delta ^{-1}\partial _z &{} 0\\ - \partial _z \Delta ^{-1} \partial _y &{} \partial _z \Delta ^{-1} \partial _x &{} 0 &{} \partial _z \Delta ^{-1}\partial _z - 1 &{} 0\\ 0&{}0&{}1&{}0&{}0\\ 0&{}0&{}-1&{}0&{}0 \end{array}\right) \left( \begin{array}{cc} &{}u\\ &{}v\\ &{}w\\ &{}{\theta _e}\\ &{}q_t\\ \end{array}\right) \end{aligned}$$

(D.7)

$$\begin{aligned} \mathscr {L}_0(\vec {v})= & {} \left( \begin{array}{ccccc} 0&{}0&{}0&{}0&{}-\partial _x \Delta ^{-1} \partial _z\\ 0&{}0&{}0&{}0&{}-\partial _y \Delta ^{-1} \partial _z\\ 0&{}0&{}0&{}0&{}1 -\partial _z \Delta ^{-1} \partial _z\\ 0&{}0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0&{}-\partial _z \end{array}\right) \left( \begin{array}{cc} &{}u\\ &{}v\\ &{}w\\ &{}{\theta _e}\\ &{}q_t\\ \end{array}\right) \end{aligned}$$

(D.8)

${\varvec{V}}_{\varvec{r}} = {\varvec{\varepsilon }}^{{\varvec{-1}}}:$

$$\begin{aligned} \mathscr {L}_*(\vec {u})= & {} \left( \begin{array}{ccccc} -\partial x \Delta ^{-1}\partial y &{} -1 + \partial x \Delta ^{-1}\partial x &{} 0 &{}\partial x \Delta ^{-1}\partial z &{} 0\\ 1 - \partial y \Delta ^{-1}\partial y &{} \partial y \Delta ^{-1}\partial x &{} 0 &{}\partial y \Delta ^{-1}\partial z &{} 0\\ - \partial z \Delta ^{-1} \partial y &{} \partial z \Delta ^{-1} \partial x &{} 0 &{} \partial z \Delta ^{-1}\partial z - 1 &{} 0\\ 0&{}0&{}1&{}0&{}0\\ 0&{}0&{}-1&{}0&{}- \partial z \end{array}\right) \left( \begin{array}{cc} &{}u\\ &{}v\\ &{}w\\ &{}{\theta _e}\\ &{}q_t\\ \end{array}\right) \end{aligned}$$

(D.9)

$$\begin{aligned} \mathscr {L}_0(\vec {u})= & {} \left( \begin{array}{ccccc} 0&{}0&{}0&{}0&{}-\partial x \Delta ^{-1} \partial z\\ 0&{}0&{}0&{}0&{}-\partial y \Delta ^{-1} \partial z\\ 0&{}0&{}0&{}0&{}1 -\partial z \Delta ^{-1} \partial z\\ 0&{}0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0&{}0 \end{array}\right) \left( \begin{array}{cc} &{}u\\ &{}v\\ &{}w\\ &{}{\theta _e}\\ &{}q_t\\ \end{array}\right) \end{aligned}$$

(D.10)

The implementation of Fourier transform $\mathscr {F}: (x,y,z,t) \rightarrow (k,l,m,\sigma )$ on the $\varepsilon ^{-1}$ balance part of abstract equation (2.2), which is $\dfrac{\partial \vec {v}}{\partial t}+\varepsilon ^{-1}\mathscr {L}_*(\vec {v})=0$, will directly give the following matrix equation

$$\begin{aligned} -i\sigma \vec {\phi } = -\tilde{A}_*\vec {\phi }. \end{aligned}$$

(D.11)

The associated matrix $\tilde{A}_*$, $\tilde{A}_0$ with respect to the dimensional form of $\varepsilon ^{-1}\mathscr {L}_*$, $\mathscr {L}_0$ are displayed below. (Note that $A_* = -{\mid \vec {k}\mid }^2 \tilde{A}_*$, $A_0 = -{\mid \vec {k}\mid }^2 \tilde{A}_0$.)

${\varvec{V}}_{\varvec{r}} = {\varvec{1:}}$

$$\begin{aligned} A_*= & {} \left( \begin{matrix} klf&{}({\mid \vec {k}\mid }^2 - k^2)f&{}0&{}-kmN_1&{}0\\ (-{\mid \vec {k}\mid }^2 + l^2)f&{}-klf&{}0&{}-lmN_1&{}0\\ lmf&{}-kmf&{}0&{}{k_h}^2 N_1&{}0\\ 0&{}0&{}-{\mid \vec {k}\mid }^2 N_1&{}0&{}0\\ 0&{}0&{}{\mid \vec {k} \mid }^2 N_2&{}0&{}0 \end{matrix}\right) \end{aligned}$$

(D.12)

$$\begin{aligned} A_0= & {} \left( \begin{matrix} 0&{}0&{}0&{}0&{}km\frac{\theta _0c_p}{L_v}N_2\\ 0&{}0&{}0&{}0&{}lm\frac{\theta _0c_p}{L_v}N_2\\ 0&{}0&{}0&{}0&{}-k^2_h\frac{\theta _0c_p}{L_v}N_2\\ 0&{}0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0&{}im{\mid \vec {k} \mid }^2 V_T \end{matrix}\right) \end{aligned}$$

(D.13)

${\varvec{V}}_{\varvec{r}} = {\varvec{\varepsilon }}^{{\varvec{-1}}}:$

$$\begin{aligned} A_*= & {} \left( \begin{matrix} klf&{}\left( {\mid \vec {k}\mid }^2 - k^2\right) f&{}0&{}-kmN_1&{}0\\ \left( -{\mid \vec {k}\mid }^2 + l^2\right) f&{}-klf&{}0&{}-lmN_1&{}0\\ lmf&{}-kmf&{}0&{}{k_h}^2 N_1&{}0\\ 0&{}0&{}-{\mid \vec {k}\mid }^2 N_1&{}0&{}0\\ 0&{}0&{}{\mid \vec {k} \mid }^2 N_2&{}0&{}im{\mid \vec {k} \mid }^2 V_T \end{matrix}\right) \end{aligned}$$

(D.14)

$$\begin{aligned} A_0= & {} \left( \begin{matrix} 0&{}0&{}0&{}0&{}km\frac{\theta _0c_p}{L_v}N_2\\ 0&{}0&{}0&{}0&{}lm\frac{\theta _0c_p}{L_v}N_2\\ 0&{}0&{}0&{}0&{}-k^2_h\frac{\theta _0c_p}{L_v}N_2\\ 0&{}0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0&{}0 \end{matrix}\right) \end{aligned}$$

(D.15)

By the incompressibility condition, notice that

$$\begin{aligned} k\hat{u}+l\hat{v}+m\hat{w}=0, \Rightarrow kl\hat{u}+l^2\hat{v}+lm\hat{w}=0, k^2\hat{u}+kl\hat{v}+km\hat{w}=0, \end{aligned}$$

(D.16)

and simple algebra presents

$$\begin{aligned} -{\mid {k}\mid }^2N_1 \hat{w} = -m^2N_1 \hat{w} - k_h^2N_1\hat{w} = km N_1 \hat{u} + lm N_1 \hat{v} - k_h^2N_1\hat{w}. \end{aligned}$$

(D.17)

Similarly, ${\mid k \mid }^2 N_2\hat{w}$ could be expressed as

$$\begin{aligned} {\mid k \mid }^2 N_2\hat{w} = -kmN_2\hat{u}-lm N_2\hat{v}+{k_h}^2 N_2\hat{w}. \end{aligned}$$

(D.18)

Complete the symmetrization for the $4 \times 4$ sub-matrix of $A_*$, giving analogous structure (see (D.19)) with previous literature (Embid and Majda 1998; Majda and Embid 1998; Remmel and Smith 2009), so as the corresponding eigenvectors. Since the last column entries of $A_*$ are different from dry case, which breaks the symmetrizing process for full matrix. In an abuse of notation, we use $\phi$ to replace $\vec {\phi }$ in following content, if there is no misunderstanding and contradiction.

For $V_r = 1$ case, new matrix $A_{s*}$ and associated eigenvalues, eigenvectors are given as:

$$\begin{aligned} {A}_{s*}= & {} \left( \begin{array}{ccccc} 0&{}m^2f&{}-lmf&{}-kmN_1&{}0\\ -m^2f&{}0&{}kmf&{}-lmN_1&{}0\\ lmf&{}-kmf&{}0&{}k^2_hN_1&{}0\\ kmN_1&{}lmN_1&{}-k^2_hN_1&{}0&{}0\\ 0&{}0&{}N_2\vert \vec {k}\vert ^2&{}0&{}0 \end{array}\right) \end{aligned}$$

(D.19)

$$\begin{aligned} \sigma= & {} 0 \text { (triple) } \quad \sigma ^2=\frac{N_1^2k^2_h+f^2m^2}{{\vert k\vert }^2}\quad \left( \sigma = {\vert \sigma ^{\pm }\vert }\right) \end{aligned}$$

(D.20)

$$\begin{aligned} {\phi }^0= & {} \dfrac{1}{\sqrt{N_1^2k^2_h+f^2m^2}}\left( \begin{array}{c} -N_1l\\ N_1k\\ 0\\ mf\\ 0 \end{array}\right) \quad \phi ^q=\left( \begin{array}{c} 0\\ 0\\ 0\\ 0\\ 1 \end{array}\right) \quad \phi ^{\pm }=\left( \begin{array}{c} \frac{m}{k_h}(\sigma k \pm ilf)\\ \frac{m}{k_h}(\sigma l \mp ikf)\\ -\sigma k_h\\ \pm iN_1k_h\\ \mp iN_2k_h\\ \end{array}\right) \end{aligned}$$

(D.21)

A special case must be considered, which is $k_h = 0$:

$$\begin{aligned} \sigma= & {} 0 \text { (triple) } \quad \sigma ^2= f^2\quad \left( \sigma = {\vert \sigma ^{\pm }\vert }\right) \end{aligned}$$

(D.22)

$$\begin{aligned} {\phi }^0= & {} \left( \begin{array}{c} 0\\ 0\\ 0\\ 1\\ 0 \end{array}\right) \quad \phi ^q=\left( \begin{array}{c} 0\\ 0\\ 0\\ 0\\ 1 \end{array}\right) \quad \phi ^{\pm }=\left( \begin{array}{c} \frac{1+i}{2}\\ \frac{1-i}{2}\\ 0\\ 0\\ 0\\ \end{array}\right) \end{aligned}$$

(D.23)

The first two eigenvectors have 0 eigenfrequencies, called slow modes, while fast modes represent the rest of two vectors with nonzero frequencies. Meanwhile, one eigenvector corresponding to 0 eigenvalue has been abandoned, since it violates the incompressibility condition. Orthogonality of the associated eigenvectors is not guaranteed. Nevertheless, one may process to analyse one of the slow modes ($\phi ^{0}$ mode also known as $PV_e$ mode) by projecting (3.6) into $\phi ^{0}$ mode in Fourier space, since $\phi ^{0}$ is perpendicular to the rest of three modes $\phi ^{q}, \phi ^{+}, \phi ^{-}$.

For $V_r = \varepsilon ^{-1}$ case, with similar argument we simply demonstrate the results of matrix $A_*$, eigenvalues and eigenvectors as follows:

$$\begin{aligned} A_*= & {} \left( \begin{array}{ccccc} klf&{}({\mid \vec {k}\mid }^2 -k^2 )f&{}0&{}-kmN_1&{}0\\ (-{\mid \vec {k}\mid }^2 + l^2)f&{}-klf&{}0&{}-lmN_1&{}0\\ lmf&{}-kmf&{}0&{}{k_h}^2 N_1&{}0\\ 0&{}0&{}-{\mid \vec {k}\mid }^2 N_1 &{}0&{}0\\ 0&{}0&{}{\mid \vec {k} \mid }^2 N_2 &{}0&{}im{\mid \vec {k} \mid }^2 V_T \end{array}\right) \end{aligned}$$

(D.24)

$$\begin{aligned} \sigma ^0= & {} 0 \text { (double) } \quad \sigma ^q=-mV_T \quad \sigma ^2=\frac{N_1^2 k^2_h+f^2m^2}{{\vert \vec {k}\vert }^2}\quad \left( \sigma = {\vert \sigma ^{\pm }\vert }\right) \end{aligned}$$

(D.25)

$$\begin{aligned} \phi ^0= & {} \dfrac{1}{\sqrt{N_1^2 k^2_h+f^2m^2}}\left( \begin{array}{c} -N_1l\\ N_1k\\ 0\\ mf\\ 0 \end{array}\right) \quad \phi ^q=\left( \begin{array}{c} 0\\ 0\\ 0\\ 0\\ 1 \end{array}\right) \quad \phi ^\pm =\left( \begin{array}{c} \frac{m}{k_h}\left( \sigma k \pm ilf\right) \\ \frac{m}{k_h}\left( \sigma l \mp ikf\right) \\ -\sigma k_h\\ \pm iN_1k_h\\ -\frac{iN_2k_h\sigma }{mV_T \pm \sigma }\\ \end{array}\right) \end{aligned}$$

(D.26)

And the special case $k_h = 0$ yields

$$\begin{aligned} \sigma= & {} 0 \text { (double) } \quad \sigma ^q=-mV_T \quad \sigma ^2= f^2\quad \left( \sigma = {\vert \sigma ^{\pm }\vert }\right) \end{aligned}$$

(D.27)

$$\begin{aligned} {\phi }^0= & {} \left( \begin{array}{c} 0\\ 0\\ 0\\ 1\\ 0 \end{array}\right) \quad \phi ^q=\left( \begin{array}{c} 0\\ 0\\ 0\\ 0\\ 1 \end{array}\right) \quad \phi ^{\pm }=\left( \begin{array}{c} \frac{1+i}{2}\\ \frac{1-i}{2}\\ 0\\ 0\\ 0\\ \end{array}\right) . \end{aligned}$$

(D.28)

It’s worth to remind reader here, under $V_r = \varepsilon ^{-1}$ and $m \ne 0$ circumstance, there is only one slow mode $\phi ^0$ since $\phi ^q$ is no longer to be slow due to the nonzero eigenvalue $\sigma ^q$.

Appendix E: Analysis of Resonant Interaction for Slow Dynamics

Based on the well-constructed eigenvectors described above, we start to build the concrete form of the average equation (3.6) in Fourier space. In the end, through the analysis of resonant triad interactions arising from bi-linear operator ($\mathscr {B}$) one could verify whether the decoupling property between slow and fast modes is still valid in the limit $\varepsilon \rightarrow 0$ under the presence water ($q_t$) and rainfall ($V_T$).

Initial condition $\bar{v}(\vec {x},t)$ in (3.2) is written in terms of the aforementioned eigenvectors $\phi ^{(\alpha )}_{(\vec {k})}$ (D.21) or (D.26) together with amplitude function $a^{(\alpha )}_{(\vec {k})}(t)$,

$$\begin{aligned} \bar{v}(\vec {x},t)=\sum _{\vec {k}\in \mathbb {Z}^3}\sum _{\alpha \in \mathscr {A}} e^{i\vec {k}\cdot \vec {x}} a^{(\alpha )}_{(\vec {k})}(t)\phi ^{(\alpha )}_{(\vec {k})}, \quad { \mathscr {A} = \{ 0, q, +, - \} }. \end{aligned}$$

(E.1)

Plugging (E.1) into $\mathscr {B}$, thus the bi-linear term could be represented explicitly

$$\begin{aligned}&\mathscr {B}(e^{-s\mathscr {L}_*}\bar{v},e^{-s\mathscr {L}_*}\bar{v}) =\nonumber \\&\quad =\sum _{\vec {k}\in \mathbb {Z}^3}\sum _{\alpha \in \mathscr {A}}\left\{ \sum _{(\vec {k}'+\vec {k}''=\vec {k})}\sum _{(\alpha ', \alpha '' \in \mathscr {A})}e^{i(\vec {k}\cdot \vec {x}-s(\sigma ^{(\alpha ')}_{(\vec {k}')}+\sigma ^{(\alpha '')}_{(\vec {k}'')}))} B^{(\alpha ', \alpha '', \alpha )}_{(\vec {k}', \vec {k}'', \vec {k})} a^{(\alpha ')}_{(\vec {k}')}(t) a^{({\alpha ''})}_{(\vec {k}'')}(t) \right\} \phi ^{(\alpha )}_{(\vec {k})}, \end{aligned}$$

(E.2)

where the coefficient B arrives to be

$$\begin{aligned} B^{(\alpha ',\alpha '',\alpha )}_{(\vec {k}',\vec {k}'',\vec {k})}=\dfrac{i}{2}\left[ (\vec {u}^{(\alpha ')}_{(\vec {k}')}\cdot \vec {k}'')(\vec {\phi }^{(\alpha '')}_{(\vec {k}'')}\cdot \vec {\phi }^{(\alpha )}_{(\vec {k})})+(\vec {u}^{(\alpha '')}_{(\vec {k}'')}\cdot \vec {k}')(\vec {\phi }^{(\alpha ')}_{(\vec {k}')}\cdot \vec {\phi }^{(\alpha )}_{(\vec {k})})\right] . \end{aligned}$$

(E.3)

Hence the quadratic contribution due to bi-linear operator $\mathscr {B}$ in the abstract averaging equation (3.6) is given as

$$\begin{aligned}&\lim \limits _{\tau \rightarrow \infty }\dfrac{1}{\tau }\int _0^{\tau } e^{s\mathscr {L}_*}\left( \mathscr {B}(e^{-s\mathscr {L}_*}\bar{v},e^{-s\mathscr {L}_F}\bar{v})\right) \mathrm{d}s = \nonumber \\&\quad =\lim \limits _{\tau \rightarrow \infty }\dfrac{1}{\tau }\int _0^{\tau } \sum _{\vec {k}\in \mathbb {Z}^3}\sum _{\alpha \in \mathscr {A}}\left\{ \sum _{\vec {k'}+\vec {k''}=\vec {k}}\sum _{\alpha ',\alpha '' \in \mathscr {A}}e^{i(\vec {k}\cdot \vec {x}-s(\sigma ^{(\alpha ')}_{(\vec {k'})}+\sigma ^{(\alpha '')}_{(\vec {k''})}-\sigma ^{(\alpha )}_{(\vec {k})}))} \times \right. \nonumber \\&\qquad \left. \times B^{(\alpha ',\alpha '',\alpha )}_{(\vec {k'},\vec {k''},\vec {k})} a^{(\alpha ')}_{(\vec {k'})}(t) a^{({\alpha ''})}_{(\vec {k''})}(t) \right\} \phi ^{(\alpha )}_{(\vec {k})}. \end{aligned}$$

(E.4)

Only three wave resonances can survive inside the fast averaging equation, and we define the set $\mathscr {S}_{\alpha ,\vec {k}}$ as survival index set:

$$\begin{aligned} \mathscr {S}_{\alpha ,\vec {k}} = \left\{ (\vec {k'}, \vec {k''}, \alpha ', \alpha '' ) \vert \vec {k'} +\vec {k''} = \vec {k}, \sigma ^{(\alpha ')}_{(\vec {k'})}+\sigma ^{(\alpha '')}_{(\vec {k''})}=\sigma ^{(\alpha )}_{(\vec {k})} \right\} . \end{aligned}$$

(E.5)

Directly projecting (3.6) onto the slow mode $\phi ^0$ will focus our attention on the analysis of slow component dynamics and its evolution equation. Verification on resonant triad interactions under the index set $\mathscr {S}_{0,\vec {k}}$ will be operated as follows (for both $V_r = 1$ and $V_r = \varepsilon ^{-1}$), which will illuminate the decoupling relationship between slow and fast components.

For $V_r = 1$ case, we turn to eigenvectors set (D.21), where $\phi ^{(0)}, \phi ^{(q)}$ are known as slow modes while $\phi ^{(+)}, \phi ^{(-)}$ are fast since previous two are associated with zero frequencies and later two own nonzero frequencies. When we confine that the resonant triad interactions involve at least one slow mode $\phi ^{(0)}$ (slow–$(*)$–$(*)$ impact), all possible resonant interactions coefficient B under the survival index set $\mathscr {S}_{0,\vec {k}}$ are

$$\begin{aligned} B^{(+,-,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(-,+,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(q,q,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(q,0,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(0,q,0)}_{(\vec {k'},\vec {k''},\vec {k})}=0 . \end{aligned}$$

(E.6)

Similar concrete form can be formulated for the linear operator $\mathscr {L}_0$ and simply yields

$$\begin{aligned} \lim _{\tau \rightarrow \infty }\dfrac{1}{\tau }\int ^{\tau }_0 e^{s\mathscr {L}_*}\mathscr {L}_0(e^{-s\mathscr {L}_*}\bar{v}(\vec {x},t))\mathrm{d}s=\sum _{\vec {k}\in \mathbb {Z}^3}\sum _{\sigma ^{(\alpha ')}_{(\vec {k})}=\sigma ^{(\alpha )}_{(\vec {k})}} L^{(\alpha ',\alpha )}_{(\vec {k})} a^{(\alpha ')}_{(\vec {k})}(t)e^{i\vec {k}\vec {x}}\phi ^{(\alpha )}_{(\vec {k})} , \end{aligned}$$

(E.7)

where $L^{(\alpha ',\alpha )}_{(\vec {k})}=\left\langle A_0(\vec {k})\phi ^{(\alpha ')}_{(\vec {k})},\phi ^{(\alpha )}_{(\vec {k})}\right\rangle$ is the coefficient for linear operator $\mathscr {L}_0$ and $A_0(\vec {k})$ is (D.13). Direct calculation gives following two inner product for $\alpha = 0$ (Note that we only need to check two cases $\alpha ' = q$ and $\alpha ' = 0$ when $\alpha = 0$ since only $\sigma ^{(0)}_{\vec {(k)}} - \sigma ^{(0)}_{\vec {(k)}} = 0$ and $\sigma ^{(0)}_{\vec {(k)}} - \sigma ^{(q)}_{\vec {(k)}} = 0$.)

$$\begin{aligned} \left\langle A_0(\vec {k})\phi ^{(0)}_{(\vec {k})}\right. , \left. \phi ^{(0)}_{(\vec {k})}\right\rangle =\left\langle A_0(\vec {k})\phi ^{(q)}_{(\vec {k})}\right. , \left. \phi ^{(0)}_{(\vec {k})}\right\rangle = 0 \Rightarrow L_{(\vec {k})}^{(q,0)} = L_{(\vec {k})}^{(0,0)} = 0. \end{aligned}$$

(E.8)

Finally for $\phi ^{0}$ mode, the explicit limiting dynamic evolution equation (derived from projecting (3.6) into $\phi ^{0}$ mode) expressed as an ODE of its amplitude $a^0_{\vec {k}}$ are given as follows (by setting $\alpha = 0$ in (E.4, E.7)),

$$\begin{aligned} \dfrac{\hbox {d}a^{(0)}_{(\vec {k})}}{\hbox {d}t}+\sum _{\begin{matrix} {k'}+\vec {k''}=\vec {k}\\ \sigma ^{(\alpha ')}_{(\vec {k'})}+\sigma ^{(\alpha '')}_{(\vec {k''})}=\sigma ^{(0)}_{(\vec {k})} \end{matrix}} B^{(\alpha ',\alpha '',0)}_{(\vec {k'},\vec {k''},\vec {k})} a^{(\alpha ')}_{(\vec {k'})} a^{(\alpha '')}_{(\vec {k''})}+\sum _{\sigma ^{(\alpha ')}_{(\vec {k})}=\sigma ^{(0)}_{(\vec {k})}} L^{(\alpha ',0)}_{(\vec {k})}a^{(0)}_{(\vec {k})}=0. \end{aligned}$$

(E.9)

We remind the reader that orthogonality is not guaranteed in previous eigenvectors (D.21), however, the reason one could still process the ODE analysis of $a^0_{\vec {k}}$ by successfully projecting (3.6) on $\phi ^0$ mode is because that $\phi ^{0}$ is perpendicular to the rest of three modes $\phi ^{q}, \phi ^{+}, \phi ^{-}$. Together with the resonant coefficient calculation showed above in (E.6) and linear term coefficient (E.8), one may observe that the slow mode ($\phi ^0$) is free of interactions with the fast modes. In other words, the amplitudes $a_{\vec {k}}^{0}$ is well determined only by itself in the limiting fast wave averaging equation (3.6):

$$\begin{aligned} \dfrac{\hbox {d}a^{(0)}_{(\vec {k})}}{\hbox {d}t}+\sum _{\begin{matrix} {k'}+\vec {k''}=\vec {k}\\ \sigma ^{(0)}_{(\vec {k'})}+\sigma ^{(0)}_{(\vec {k''})}=\sigma ^{(0)}_{(\vec {k})} \end{matrix}} B^{(0,0,0)}_{(\vec {k'},\vec {k''},\vec {k})} a^{(0)}_{(\vec {k'})} a^{(0)}_{(\vec {k''})}=0 . \end{aligned}$$

(E.10)

An inversion transformation of the Fourier-space equation for slow mode $\phi ^0$ leads to the conservation of equivalent potential voriticity. Technically speaking, the fast-wave-averaging equation for $PV_e$ in purely saturated region with $V_r = 1$ is given by

$$\begin{aligned} \dfrac{D}{D t}PV_e = \left( \dfrac{\partial }{\partial t} + \vec {u}_{(PV_e)} \cdot \nabla \right) PV_e = 0, \end{aligned}$$

(E.11)

implying that slow mode ($PV_e$ or $\phi ^0$) evolves independently from fast mode (waves or $\phi ^{\pm }$) under the presence of water and rainfall. The subscript $(PV_e)$ indicates that a variable has been computed by inverting from $(M,PV_e,W_1,W_2,u_m,v_m)$ to $(\vec {u},\theta _e,q_t)$ using $(PV_e)$ only. From the perspective of Fourier space, one may treat $\vec {u}_{(PV_e)}$ as the contribution only from the entries in slow mode $\phi ^0$.

For $V_r = \varepsilon ^{-1}$ case, eigenvectors set (D.26) will be used to process analysis. In contrast with $V_r = 1$ case, only one mode $\phi ^0$ with zero eigenvalue remains to be slow. Similar algebra states the following resonant interactions coefficient B under the survival index set $\mathscr {S}_{0,\vec {k}}$ and linear term coefficient L as follows

$$\begin{aligned}&B^{(+,-,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(-,+,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(+,q,0)}_{(\vec {k'},\vec {k''},\vec {k})}=B^{(q,+,0)}_{(\vec {k'},\vec {k''},\vec {k})}=0 , \end{aligned}$$

(E.12)

$$\begin{aligned}&\left\langle A_0(\vec {k})\phi ^{(0)}_{(\vec {k})}\right. , \left. \phi ^{(0)}_{(\vec {k})}\right\rangle = 0 \Rightarrow L_{(\vec {k})}^{(0,0)} = 0. \end{aligned}$$

(E.13)

Hence, in the remarkable resonant triad interactions only slow-slow-slow impact survives. The possibility of slow-fast-fast has been eliminated by (E.12), meanwhile, slow-fast-slow, slow-slow-fast aren’t counted since no resonant interaction is generated from them ($(\vec {k'},\vec {k''},slow,fast)\not \in \mathscr {S}_{0,\vec {k}}$). In conclusion, $V_r = \varepsilon ^{-1}$ gives the same result as (E.10) and (E.11).