Kerry Emanuel Lorenz Center, Department of Earth, Atmospheric, and Planetary
Sciences,
Massachusetts Institute of Technology, Cambridge, Massachusetts

# On the Effect of Surface Friction and Upward Radiation of Energy on
Equatorial Waves
Submitted for peer review.

###### Abstract

In theoretical models of tropical dynamics, the effects of both surface
friction and upward wave radiation through interaction with the stratosphere
are oft-ignored, as they greatly complicate mathematical analysis. In this
study, we relax the rigid-lid assumption and impose surface drag, which allows
the barotropic mode to be excited in equatorial waves. In particular, a
previously developed set of linear, strict quasi-equilibrium tropospheric
equations is coupled with a dry, passive stratosphere, and surface drag is
added to the troposphere momentum equations. Theoretical and numerical model
analysis is performed on the model in the limits of an inviscid surface
coupled to a stratosphere, as well as a frictional surface under a rigid-lid.
This study confirms previous research that shows the presence of a
stratosphere strongly shifts the growth rates of fast propagating equatorial
waves to larger scales, reddening the equatorial power spectrum. The growth
rates of modes that are slowly propagating and highly interactive with cloud-
radiation are shown to be negligibly affected by the presence of a
stratosphere. Surface friction in this model framework acts as purely a
damping mechanism and increases the poleward extent of the equatorial waves
through barotropic vorticity generation. Numerical solutions of the coupled
troposphere-stratosphere model with surface friction also show that the
barotropic mode can be tropospherically trapped when excited by surface
friction but in the presence of a highly stratified stratosphere. The
superposition of phase-shifted barotropic and first baroclinic modes is also
shown to lead to an eastward vertical tilt in the dynamical fields of Kelvin-
wave like modes.

##  1 Introduction

Reduced models of the tropical atmosphere have found much success in
replicating many characteristics of the tropical circulation. Of prominence
are theoretical models that reduce the linear, primitive equations to the
shallow water equations by use of only the first baroclinic mode (Matsuno,
1966; Gill, 1980). Indeed, remarkable evidence of the linear and neutral
equatorial waves that arise from the Matsuno-Gill model have been documented
in the tropics (Wheeler and Kiladis, 1999). The first baroclinic mode has also
been used extensively in simple models of the tropics, from studies of
intraseasonal oscillations (Emanuel, 1987; Sobel et al., 2001) to steady
circulations (Neelin and Held, 1987; Emanuel et al., 1994; Neelin and Zeng,
2000; Sobel and Bretherton, 2000), among many others.

While the first baroclinic mode is certainly a prominent feature of the
tropical atmosphere (Xu and Emanuel, 1989), observational data and analysis
have suggested the existence of another vertical mode, the second (stratiform)
baroclinic mode (Mapes and Houze, 1995; Straub and Kiladis, 2002). These
observations have lead to a plethora of theoretical studies that analyze how
the interaction between the first and second baroclinic modes can lead to
instability in equatorial waves (Mapes, 2000; Kuang, 2008b).

Perhaps curiously left behind is the barotropic mode, even though it does
indeed survive the strict quasi-equilibrium assumption; in fact strict quasi-
equilibrium eliminates all but the barotropic and first baroclinic modes
(Emanuel, 1987; Neelin and Zeng, 2000). The barotropic mode can be excited in
linear models of equatorial waves through coupling to the first baroclinic
mode via surface friction (Wang, 1988; Wang and Rui, 1990; Moskowitz and
Bretherton, 2000) or removal of the rigid-lid assumption (Yano and Emanuel,
1991). In this study, we build on the linear, strict quasi-equilibrium model
first formulated in Khairoutdinov and Emanuel (2018) (henceforth, KE18) and
further analyzed in Emanuel (2020) (henceforth, E20), and investigate how
excitation of the barotropic mode through both surface friction and coupling
to the stratosphere affects the growth, structure, and propagation of
equatorial waves.

Baked into the modal decomposition of many simple models of the tropics is the
rigid lid assumption, since it dramatically simplifies analytic solutions.
However, in reality the tropopause does not act as a rigid lid on the
troposphere. While it is true the stratosphere has a larger stratification
than the troposphere, the stratification in the stratosphere is not infinite.
A “leaky lid” analogy is more accurate, as wave energy can radiate to the
stratosphere.

A few studies have investigated the impact of a stratification jump at the
tropopause on the evolution of a wave in the atmosphere. Yano and Emanuel
(1991) extends the tropical intraseasonal model introduced in Emanuel (1987)
by imposing a dry, passive stratosphere above the troposphere, and found that
adding a stratosphere strongly damps the smallest scale $v=0$ waves, shifting
the growth rates to the larger scales. By imposing wave-radiation boundary
conditions, other studies have found similar effects: that the effect of the
stratosphere is a dampening one (Moskowitz and Bretherton, 2000; Kuang,
2008a). Chumakova et al. (2013) further investigates the leaky lid effect by
deriving a set of vertical, dissipative modes using the 2-D, linear,
Boussinesq equations, overlaying a stratosphere with buoyancy frequency
$N_{2}$ over a troposphere with buoyancy frequency $N_{1}$, where
$N_{2}>N_{1}$. In their model, a new barotropic-like mode appears with a fast
damping time scale. However, the vertical modes of Chumakova et al. (2013) are
not orthogonal; they have also been criticized as unphysical, since the energy
is unbounded with height, and they also do not admit steady state solutions to
steady state heating (Edman and Romps, 2017). Regardless of the exact
specifics in applying a leaky-lid above a troposphere, the inclusion of a
stratosphere tends to shift growth rates of unstable modes to larger scales
and allows the barotropic mode to be excited.

There have also been many studies that have investigated the role surface
friction plays in modifying equatorial waves, and more prominently, the
Madden-Julian Oscillation (MJO). These theories based on CISK (conditional
instability of the second kind), in which surface friction was postulated to
act, through its induced moisture convergence, as a destabilizing mechanism
for convectively coupled equatorial waves (Wang, 1988; Wang and Rui, 1990;
Moskowitz and Bretherton, 2000). Wang (1988) and Wang and Rui (1990)
formulated a 2.5 layer “frictional WAVE-CISK” model (2 tropospheric layers and
a thin frictional boundary layer) in which the barotropic mode can be excited
through surface friction. The surface friction acts to induce vertical motion
at the top of the boundary layer, which can amplify wave-disturbances if
correlated with temperature anomalies. CISK theories have received much
criticism [see Emanuel et al. (1994) and Neelin and Yu (1994)], since they
violate the convective statistical equilibrium hypothesis, where the rate of
production of convectively available potential energy (CAPE) by the large-
scale environment is very nearly balanced by its consumption via convection
(Arakawa and Schubert, 1974). In the context of surface friction, this means
that surface convergence is a by-product and not a driver of convection.
Indeed, observations in the tropics support the statistical equilibrium model
for convection in the tropics (Betts, 1982; Xu and Emanuel, 1989). Further,
numerical simulations of large-scale equatorial waves do not support the idea
that surface friction acts as a destabilizing mechanism for large-scale
equatorial waves (Chao and Chen, 2001).

Surface friction and boundary layer convergence have also been cited as one
mechanism for moistening of the lower troposphere by shallow upward motion
east of the MJO center (Wang, 1988; Hsu and Li, 2012; Adames and Wallace,
2014). In theoretical models, this mechanism has been shown to influence the
propagation speed and growth rates of the MJO through modulation of the gross
moist stability (Sobel and Maloney, 2013; Adames and Kim, 2016).

In this study, we will show how surface friction and interaction with the
stratosphere through upward wave energy radiation modifies the characteristics
of equatorial waves. The rigid-lid assumption is removed by explicitly
coupling a passive and dry stratosphere on top of a convecting troposphere,
and the energy density of solutions is enforced to decay with height. Surface
drag is imposed on a thin boundary layer at the surface. In particular, we
will focus on how the two aforementioned mechanisms can excite the barotropic
mode. More uniquely, the barotropic mode in this study does not separate
convective heating from large-scale thermodynamics, and thus it does not
violate the strict quasi-equilibrium hypothesis.

The paper is organized as follows. Section 2 describes the linear model.
Section 3 presents the solutions of the linear model. The paper concludes with
a discussion and summary in section 4.

##  2 Linear model

In this section, we describe and formulate the governing equations of our
linear model. The model formulation is separated into two parts, section 22.1
which describes the tropospheric model, and section 22.2 which describes the
stratospheric model.

###  2.1 Strict quasi-equilibrium troposphere

KE18 and E20 derive and analyze an equation set for a linear system that
describes the dynamics and thermodynamics of an atmosphere that maintains a
vertically constant saturation moist entropy $s^{*}$ in the free troposphere.
However, in both of those studies, the authors assume a rigid-lid and
frictionless surface. Hence, in their model, only the baroclinic mode can be
excited, and upward radiation into the stratosphere is absent.

Here we derive nearly equivalent dynamics, but remove the rigid-lid hypothesis
and include the barotropic mode. To begin, we first apply a Galerkin
decomposition of the vertical modes of the troposphere and truncate all modes
except the first two basis functions, $V_{0}$ and $V_{1}$, which are defined
as the barotropic and baroclinic modes, respectively (Neelin and Zeng, 2000).
Mathematically, they are:

$\displaystyle V_{0}$ $\displaystyle=1$ (1) $\displaystyle V_{1}$
$\displaystyle=\frac{\overline{T}(p)-[\overline{T}]}{T_{b}-[\overline{T}]}$
(2)

where $\overline{T}$ is the basic state temperature, $T_{b}$ is the boundary
layer temperature, $[\overline{T}]$ is the pressure-weighted vertical average
of temperature. The operator $[\>\bullet\>]=\frac{1}{\Delta
p}\int_{p_{t}}^{p_{s}}\bullet\>\>dp$ is the pressure-weighted vertical average
along a moist adiabat, where $p_{s}$ is the surface pressure, $p_{t}$ is the
tropopause pressure, and ${\Delta p=p_{s}-p_{t}}$. As is standard for vertical
modes, the basis functions are orthogonal, or
$\int_{p_{s}}^{p_{t}}V_{0}V_{1}dp=0$. Furthermore, note that $[V_{1}]=0$. From
this vertical mode decomposition, we assume separable dependencies between the
horizontal and vertical modes as follows:

$\displaystyle\phi(x,y,p,t)$
$\displaystyle=\phi_{0}(x,y,t)V_{0}+\phi_{1}(x,y,t)V_{1}(p)$

and likewise for the other prognostic variables.

For a strict quasi-equilibrium troposphere in which the saturation moist
entropy $s^{*}$ is constant with height, linearized geopotential perturbations
are directly connected to $s^{*}$ perturbations (Emanuel, 1987).

$\frac{\partial\phi^{\prime}}{\partial p}=-\bigg{(}\frac{\partial T}{\partial
p}\bigg{)}_{s^{*}}s^{*\prime}$ (3)

where prime superscripts indicate perturbation quantities. The above may be
directly integrated from the surface upwards to yield:

$\phi^{\prime}(p)=\phi^{\prime}_{b}+s^{*\prime}(\overline{T}_{b}-\overline{T}(p))$
(4)

where $\phi^{\prime}_{b}$ is the geopotential in the boundary layer. When non-
dimensionalized (see Appendix A for details), Equation (4) yields:

$\phi^{\prime}(p)=\phi_{b}^{\prime}+(1-V_{1})s^{*\prime}$ (5)

Note, the geopotential can be separated into its barotropic and baroclinic
components:

$\displaystyle\phi_{0}(x,y,t)$
$\displaystyle=(\phi_{b}^{\prime}+s^{*\prime})V_{0}$ (6)
$\displaystyle\phi_{1}(x,y,p,t)$ $\displaystyle=-s^{*\prime}\,V_{1}$ (7)

The pressure-weighted vertical average operator is also applied to Equation
(4) to give:

$\phi^{\prime}_{b}=[\phi]^{\prime}+s^{*\prime}([\overline{T}]-\overline{T}_{b})$
(8)

which in non-dimensional form is:

$\phi^{\prime}_{b}=\phi^{\prime}_{0}-s^{*\prime}$ (9)

Unlike the purely baroclinic motions described in Khairoutdinov and Emanuel
(2018), the geopotential now contains contributions from the barotropic mode.
Note that Equation (5) can be evaluated at the tropopause and combined with
Equation (9) to obtain the tropopause geopotential
$\phi^{\prime}_{\text{tp}}$, which will be required to couple the system to
the stratosphere:

$\phi^{\prime}_{\text{tp}}=\phi^{\prime}_{0}-V_{1}(p_{t})s^{*\prime}$ (10)

where $p_{t}$ is the non-dimensional tropopause pressure.

Next, we formulate the full equation of motion on an equatorial $\beta$-plane,
adding in surface friction, which is represented as applying drag on an
infinitesimally small boundary layer.

$\displaystyle\frac{D\bm{V}}{Dt}$
$\displaystyle=-\nabla\phi-\hat{k}\times\beta
y\bm{V}-\delta(p-p_{s})\frac{C_{d}}{h_{b}}|\bm{V}|\bm{V}$ (11)

where $\bm{V}$ is the vector wind, $\beta$ is the meridional gradient of the
Coriolis force, $\delta$ is the Dirac delta function, $p$ is pressure, $C_{d}$
is the drag coefficient, and $h_{b}$ is the boundary layer depth. Note the
surface stress is parameterized using the bulk aerodynamic drag formula.
Linearizing around surface easterlies, non-dimensionalizing according to
details in Appendix A, substituting in Equation (5), and dropping all primes
of perturbation quantities, we obtain:

$\displaystyle\frac{\partial u}{\partial t}$
$\displaystyle=-\frac{\partial\phi_{b}}{\partial x}+(1-V_{1})\frac{\partial
s^{*}}{\partial x}+yv-2Fu\delta(p-p_{s})$ (12)
$\displaystyle\frac{1}{\delta_{x}}\frac{\partial v}{\partial t}$
$\displaystyle=-\frac{\partial\phi_{b}}{\partial y}+(1-V_{1})\frac{\partial
s^{*}}{\partial y}-yu-\frac{F}{\delta_{x}}v\delta(p-p_{s})$ (13)

where $F$ is the non-dimensional surface friction coefficient and $\delta_{x}$
represents the magnitude of zonal geostrophy [corresponding to $\delta$ in
KE18]. Finally, we project the linearized horizontal momentum equations onto
the barotropic and baroclinic modes:

$\displaystyle\frac{\partial u_{0}}{\partial t}$
$\displaystyle=-\frac{\partial\phi_{0}}{\partial x}+yv_{0}-2F(u_{0}+u_{1})$
(14) $\displaystyle\frac{1}{\delta_{x}}\frac{\partial v_{0}}{\partial t}$
$\displaystyle=-\frac{\partial\phi_{0}}{\partial
y}-yu_{0}-\frac{F}{\delta_{x}}(v_{0}+v_{1})$ (15) $\displaystyle\frac{\partial
u_{1}}{\partial t}$ $\displaystyle=\frac{\partial s^{*}}{\partial
x}+yv_{1}-2F(u_{0}+u_{1})$ (16)
$\displaystyle\frac{1}{\delta_{x}}\frac{\partial v_{1}}{\partial t}$
$\displaystyle=\frac{\partial s^{*}}{\partial
y}-yu_{1}-\frac{F}{\delta_{x}}(v_{0}+v_{1})$ (17)

Next, we enforce mass continuity through the continuity equation in pressure
coordinates:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}+\frac{\partial\omega}{\partial p}=0$ (18)

where $\omega$ is the pressure vertical velocity. Integrating the continuity
equation from the surface to the tropopause and using the fact that
$[V_{1}]=0$:

$\int_{p_{s}}^{p_{t}}\nabla_{H}\cdot\bm{V}=\frac{\partial u_{0}}{\partial
x}+\frac{\partial v_{0}}{\partial y}=\omega(p_{t})$ (19)

where we have used zero vertical velocity condition at the lower boundary, and
$\omega(p_{t})$ is the tropopause vertical velocity. Equation (19) shows that
the tropopause vertical velocity is only a function of the divergence of the
barotropic mode, as the baroclinic mode is zero at the tropopause, by
definition.

Finally, the thermodynamic equations in the troposphere link the dynamics to
the thermodynamics, and are only slightly modified from KE18 and E20 in that
horizontal diffusion is removed and $\kappa=1$:111The second author discovered
that $\kappa$, an additional non-dimensional coefficient that scales the
cloud-radiative feedback and was defined in KE18 and E20, must be equal to 1
for consistency of the non-dimensional scaling.

$\displaystyle\frac{\partial s^{*}}{\partial t}$
$\displaystyle=(1+C)s_{m}-w-\alpha u_{b}-\chi s^{*}$ (20)
$\displaystyle\gamma\frac{\partial s_{m}}{\partial t}$
$\displaystyle=-Ds^{*}-\alpha u_{b}-Gw+Cs_{m}$ (21)

where $s_{m}$ is a characteristic moist entropy of the free troposphere,
$w=-\frac{\partial u_{b}}{\partial x}-\frac{\partial v_{b}}{\partial y}$ is a
proxy for the mid-level vertical velocity based on the boundary layer zonal
velocity, $u_{b}=u_{0}+u_{1}$, and boundary layer meridional velocity
$v_{b}=v_{0}+v_{1}$. The non-dimensional coefficients $C$, $\alpha$, $\chi$,
$D$, $G$, and $\gamma$ are described and formulated in detail in KE18.
Briefly, $C$ represents the strength of cloud radiative feedback, $\alpha$ is
the wind-induced surface heat exchange (WISHE) feedback parameter, $\chi$ is
boundary layer damping, $D$ is entropy damping, $G$ is the gross moist
stability, and $\gamma$ modifies the time scale of tropospheric entropy.

Equations (14) - (17), (19)- (21) formulate the tropospheric system, where
$u_{0}$, $v_{0}$, $u_{1}$, $v_{1}$, $\phi_{0}$, $s$, $s_{m}$ are the unknown
variables. Note that the linear system is not complete: additional
specification of the vertical velocity at the tropopause is required to
complete the system. As mentioned previously, studies that assume a rigid-lid
set the tropopause velocity to be zero. Other studies parameterize the
tropopause dynamics using a wave-radiation upper boundary condition (Moskowitz
and Bretherton, 2000; Kuang, 2008a). In this study, we take a different
approach and couple the vertical velocity to an explicit stratosphere model,
as will be derived in section 22.2.

###  2.2 Coupling to the stratosphere

In this section, we couple a dry, passive stratosphere to the strict quasi-
equilibrium troposphere described in section 22.1. We choose to represent a
dry and passive stratosphere using the linearized, inviscid primitive
equations in log-pressure coordinates and in hydrostatic balance [see Chapter
3 of Andrews et al. (1987)]:

$\displaystyle\frac{\partial u^{\prime}_{s}}{\partial
t}=-\frac{\partial\phi^{\prime}_{s}}{\partial x}+\beta yv^{\prime}_{s}$ (22)
$\displaystyle\frac{\partial v^{\prime}_{s}}{\partial
t}=-\frac{\partial\phi^{\prime}_{s}}{\partial y}-\beta yu^{\prime}_{s}$ (23)
$\displaystyle\frac{\partial u^{\prime}_{s}}{\partial x}+\frac{\partial
v^{\prime}_{s}}{\partial
y}+\frac{1}{\rho_{s}}\frac{\partial(\rho_{s}w^{*\prime}_{s})}{\partial
z^{*}}=0$ (24) $\displaystyle\frac{\partial}{\partial
t}\frac{\partial\phi^{\prime}_{s}}{\partial z^{*}}+w^{*\prime}_{s}N^{2}=0$
(25)

where subscripts of $s$ denote quantities in the stratosphere, $w^{*}_{s}$ is
the log-pressure vertical velocity, $N^{2}$ is the buoyancy frequency,
$\rho_{s}$ is the basic state density, and the log-pressure vertical
coordinate ${z^{*}\equiv-H\ln(p/p_{t})+1}$ is defined such that $z^{*}=1$ is
the bottom boundary, or the tropopause. Equations (22)-(25) are non-
dimensionalized according to notation shown in Appendix A, with the additional
specification that the non-dimensional density decays exponentially with a
scale height $H_{s}$. The resulting, non-dimensional equations are shown in
Equations (26) - (30), with primes removed from perturbation quantities.

$\displaystyle\frac{\partial u_{s}}{\partial
t}=-\frac{\partial\phi_{s}}{\partial x}+yv_{s}$ (26)
$\displaystyle\frac{1}{\delta_{x}}\frac{\partial v_{s}}{\partial
t}=-\frac{\partial\phi_{s}}{\partial y}-yu_{s}$ (27)
$\displaystyle\frac{\partial u_{s}}{\partial x}+\frac{\partial v_{s}}{\partial
y}+\frac{1}{\rho_{s}}\frac{\partial(\rho_{s}w^{*}_{s})}{\partial z^{*}}=0$
(28) $\displaystyle\frac{\partial}{\partial t}\frac{\partial\phi_{s}}{\partial
z^{*}}+w^{*}_{s}S=0$ (29)
$\displaystyle\rho_{s}=\exp\Big{(}\frac{H}{H_{s}}(1-z^{*})\Big{)}$ (30)

where $S$ is a non-dimensional stratospheric stratification.

It is important to note that these equations form a complete system by
themselves. The stratospheric linear system admits neutral equatorial wave
solutions under a rigid-lid upper boundary condition. However, under a upward
wave radiation boundary condition, all of the solutions decay exponentially in
time since there is no forcing in the stratosphere model and wave energy
escapes upwards. Growing solutions that satisfy the upward wave radiation
boundary condition in the stratosphere are possible, however, if there is
mechanical forcing from the troposphere via the tropopause. In order to
investigate these kinds of solutions, we must couple the troposphere system
with the stratosphere system in a consistent fashion. Classical coupling
conditions require continuity of normal stress across the interface, and
continuity of normal displacement to the fluid interface. Since the free-
troposphere is modeled as inviscid, the first condition simplifies to
continuity of pressure:

$\phi_{s}(x,y,z^{*}=1,t)=\phi(x,y,p=p_{t},t)$ (31)

Since there is no imposed shear across the tropopause, the second condition
implies continuity of vertical velocity:

$w^{*}_{s}(x,y,z^{*}=1,t)=-B\omega(x,y,p=p_{t},t)$ (32)

where $B=((p_{s}-p_{t})/p_{t})(H_{s}/H)>0$ is a conversion coefficient between
pressure coordinates and log-pressure coordinates.

###  2.3 Full, linear model

The troposphere system [Equations (14) - (17), (19) - (21)] is coupled to the
stratosphere system [Equations (26) - (30)] through the two matching
conditions [Equations (31) and (32)]. Altogether, these formulate a complete
linear system, in which growing solutions whose energy decays to zero as
$z^{*}\rightarrow\infty$ represent modes that grow in the troposphere and
propagate vertically into the stratosphere.

##  3 Solutions

The full linear model is a complex system that cannot be easily solved
theoretically. However, analyzing the model in the limits of (1) an inviscid
surface with coupling to the stratosphere, and (2) a frictional surface under
a rigid-lid, allows us to isolate the impacts of both mechanisms. Solutions of
the full model, with active surface friction and stratosphere coupling, are
then analyzed to illuminate their combined effects.

###  3.1 Leaky modes

In this section, we focus first on solutions of the purely leaky modes, with
no surface friction ($F=0$). The solutions are analyzed separately: modes
where $v=0$ and higher order meridional modes where there is non-zero
meridional velocity. In what follows, unless otherwise stated, the primes are
dropped from the linear perturbation variables.

####  3.1.1 $v=0$ modes

Although the full linear model is extremely complex, restricting the solutions
to $v=0$ allows for tractable analytical insight. In the troposphere, we
assume solutions of the form:

$u_{0}=\hat{U}_{0}(y)\exp(ikx+\sigma t)$ (33)

where $k$ is the zonal wave number, $\sigma$ is the complex growth rate, and
capitalized variables with hat notations are the meridional structure
functions. Equivalent forms are assumed for $\phi_{0}$, $u_{1}$, $s$, $s_{m}$.
Solutions in the stratosphere are assumed of the form:

$u_{s}=\frac{\hat{U_{s}}(y)}{\sqrt{\rho_{0}}}\exp(ikx+\sigma t+im(z^{*}-1))$
(34)

where $m$ is the complex vertical wavenumber. Equivalent forms are also
assumed for $\phi_{s}$, and $w^{*}_{s}$. As in KE18, the meridional structure
$Y$ of the $v=0$ modes in the troposphere can be derived by combining
Equations (14) - (17):

$\frac{\partial}{\partial
y}\Big{(}\hat{S}-\hat{\Phi}_{0}\Big{)}=\frac{iky}{\sigma}\Big{(}\hat{S}-\hat{\Phi}_{0}\Big{)}$
(35)

Next, we combine Equations (14), (16), (20), and (21) to eliminate $u_{0}$,
$u_{1}$, and $s_{m}$ and obtain the relationship between the meridional
function of the saturation moist entropy and the barotropic geopotential:

$\displaystyle\hat{S}=\lambda\hat{\Phi}_{0}$ (36)

where

$\displaystyle\lambda$ $\displaystyle=\frac{ika_{2}+k^{2}a_{3}}{\sigma
a_{1}+ika_{2}+a_{3}k^{2}}$ (37) $\displaystyle a_{1}$ $\displaystyle\equiv
D(1+C)+(\chi+\sigma)(\gamma\sigma-C)$ (38) $\displaystyle a_{2}$
$\displaystyle\equiv\alpha(\gamma\sigma+1)$ (39) $\displaystyle a_{3}$
$\displaystyle\equiv\gamma\sigma+(G-1)C+G$ (40)

Combining Equations (35) and (36) gives us the meridional structure $Y$ of the
troposphere portion of the $v=0$ modes:

$Y=\exp\Big{(}\frac{ik}{2\sigma}y^{2}\Big{)}$ (41)

which is equivalent in form to the meridional structure of the modes in the
rigid-lid case. Furthermore, only solutions with an eastward phase speed
satisfy the boundary conditions in $y$.

Next, we move on to solving the portion of the mode that exists in the
stratosphere. With the solution form shown in Equation (34), Equations (26) -
(29) reduce to:

$\displaystyle\sigma\hat{U}_{s}+ik\hat{\Phi}_{s}=0$ (42) $\displaystyle
y\hat{U}_{s}+\frac{\partial\hat{\Phi}_{s}}{\partial y}=0$ (43) $\displaystyle
ik\hat{U}_{s}+\frac{m^{2}\sigma}{S}\hat{\Phi}_{s}=0$ (44)

where $\hat{W}_{s}(y)=-\frac{im\sigma}{S}\hat{\Phi}_{s}(y)$, assuming
$m^{2}>>(H/2H_{s})^{2}$, which is a vertical short-wave approximation.
Equations (42) and (44) combine into the well-known dispersion relation for
the Kelvin-wave:

$\sigma=\pm i\frac{\sqrt{S}k}{m}$ (45)

As shown in Yano and Emanuel (1991), the associated group velocity is:

$c_{g,z}=\frac{k\sqrt{S}}{|m|^{4}}\Big{(}\text{Real}(m)^{2}-\text{Imag}(m)^{2}\Big{)}$
(46)

which indicates that the vertical group velocity increases with zonal
wavenumber.

Next, Equations (42) and (43) combine to give a meridional structure $Y$ that
is equivalent to that of the troposphere shown in Equation (41). For solutions
that obey the meridional boundary conditions, or that the mode amplitudes go
to zero as $y\rightarrow\pm\infty$, we must have that $\sigma_{i}<0$
(equivalent to an eastward phase speed), so we choose the positive root.

We now apply the matching conditions to derive the dispersion relation. The
continuity of pressure condition demands that:

$\phi_{s}(x,y,z^{*}=1,t)=\phi_{0}-V_{1}(p_{t})s^{*}$ (47)

Equations (14) and (19) combine to give the vertical velocity at the
tropopause in the troposphere, which must be equal to the vertical velocity at
the tropopause in the stratosphere:

$\omega(p_{t})=\frac{k^{2}}{\sigma}\phi_{0}=-\frac{1}{B}w^{*}_{s}=\frac{1}{B}\frac{im\sigma}{S}\hat{\phi}_{s}$
(48)

where we note the equivalence of a rigid-lid ($\omega(p_{t})=0$) to the
absence of the barotropic mode. Combining with Equations (10), (36), (45), and
(47), Equation (48) reduces to the dispersion relation:

$\underbrace{\sigma a_{1}+a_{4}}_{\text{rigid
lid}}+\underbrace{\sigma(kB\sqrt{S})^{-1}(\sigma a_{1}+\nu
a_{4})}_{\text{correction}}=0$ (49)

where $\nu=1-V_{1}(p_{t})$ and $a_{4}=ika_{2}+a_{3}k^{2}$. As shown in the
underbraces, Equation (49) is written in the form of the rigid-lid dispersion
relation, plus a quartic-order correction term whose magnitude is inversely
propotional to the square root of the stratosphere stratification, $S$. It is
clear that in the limit of $S\rightarrow\infty$, the dispersion relation
reduces to that of the rigid-lid case. The quartic dispersion relation is
solved numerically and the solutions are checked rigorously to satisfy the
governing equations, boundary conditions, and matching conditions.

Before examining the solutions of the leaky $v=0$ modes, we first briefly
derive the functional form of $\omega$ in the troposphere, in the case where
$v=0$. We start by taking the time derivative of the non-dimensional
continuity equation and substitute the zonal velocity using Equations (14) and
(16):

$\sigma\frac{\partial\omega}{\partial p}=k^{2}\phi$ (50)

Taking a derivative in pressure allows us to substitute in Equation (3) and
yields:

$\sigma\frac{\partial^{2}\omega}{\partial p^{2}}=-k^{2}\bigg{(}\frac{\partial
T}{\partial p}\bigg{)}_{s^{*}}s^{*}$ (51)

Integrating in pressure once returns:

$\sigma\frac{\partial\omega}{\partial
p}=-k^{2}s^{*}\overline{T}(p)\Big{|}_{s^{*}}+C$ (52)

where $C$ is an integration constant that must be determined through the
boundary conditions. We integrate from the surface to the tropopause:

$\sigma(\omega(p_{t})-\omega(p_{s}))=(p_{s}-p_{t})(-k^{2}s^{*}[\overline{T}]+C)$
(53)

which allows us to determine $C$:

$C=\sigma\frac{\omega(p_{t})-\omega(p_{s})}{p_{s}-p_{t}}+k^{2}s^{*}[\overline{T}]$
(54)

Substituting for $C$ into Equation (52), setting $\omega(p_{s})=0$ as earlier,
and integrating once in pressure yields the non-dimensional vertical structure
of $\omega$:

$\omega(p)=\underbrace{\frac{p_{s}-p}{p_{s}-p_{t}}\omega(p_{t})}_{\text{barotropic}}-\underbrace{\frac{k^{2}}{\sigma}s^{*}\int_{p_{s}}^{p_{t}}V_{1}\>dp}_{\text{baroclinic}}$
(55)

As shown in the underbraces, the vertical structure of the pressure vertical
velocity is a sum of the barotropic and baroclinic modes. The barotropic
component changes linearly with pressure, and the baroclinic mode is zero at
the surface/tropopause and peaks in the mid-troposphere. The superposition of
the barotropic and baroclinic mode can lead to a vertical tilt in the vertical
velocity profile that depends on the phase lag between the two modes. It is
worth noting that by definition, the baroclinic mode cannot interact with the
stratosphere; in this model, it is only through the excitement of the
barotropic mode that waves can radiate energy into the stratosphere.

Figure 1: Non-dimensional (a) frequency, (b) growth rate, (c) phase speed and
for the fastest growing $v=0$ mode. Vertical group velocity of the leaky mode
is shown in black in the bottom panel, while rigid-lid solutions are shown in
blue and leaky-lid solutions shown in red. Non-dimensional parameters are
$\alpha=1.5$, $\chi=0.5$, $C=0$, $\gamma=2$, $D=0.5$, $G=0.1$, $S=100$.

Figure 1 compares the non-dimensional frequency, growth rates, phase speeds,
and vertical group velocity of the leaky $v=0$ modes to those of the rigid-lid
modes, with non-dimensional coefficients $\alpha=1.5$, $\chi=0.5$, $C=0$,
$\gamma=2$, $D=0.5$, $G=0.1$, $S=100$. These parameters somewhat reflect
Earth-like conditions in the tropics, with the exception of no cloud radiation
interaction, and were specifically chosen to examine the branch of solutions
that closely resemble the classical Kelvin-wave solutions but are instead
amplified through the WISHE feedback. Figure 1a shows that the frequencies of
the leaky waves are only slightly larger than their rigid-lid counterparts,
and as such the phase speeds are slightly faster (Figure 1c). The modification
of the growth rates tell a much different story: Figure 1b shows strong
damping of the growth rate of smaller-scale waves, shifting the power spectrum
towards larger scale waves, which is consistent with the results of Yano and
Emanuel (1991). This reddening of the power spectrum can be physically
explained through the large vertical group velocities of the smaller scale
waves, as shown in Figure 1c and indicated in Equation (46). The smaller scale
waves propagate their energy very quickly into the stratosphere, dampening
their growth rate. For instance, the $k=10$ wave has a vertical group velocity
that nearly exceeds its zonal phase speed.

Figure 2: Eigenfunction of the $v=0$, $k=1$ mode corresponding to the
parameters detailed in Figure 1, on a (left) horizontal cross-section at the
surface, and (right) vertical cross-section centered on the equator. Contours
indicate the saturation entropy or geopotential, where solid (dashed) lines
indicate positive (negative) perturbations. Arrows indicate zonal velocity
perturbations, and color shading indicates vertical velocity or $\omega$
perturbations, where indicated. Magenta line outlines the tropopause.

Figure 2 shows a summary of the $k=1$ eigenfunction using identical non-
dimensional parameters as detailed in Figure 1. The horizontal cross-section
at the surface shows a Kelvin-wave pattern with surface easterlies maximizing
east of the maximum temperature anomaly. Vertical velocity maximizes west of
the peak temperature anomalies, and the WISHE feedback is responsible for wave
growth and enhanced eastward propagation. The horizontal structure is not
significantly modified from that of the rigid-lid. The vertical structure of
the mode exhibits a bit more complexity. Unlike the vertical structure of the
rigid-lid model, which is purely baroclinic, the superposition of the
barotropic and baroclinic modes leads to an eastward vertical tilt, as shown
in Figure 2 and Equation (55). The mode has a downward phase propagating,
eastward tilted component in the stratosphere, which is consistent with a
Kelvin-wave that has upward vertical energy propagation. As indicated by the
dispersion relation, the vertical wavelength in the stratosphere is controlled
by the horizontal wavenumber and stratification; a stronger stratification or
shorter horizontal wavelength decreases the vertical wavelength.

The tropospheric system also allows for solutions of significantly slower
propagating modes when cloud-radiation interaction is turned on, as shown in
KE18 and E20. We select the non-dimensional parameters $\alpha=1$, $\chi=1$,
$C=2.5$, $\gamma=2$, $D=1$, $G=0.02$, $S=100$, in order to obtain the slow
modes as the fastest growing solutions. Note that since we eliminated
horizontal diffusion of MSE as compared to the model described in KE18, the
growth rates no longer peak at low wavenumbers. Figure 3 shows the frequency,
growth rate, phase speed, and vertical group velocities for the $v=0$ slow
modes that can interact with the stratosphere. Aside from a small modification
of the wave properties at the highest wave numbers, these slow modes are not
significantly affected by the presence of a stratosphere. Their vertical group
velocities are almost negligibly small, and the corresponding vertical
wavelengths in the stratosphere are extremely short. These slow modes are
trapped in the troposphere and do not leak much energy into the stratosphere.
Thus, these results suggest that for slow propagating modes, such as the MJO,
a rigid-lid assumption is an accurate approximation with regards to the growth
rates.

t

Figure 3: Analogous to Figure 1, but for the fastest-growing mode with non-
dimensional parameters $\alpha=1$, $\chi=1$, $C=2.5$, $\gamma=2$, $D=1$,
$G=0.02$, $=100$.

Next, we investigate the barotropic mode magnitude and phase tilt across a
range of non-dimensional parameters. The phase lag and amplitude relationship
between the baroclinic and barotropic zonal velocities is:

$u_{0}=\Big{(}\frac{kB\sqrt{S}}{\sigma\lambda}-\nu\Big{)}u_{1}$ (56)

Equation (56) shows that the amplitude of the barotropic mode decreases with
increasing zonal wavenumber and stratospheric stratification. For Earth-like
parameters, the amplitude of the barotropic mode for the $k=1$ solution is
around an order of magnitude smaller than the amplitude of the baroclinic
mode. The phase lag is determined by the complex coefficient in parentheses in
Equation (56).

Figure 4: (a, c) Barotropic zonal velocity to total zonal velocity ratio
[$u_{0}$ / ($u_{0}+u_{1}$)], and (b, d) phase lead of the barotropic mode with
respect to the baroclinic mode, in radians, as functions of stratosphere
stratification, $S$, and horizontal wavenumber. (a, b) Top row generated using
the non-dimensional parameters described in Figure 1 (WISHE-modified Kelvin
waves), and (c, d) bottom row generated using the non-dimensional parameters
described in Figure 3 (slow $v=0$ modes). White area indicates non-growing
modes.

Figure 4a and Figure 4c show the barotropic zonal velocity to total zonal
velocity ratio, which we define as $u_{0}/(u_{0}+u_{1})$, as well as the phase
lead of the barotropic mode to the baroclinic mode, for a wide range of
stratosphere stratification $S$ (an Earth-like range of stratosphere
stratification is $S\approx 20-200$). In general, we observe decreasing
amplitude of the barotropic mode with increasing stratosphere stratification
and zonal wavenumber, as suggested by Equation (56). For both the classical
WISHE-driven Kelvin waves (Figure 4, top row) and radiatively destabilized
slow modes (Figure 4, bottom row), the barotropic mode is around an order of
magnitude smaller than the baroclinic mode for Earth-like parameters.
Furthermore, the barotropic mode leads the baroclinic mode by approximately
0.1-0.2 radians. The phase lead of the barotropic mode is always positive,
indicating that the eastward phase tilt shown in Figure 2 is robust across a
wide-range of scales and $S$. In addition, an exploration of the non-
dimensional parameter space did not find any solutions where the barotropic
mode lagged the baroclinic mode (not shown). This result is consistent with
physical intuition. For the WISHE-driven $v=0$ modes, the baroclinic
easterlies peak east of the maximum temperature perturbation. If the
barotropic zonal velocities are instead westward of the baroclinic zonal
velocities, then the total surface easterlies could peak west of the maximum
temperature anomalies, which is inconsistent with an WISHE-driven eastward
propagating mode. Note, this is inconsistent with observational analyses of
convectively coupled Kelvin waves, in which a westward tilt with height is
observed (Straub and Kiladis, 2002). This may indicate the importance of the
second baroclinic mode for Kelvin waves (Mapes, 2000; Kuang, 2008a).

Finally, it is worth noting that the $u_{0}$ magnitude is larger than the
$u_{1}$ magnitude for extremely leaky stratospheres (Figure 4); though this
may not be relevant to Earth’s atmosphere, a barotropic mode that exceeds its
baroclinic counterpart leads to top-heavy vertical velocity profiles and fast
propagation into the stratosphere.

####  3.1.2 Higher order modes

Without the $v=0$ approximation, the solution set is extremely complex. In the
rigid-lid case, the meridional structure of the eigenmodes of the full
equations are parabolic cylinder functions of degree $n$, where $n$ is the
order of the Hermite polynomial. A brief mathematical analysis of the
(intractable) leaky-lid solutions suggests that the meridional structure
functions of the linear solutions are sums of parabolic cylinder functions.
Note that by convention, we denote $n=-1$ solutions as those where $v=0$.

For these higher order modes, we instead solve the linear problem using
numerical code, which is only appropriate for finding growing modes. The
troposphere system [Equations (14) - (17), (20) - (21)] is discretized in $y$,
while the stratosphere system [Equations (26) - (29)] are discretized in $y$
and $z$. Solutions are assumed to have a zonal structure of the form
$\exp(ikx)$. In order to integrate the linear system in time, the stratosphere
equations need to be transformed to a set of linear prognostic equations. The
mass continuity equation [Equation (28)] is first integrated from the lower
boundary in $z^{*}$, to obtain:

$\rho_{s}w^{*}_{s}(y,z^{*})=w^{*}_{s}(y,z^{*}=1)-\int_{z^{*}=1}^{z}\Big{[}\rho_{s}\Big{(}iku_{s}(y,z^{*})+\\\
\frac{\partial}{\partial y}v_{s}(y,z^{*})\Big{)}\Big{]}dz^{*}$ (57)

where the tropopause velocity $w^{*}(z^{*}=1)$ is equal to the vertical
velocity at the tropopause in the troposphere equations, as required from the
vertical velocity matching condition. Given the vertical velocity, Equation
(29) can be used to calculate the geopotential:

$\frac{\partial}{\partial t}\frac{\partial}{\partial
z}\phi(y,z^{*})=-w^{*}(y,z)S$ (58)

Integrating from the upper boundary downwards gives the prognostic equation
for $\phi_{s}$:

$\frac{\partial}{\partial
t}\phi_{s}(y,z^{*})=-\int^{z}_{\infty}w^{*}_{s}(y,z^{*})S\>dz^{*}+C_{0}$ (59)

where $C_{0}$ must be determined by the upper boundary condition. An upwards
wave-radiation boundary condition could be applied to determine $C_{0}$, as in
Moskowitz and Bretherton (2000), but such a condition does not exist for the
slower modes that have interactive cloud radiation and no clear analog in the
stratosphere. While it is not necessary that $\phi_{s}=0$ as
$z\rightarrow\infty$ (as long as the energy density goes to zero), we have
included a Newtonian damping in a sizeable layer extending from the upper
boundary, and as such, we set $C_{0}=0$. Another advantage of the Newtonian
damping in the upper part of the domain is to eliminate any nonphysical
downward propagating modes. For realistic values of $S$, however, the
amplitude of the modes typically decay to zero at the top of the numerical
domain. Finally, $\phi_{s}(y,z^{*}=1)$ is connected to the troposphere through
the continuity of pressure matching condition [Equation (47)].

Spatial derivatives are approximated using fourth-order central differences,
and the system is forward time-stepped using fourth order Runge-Kutta. After
specifying non-dimensional parameters, the corresponding rigid-lid solution is
used to initialize the troposphere domain, while the stratosphere is
initialized at rest. Since the unbalanced wave must undergo rapid gravity-wave
adjustment, several dampening mechanisms are used to eliminate spurious
gravity-wave energy. A spectral filter is applied at each time step to
eliminate small-scale noise. A strong sponge-layer is included along the edges
and top of the domain to eliminate reflection of gravity waves, downward
propagating waves, and spurious noise. For details on the mathematical form of
the full numerical system, see Appendix B. The system is integrated for a long
period of time, during which the domain is periodically rescaled by a constant
to prevent numerical overflow. We then isolate the growing mode of interest
and infer both the complex growth rate and meridional/vertical structures.
Although we initially assume that the damping at the upper boundary does not
significantly affect the stratospheric solution, the corresponding growth rate
and the meridional structures of all prognostic variables are rigorously
checked to satisfy the governing equations, boundary conditions, and matching
conditions. As a partial test of correctness of the numerical code, the
numerical solutions for the $v=0$ ($n=-1$) modes were cross-referenced with
the $v=0$ analytic solutions.

The numerical solution is robust for linear solutions that are the fastest
growing equatorially symmetric or asymmetric mode. However, for solutions that
are not the fastest growing mode, initial unbalanced energy and/or numerical
error that projects onto the fastest growing mode symmetric/asymmetric mode
will cause the slower growing mode to be overtaken by the fastest growing mode
before the domain energy is concentrated solely in the mode of interest. For
instance, the second fastest growing equatorially symmetric mode can be
overtaken by the fastest growing equatorially symmetric mode. In these cases,
the meridional/vertical structure of the fastest growing mode is isolated,
projected onto the domain output, and removed from the fields. If the second
fastest growing mode does not have a weak growth rate, the remaining, filtered
fields will contain the mode of interest for a long enough period of time to
infer the complex growth rate and structure.

##### WISHE classical modes

We first compare the leaky wave solutions to the rigid-lid solutions for the
WISHE-destabilized classical modes of the equatorial waveguide. This is done
by choosing $\alpha=3.5$, $\chi=0.5$, $C=0$, $\gamma=1$, $D=2.5$, $G=0.25$,
$\delta_{x}=15$, $S=75$; note that $C=0$ eliminates the slower propagating
modes observed in KE18 and E20. Figure 5 shows the non-dimensional growth
rates and phase speeds for select meridional orders $n=-1,0,1,2$ leaky
equatorial waves. The growth rates of the Kelvin, westward and eastward mixed
Rossby-gravity, and inertia-gravity waves are clearly dampened in relation to
their rigid-lid counterparts. The strength of dampening, in a percent relative
sense, grows stronger with increasing wavenumber for the eastward modes. This
is in contrast to the almost negligible dampening of the westward propagating
modes as $|k|$ increases. The growth rates of the westward propagating $n=0$
and $n=1$ waves are not strongly affected by the stratosphere. The phase
speeds of the leaky-waves are almost negligibly faster than their rigid-lid
counterparts, with the largest, though still slight, modifications observed
for the $n=-1$ solutions. The growth rates for smaller scale eastward
propagating waves are further dampened with reduced stratification in the
stratosphere (compare $S=75$ to $S=25$). The effect is quite pronounced; in
the case shown in Figure 5, the smallest scale rigid-lid modes possess the
largest growth rates, but imposing a leaky stratosphere greatly shifts the
growth rates towards larger scales. For the $S=25$ case, the growth rate peaks
at $k=3$ for the $n=-1,0$ modes, and $k\approx 5$ for the $n=1,2$ modes. These
results are consistent with the qualitative behavior shown in the mathematical
analysis of the leaky $v=0$ modes. As such, the solutions suggest that the
stratosphere acts as a reddener of the equatorial power spectrum, especially
for eastward propagating waves.

Figure 5: (Left) Non-dimensional growth rate for the $n=-1,0,1,2$ modes of the
(solid) rigid-lid system, the (dashed) leaky-lid system with $S=75$, and the
(dot-dashed) leaky-lid system with $S=25$. (Right) Non-dimensional phase speed
for the same modes, in the (solid) rigid-lid and (black-dashed) leaky-lid
system with $S=75$. $k<=3$ for $n=2$, and $k<=0$ for $n=1$ for $S=25$
solutions are omitted since they do not grow rapidly enough to infer the
complex growth rate. Non-dimensional parameters selected for these modes are
$\alpha=3.5$, $\chi=0.5$, $C=0$, $\gamma=1$, $D=2.5$, $G=0.25$,
$\delta_{x}=15$.

##### Slow modes

We next switch to a non-dimensional parameter set in which cloud radiative
feedbacks are turned on, such that the slow modes analyzed in E20 are the
fastest growing modes. In particular, we choose $\alpha=1$, $\chi=1$, $C=2.5$,
$\gamma=2$, $D=1$, $G=0.02$, $\delta_{x}=30$, $S=100$. Although the westward
propagating slow modes do grow in time under the rigid-lid model, they take
too long to reach a steady state and are quickly overtaken by the eastward
propagating modes. We thus focus on the eastward propagating slow modes, which
have a much clearer analog in the real atmosphere (i.e. the MJO).

Figure 6: (Left) Non-dimensional growth rate for the $n=-1,1,2$, eastward
propagating slow modes of the (solid) rigid-lid system, and the (dashed)
leaky-lid system with $S=100$ inferred from the numerical solution. (Right)
Non-dimensional frequency for the same modes, in the (solid) rigid-lid and
(dashed) leaky-lid system with $S=100$. Black dots indicate the analytic
solutions of the growth rate and frequencies of the $n=-1$ slow modes. All
$n=0$ solutions are not growing. Non-dimensional parameters selected for these
modes are $\alpha=1$, $\chi=1$, $C=2.5$, $\gamma=2$, $D=1$, $G=0.02$,
$\delta_{x}=30$.

Figure 6 shows the non-dimensional growth rates and frequencies for the
eastward propagating $n=-1,1,2$ slow modes. We first compare the numerical
$n=-1$ solutions to the analytic solutions (Figure 6, black dots). We observe
that the numerical model slightly dampens the growth rate, and only negligibly
modifies the frequencies. This slight damping could be attributed to the
spectral filter and/or numerical error. Assuming that the numerical error
affects the wave characteristics of the $n=1,2$ slow modes the same way, we
can infer that like the growth rates and frequencies of the leaky $n=-1$ slow
modes, the growth rates and phase speeds of the higher order slow modes are
not affected by the presence of the stratosphere.

Perhaps of greater interest is how the slow modes interact with the
stratosphere through the barotropic mode. We focus on the eastward propagating
$n=1$, $k=2$ mode which has a horizontal structure that closely resembles that
of the observed MJO (Emanuel, 2020). Figure 7 shows the horizontal and
vertical structure of the $n=1$, $k=2$ mode, where we observe surface cyclonic
gyres westward and poleward of the area of maximum ascent. Strong westerlies
are also observed westward of the maximum ascent region. The boundary layer
horizontal structure very closely resembles that of the rigid-lid solution
(not shown). At the tropopause (Figure 7, top right), the horizontal structure
is almost nearly the opposite of the horizontal structure at the surface,
indicating the prominence of the first baroclinic mode. However, the vertical
velocity at the tropopause is non-zero, and the pattern extends poleward into
the cyclonic/anti-cyclonic gyres, an indication of the presence of the
barotropic mode. If we move further up into the stratosphere, at around 20 km
(Figure 7, bottom left), the signature of the equatorial portion of the
tropospheric wave disappears, and the poleward cyclonic/anti-cyclonic gyres
become the prominent pattern. These gyres have a westward tilt with height, as
indicated in the vertical cross section at $y=2.5$ (Figure 7, bottom right),
though their amplitudes decay exponentially with height and are much smaller
with respect to the amplitudes of the corresponding gyres in the troposphere.
At the equator, the $n=1$, $k=2$ slow mode has a very weak eastward tilt with
height (not shown), since the barotropic mode leads the baroclinic mode. Like
the $v=0$ mode shown in Figure 2, in the stratosphere, the slow mode also
tilts eastward with height at the equator (not shown). This is consistent with
the stratospheric eastward tilt with height of the zonal wind on the equator
of observational MJO composites [see fig. 3 of Kiladis et al. (2005)].

Figure 7: Eigenfunction of the $n=1$, $k=2$ mode corresponding to the
parameters detailed in Figure 6, on a horizontal cross-section at the (top
left) the surface, (top right) tropopause [15 km], (bottom left) 18 km.
Vertical cross section eigenfuction at $y=2.5$ is shown on the bottom right.
Contours indicate the geopotential, where solid (dashed) lines indicate
positive (negative) perturbations. Arrows indicate wind perturbations, and
color shading indicate vertical velocity perturbations at the level indicated,
except for the boundary layer cross-section, where color shading indicates
mid-level vertical velocity.

###  3.2 Frictionally modified modes

In this section, we consider how surface friction can also excite the
barotropic mode. First, we simplify the linear system in the rigid-lid limit,
but now with non-zero surface friction. From Equation (19), if
$\omega(p_{t})=0$, it follows that the barotropic mode is non-divergent.
However, the solution cannot be solved in terms of a non-divergent stream
function $\psi_{0}$, since $\delta_{x}$ scales the meridional momentum
equation but not the mass continuity equation. Instead, we take the divergence
of Equations (14) and (15), yielding:

$\Big{(}\delta_{x}\frac{\partial^{2}}{\partial
y^{2}}-k^{2}\Big{)}\phi_{0}=y\Big{(}ikv_{0}-\delta_{x}\frac{\partial
u_{0}}{\partial
y}\Big{)}-\delta_{x}u_{0}-2Fik(u_{0}+u_{1})-F\frac{\partial}{\partial
y}\Big{(}v_{0}+u_{1}\Big{)}$ (60)

From Equation (60), we see that a barotropic mode at rest cannot be excited if
there is no surface friction ($F=0$). Thus, surface friction acts to couple
the barotropic mode with the baroclinic mode; since the forcing is primarily
in the first baroclinic mode; it is the first baroclinic mode that excites the
barotropic mode.

The barotropic geopotential $\phi_{0}$ is solved by inverting Equation (60).
Equations (14) - (17), (20), (21) form the complete linear system. The system
is solved numerically by discretizing in $y$. Spatial derivatives are again
approximated using fourth-order central differences, and the system is forward
time-stepped using fourth order Runge-Kutta. As before, the frictionless,
rigid-lid solution is used to initialize the domain, after which we integrate
the system for a long period of time and then isolate the growing mode of
interest.

Figure 8: (Left) Non-dimensional growth rate for the $n=-1,0,1$ modes of the
(solid) rigid-lid system, the (dashed) surface friction + rigid-lid system
with $F=0.1$, and the (dot-dashed) surface friction + rigid-lid system
$F=0.25$. $n=-1$ modes differ from the traditional definition of $v=0$, in
that the meridional wind is allowed to be non-zero but small. (Right) Non-
dimensional phase speed for the same modes, in the (solid) rigid-lid and
(black-dashed) surface friction + rigid-lid system with $F=0.1$. All $n=2$,
and $n=1$ for $k<=0$ solutions are omitted since they do not grow rapidly
enough to infer the complex growth rate. Non-dimensional parameters selected
for these modes are $\alpha=3.5$, $\chi=0.5$, $C=0$, $\gamma=1$, $D=2.5$,
$G=0.25$, $\delta_{x}=15$.

First, we investigate the complex growth rates as a function of the non-
dimensional surface friction, $F$, for an Earth-like range is $F\approx
0.1-0.4$. Figure 8 shows the growth rate of the $n=-1$, $n=0$, and $n=1$ modes
for the same parameter choices as in Figure 5, which selects for WISHE-
modified Matsuno modes. It is important to note that the $n=-1$ modes in this
section differ from the traditional definition of $v=0$, in that the
meridional wind is allowed to be non-zero but small, since $v_{0}=v_{1}=0$
solutions do not exist under a fixed surface and rigid-lid assumption. Figure
8 indicates that the surface friction acts as a damping effect on all
wavelengths, and the strength of the damping is nearly constant across all
wavelength. This is expected as the aerodynamic drag law acts on a fixed
damping time scale. The phase speeds are not significantly changed from the
phase speeds of the rigid-lid modes.

Figure 9: (a) Meridional structure of 15 times the barotropic zonal wind and
the baroclinic zonal wind for the $n=-1$, $k=1$ mode with $F=0.25$ and the
parameter set chosen in Figure 8. (b) Same as (a) but for the barotropic and
baroclinic meridional wind. (c) Horizontal structure of the boundary layer
(shading) barotropic vorticity and (contour) baroclinic vorticity. (d)
Equatorial vertical cross sections of the (shading) pressure vertical velocity
and (contours) zonal wind perturbations for the same mode. Solid (dashed)
contours are positive (negative) zonal wind anomalies.

As illustrated in Equation (60), surface friction can lead to an excitation of
the barotropic mode. Figure 9a and Figure 9b show the meridional structure of
the $n=-1$, $k=1$ barotropic zonal and meridional velocities as compared to
their baroclinic counterparts. Like the leaky-lid case, the barotropic zonal
wind magnitude is around an order of magnitude smaller than the baroclinic
mode. However, while the baroclinic mode magnitudes are primarily confined
close to the equator (around $y=-5$ to $y=5$), the barotropic mode wind
velocities have long tails that extend far away from the equator. These long
tails can be understood from the barotropic mode vorticity equation, which can
be solved for by taking the curl of Equations (14) and (15):

$\sigma\zeta_{0}+v_{0}=-F\bigg{[}\frac{1}{\delta_{x}}\frac{\partial
v_{0}}{\partial x}-2\frac{\partial u_{0}}{\partial
y}+\frac{1}{\delta_{x}}\frac{\partial v_{1}}{\partial x}-2\frac{\partial
u_{1}}{\partial y}\bigg{]}$ (61)

where the barotropic mode vorticity
$\hat{\zeta}_{0}=\hat{\nabla}_{H}\times\vec{v}_{0}$, and
$\hat{\nabla}_{H}=[\frac{1}{\delta_{x}}\frac{\partial}{\partial
x},\frac{\partial}{\partial y}]$. From Equation (61), it is evident that the
barotropic mode vorticity can be excited by not only the baroclinic mode
vorticity, but also by itself. The generation of barotropic vorticity far from
the equator [see Figure 9c] leads to non-negligible barotropic winds far from
the equator. While further analysis of this effect is out of the scope of this
study, the barotropic mode might be an important teleconnection mechanism
between the tropics and the extratropics (Horel and Wallace, 1981). An
important caveat to note, however, is that a rigid-lid barotropic mode is
completely trapped in the troposphere; the addition of a leaky stratosphere
may limit the poleward extent of the barotropic mode, as will be examined in
the next section.

Under a rigid lid the non-divergent barotropic mode cannot be associated with
vertical velocity perturbations, though a vertical tilt can still exist in the
horizontal wind fields from the superposition of the barotropic and baroclinic
horizontal winds. Figure 9d shows the vertical cross section on the equator of
the $n=-1$, $k=1$ mode. The vertical structure of the vertical velocity is
only first baroclinic, despite the zonal wind exhibiting a slight eastward
vertical tilt. This eastward vertical tilt in the zonal wind field is
consistent with the eastward tilt observed in the purely leaky $n=-1$ modes.
This is because the barotropic mode leads the baroclinic mode, as can also be
seen in the horizontal structure of the boundary layer vorticity decomposition
in Figure 9c.

Unfortunately the numerical model is quite unstable for the slow modes in the
rigid-lid, surface friction limit. Despite strong sponge layers, instabilities
unrelated to the slow modes (spurious noise, gravity waves) develop quickly in
the numerical domain, precluding inferral of the complex growth rate and
horizontal structure. The behavior of the barotropic mode in the slow modes
will be discussed in the following section using the fully coupled, surface
friction model, albeit with a large stratosphere stratification, instead of a
rigid-lid.

These results show that surface friction act strictly as a damping mechanism
in our framework. This is because the heating associated with frictional
convergence is not divorced from the large-scale thermodynamics. Surface
friction does, however, act to modify both the horizontal and vertical
structure of the equatorial waves, as shown through the long poleward tails of
the barotropic mode and vertical tilt in the zonal wind fields.

###  3.3 Leaky, frictionally modified modes

Figure 10: (a, c) Non-dimensional growth rate and (b, d) non-dimensional phase
speed for the $n=-1,0,1$ modes of the (solid) rigid-lid system, the (dashed)
surface friction + rigid-lid system with $S=150$ and $F=0.1$, and the (dot-
dashed) surface friction + rigid-lid system with $S=75$ and $F=0.1$. Non-
dimensional parameters selected for (a, b) are $\alpha=3.5$, $\chi=0.5$,
$C=0$, $\gamma=1$, $D=2.5$, $G=0.25$, $\delta_{x}=15$, and for (c, d) are
$\alpha=1$, $\chi=1$, $C=2.5$, $\gamma=2$, $D=1$, $G=0.02$, $\delta_{x}=30$.
Only $S=75$ case is shown for (c, d).

While examining the excitation of the barotropic mode in the limits of surface
friction under a rigid lid and a frictionless surface under a leaky lid were
both useful exercises to isolate their respective effects on tropical waves,
in the real world, both surface friction and leakage of energy to the
stratosphere can act simultaneously to modify equatorial waves through
excitation of the barotropic mode. In particular, the long-tails of the
tropospherically trapped barotropic mode observed in the rigid-lid, surface
friction model may behave differently if the barotropic mode can leak energy
into the stratosphere.

To understand the extent to which both of these mechanisms can interact, we
run the full troposphere-stratosphere numerical model with non-zero surface
friction, using the same methods to infer the complex growth rates and
eigenmodes. We first examine their combined effect on the WISHE-driven
classical Matsuno modes through the same parameter set as used in Figures 5
and 8. In addition, we focus on the eastward propagating modes since the
westward propagating modes were not strongly affected by the presence of the
stratosphere.

Figure 10a and Figure 10b show the non-dimensional growth rate and phase
speeds for the WISHE destabilized Matsuno modes, respectively. Since there is
both energy leakage into the stratosphere and surface friction, the growth
rates are strongly dampened from the equivalent modes in the rigid-lid,
inviscid limit. The frequencies of the WISHE Matsuno modes are not greatly
modified, though all of the $n=0$ modes and smaller scale $n=-1$ modes ($k\geq
5$) have slower phase speeds than their rigid-lid counterparts. Comparing with
the growth rates in the stratosphere-only case shown in Figure 5, we can see
that the damping effects of the stratosphere and surface friction are
approximately additive.

The more interesting question, perhaps, is if the barotropic mode behaves
differently when it can be excited by surface friction but also interact with
the stratosphere. To examine this, we decompose the horizontal divergence of
the barotropic mode into $\frac{\partial u}{\partial x}$, $\frac{\partial
v}{\partial y}$, and the sum of both, which is equivalent to the pressure
vertical velocity at the tropopause. Figure 11 shows the barotropic mode
velocities and the barotropic horizontal divergence decomposition for the
$n=0$, $k=1$ WISHE Matsuno mode, for a realistic stratosphere stratification
($S=75$), as well as a highly stratified one ($S=500$). In the case with a
highly stratified stratosphere, we see the long-tail feature of the barotropic
mode velocities, as also seen in the surface friction under a rigid-lid model.
However, a key differing feature is that the barotropic mode velocities are no
longer completely non-divergent, as evidenced in Figure 11b. Further, the
horizontal divergence of the barotropic mode (or, equivalently, the tropopause
vertical velocity) decays to zero very quickly polewards, around $|y|=3$. This
is because near the equator, $\frac{\partial u_{0}}{\partial x}$ and
$\frac{\partial v_{0}}{\partial y}$ have the same sign, evidence that the
barotropic mode is exciting the stratosphere near the equator. In contrast,
despite the barotropic velocities being small (though non-zero) polewards of
the equator, the horizontal divergence is zero because $\frac{\partial
u_{0}}{\partial x}$ and $\frac{\partial v_{0}}{\partial y}$ have opposite
signs and almost exactly cancel. This is evidence of a tropospherically
trapped barotropic mode. When we reduce the stratosphere stratification to
$S=75$, as shown in Figure 11c and Figure 11d, in effect allowing the
tropopause to be more leaky, we see that the polewards extent of the
barotropic mode velocities is greatly reduced, and hence the poleward extent
of $\frac{\partial u_{0}}{\partial x}$ and $\frac{\partial v_{0}}{\partial y}$
is also reduced. The horizontal divergence still decays to zero around
$|y|=3$, as in the case where $S=500$. Thus, when the troposphere is allowed
to leak more energy into the stratosphere, the barotropic mode is no longer
trapped in the troposphere, and the long-tail of the barotropic mode
velocities is reduced. Finally, while we have taken a large stratification
($S=500$) to emphasize the characteristics of the trapped barotropic mode, the
extent of the trapping of the barotropic mode increases with $S$. Thus, in
general, the polewards extent of the barotropic mode increases when the
stratosphere stratification increases.

Figure 11: (a, c) Barotropic mode zonal and meridional velocities for the
WISHE modified Matsuno $n=0$, $k=1$ mode using the same parameter set as in
Figure 8, $\alpha=3.5$, $\chi=0.5$, $C=0$, $\gamma=1$, $D=2.5$, $G=0.25$,
$\delta_{x}=15$, $F=0.1$, but with (a) $S=500$ and (c) $S=75$. (b, d)
$2\frac{\partial u_{0}}{\partial x}$, $2\frac{\partial v_{0}}{\partial y}$
(scaled by 2) and $\nabla\cdot\vec{v_{0}}$ for the same mode but with (b)
$S=500$ and (d) $S=75$.

Figure 10c and Figure 10d show the non-dimensional growth rate and phase
speeds for the slow modes, for which we obtain by choosing the same parameter
set as used in Figure 3 and 6. Note, the $n=0$ waves do not grow in time. The
growth rates are not modified greatly from the rigid lid solutions; surface
friction seems to have a weaker damping effect on the slower propagating
modes. The damping effect of the stratosphere on the growth rates is small but
the greatest at the largest scales, and almost negligible for the smaller
scale waves, as before. However, the frequencies/phase speeds for the slow
modes are reduced, by as much as 30% for the $k=1$ wave. The damping of the
frequency is much greater for the $n=1$ slow propagating modes than the $n=-1$
modes, and much stronger for larger scale $n=1$ waves. It is also worth noting
that the $S=150$ case is not shown for the slow modes since the lines are
nearly indistinguishable from the $S=75$ case.

The behavior of the poleward extent of the barotropic mode for the slow modes
is similar to that for the WISHE-modified Matsuno modes. Figure 12 shows that
the barotropic mode associated with the $n=1$, $k=2$ slow mode becomes non-
divergent away from the equator and trapped in the troposphere. However, one
key difference is that the slow modes do not leak much energy into the
stratosphere to begin with, and thus decreasing the stratosphere
stratification does not completely eliminate the long-tailed behavior of the
barotropic velocities. This is evidenced by the small magnitude of the
horizontal divergence of the barotropic mode (Figure 12b, Figure 12d),
regardless of the stratosphere stratification, as $iku_{0}$ and
$\frac{\partial v_{0}}{\partial y}$ nearly cancel close to the equator
($|y|<3$), and exactly cancel for $|y|>3$.

Figure 12: Analogous to Figure 11 but for the slow propagating $n=1$, $k=2$
mode using the parameter set $\alpha=1$, $\chi=1$, $C=2.5$, $\gamma=2$, $D=1$,
$G=0.02$, $\delta_{x}=30$, $F=0.1$.

What are the physical parameters that control the non-dimensional, $S$, and
how can they vary across different equatorial waves? The non-
dimensionalization of $S$ [see Appendix A] suggests that the stratosphere
buoyancy frequency, $N^{2}$ and the meridional length scale, $L_{y}$, are
quantities that could potentially lead to large variations in $S$. $L_{y}$
depends on a multitude of factors that can vary greatly in the tropics, such
as the troposphere dry stratification, moist adiabatic lapse rate, and the
precipitation efficiency. The largest influencing factor, however, is the
precipitation efficiency, $\epsilon_{p}$: waves with larger $\epsilon_{p}$
experience a much greater stratosphere stratification. For instance, for
Earth-like parameters and a $\epsilon_{p}=0.5$, $S\approx 20$, while for
$\epsilon_{p}=0.95$, $S\approx 200$, an order of magnitude difference. While
our definition of $\epsilon_{p}$ is an egregiously simple parameterization of
cloud microphysics, the impact of $\epsilon_{p}$ on $S$ links cloud
microphysical properties to the behavior of the barotropic mode.

##  4 Discussion and summary

In this study, we extended a previously developed linear framework to include
two mechanisms that can excite the barotropic mode in equatorial waves,
surface friction, and coupling to the stratosphere. We first analyzed the
modification of equatorial waves in the two separate limits of (1) coupling to
the stratosphere with an inviscid surface, and (2) surface friction under a
rigid-lid. Using a combination of theoretical solutions for the $v=0$ mode and
numerical solutions for higher order meridional modes, we found that the
presence of a stratosphere leads to upward wave energy propagation that
strongly dampens the growth rate of smaller-scale waves. This effect is
consistent across the WISHE-modified Matsuno modes explored in this study,
though it does not appreciably affect the slow propagating modes that are
destabilized from cloud-radiative feedbacks. The barotropic mode is found to
lead the baroclinic mode for the $n=-1$ eastward propagating modes, leading to
a slight eastward tilt in the vertical that arises from the superposition of
the two modes. The eastward tilt was found to be robust across horizontal
scales and stratosphere stratification. In the limit of surface friction under
a rigid-lid, we found that surface friction acts only to dampen growth rates
by nearly a constant across all waves and zonal wavenumbers, which is
reasonable given that surface friction acts on a constant time scale of
damping in our framework. We also found that frictional excitation of the
barotropic mode leads to long tails of non-divergent barotropic zonal and
meridional velocities away from the equator. The barotropic mode is also found
to lead the baroclinic mode for the $n=-1$ eastward propagating modes, leading
to a vertical tilt in the zonal wind field, but not the vertical velocity
field since the barotropic mode is non-divergent under a rigid lid.

The combined effects of surface friction and coupling to the stratosphere were
analyzed using numerical solutions of the full linear model. We found that for
the growth rates of the waves, the damping effects of surface friction and
coupling to the stratosphere were approximately additive. Under a large non-
dimensional stratosphere stratification, both the WISHE-modified Matsuno modes
and slow modes exhibit tropospherically trapped barotropic modes away from the
equator. Thus, the non-divergent barotropic velocities extend far away from
the equator and are small but non-zero. When the non-dimensional stratosphere
stratification is reduced, the poleward extent of the barotropic mode is
greatly reduced. These results indicate that the dimensional variables that
influence the non-dimensional stratosphere stratification, such as the
buoyancy frequency in the stratosphere and precipitation efficiency, play key
roles in modulating the poleward extent of the barotropic mode. While this
study was restricted to theoretical analysis of the barotropic mode, future
work will attempt to find evidence of the barotropic mode in both
observational and numerical data.

This work models the interaction of equatorial waves with a zero mean flow
stratosphere, and is the first basic step to illuminate how a dynamically dry
and passive fluid influences the dynamics of a moist and convecting fluid
underneath. Extension of the framework developed in this study to a non-zero
mean flow in the stratosphere will be the subject of future work, which would
allow for an investigation into the extent to which upwards wave radiation can
explain the observed relationship between the MJO and Quasi-Biennial
Oscillation (QBO) (Yoo and Son, 2016).

Finally, it is worth discussing some of the short-comings of the modeling
framework. Our application of surface friction in an infinitesimally small
boundary layer with zero vertical velocity is exceptionally crude compared to
real-world frictional boundary layers. In addition, the strict-quasi
equilibrium approximation may not be as accurate for extremely short or high
frequency waves (Ahmed et al., 2021), as evidenced by the presence of what
looks like a second baroclinic mode in observations of convectively coupled
Kelvin waves (Straub and Kiladis, 2002). Regardless, both surface friction and
upward radiation of wave energy are shown to important mechanisms that
influence the horizontal structure, vertical structure, and growth rates of
equatorial waves, through excitation of the barotropic mode.

###### Acknowledgements.

The authors gratefully acknowledge the support of the National Science
Foundation through grant NSF ICER-1854929. Model source code, instructions,
and code to generate the figures from outputs of the numerical model are
available at github.com/linjonathan. All code and figures were generated using
Matlab. [A] Non-dimensionalization Here, we define the non-dimensional
scalings for the variables that appear in the full linear model. The scalings
for the tropospheric quantities $s^{*\prime}$, $s_{m}^{\prime}$, $\chi$,
$\alpha$, $\gamma$, $D$ and $G$ are identical to those described in the
appendix of KE18.

$\displaystyle x$ $\displaystyle\rightarrow ax$ (62) $\displaystyle y$
$\displaystyle\rightarrow L_{y}y$ (63) $\displaystyle p$
$\displaystyle\rightarrow(p_{s}-p_{t})p$ (64) $\displaystyle z^{*}$
$\displaystyle\rightarrow Hz^{*}$ (65) $\displaystyle t$
$\displaystyle\rightarrow\frac{a}{\beta L_{y}^{2}}t$ (66) $\displaystyle
L_{y}^{4}$
$\displaystyle\rightarrow\frac{\Gamma_{d}}{\Gamma_{m}}\frac{d\overline{s}_{d}}{dz}(T_{b}-[\overline{T}])H\frac{1-\epsilon_{p}}{\beta^{2}}$
(67) $\displaystyle u^{\prime}$
$\displaystyle\rightarrow\frac{aC_{k}|\overline{\textbf{V}}|}{H}u^{\prime}$
(68) $\displaystyle v^{\prime}$
$\displaystyle\rightarrow\frac{L_{y}C_{k}|\overline{\textbf{V}}|}{H}v^{\prime}$
(69) $\displaystyle w^{*\prime}$ $\displaystyle\rightarrow
C_{k}|\overline{\textbf{V}}|w^{*\prime}$ (70) $\displaystyle\phi^{\prime}$
$\displaystyle\rightarrow\frac{a\beta
L_{y}^{2}C_{k}|\overline{\textbf{V}}|}{H}\phi^{\prime}$ (71)
$\displaystyle\omega$
$\displaystyle\rightarrow\frac{C_{k}\overline{\textbf{V}}(p_{s}-p_{t})}{H}\omega$
(72) $\displaystyle N^{2}$
$\displaystyle\rightarrow\frac{\beta^{2}L_{y}^{4}}{H^{2}}S$ (73)
$\displaystyle F$
$\displaystyle\rightarrow\frac{aC_{d}|\overline{\textbf{V}}|}{\beta
L_{y}^{2}h_{b}}$ (74)

where most dimensional parameters are described in the main text and in KE18.
Parameters not defined in this text are the mean radius of the Earth $a$, the
dry adiabatic lapse rate $\gamma_{d}$, the moist adiabatic lapse rate
$\gamma_{m}$, the dry entropy stratification $\frac{d\overline{s}_{d}}{dz}$,
the precipitation efficiency $\epsilon_{p}$, and the enthalpy exchange
coefficient $C_{k}$. The terms on the left of the arrow are the dimensional
quantities, and those on the right are the non-dimensional quantities.

[B]

Numerical Model The full mathematical description of the numerical system
(including damping terms) is below:

$\displaystyle\frac{\partial u_{0}}{\partial t}$
$\displaystyle=-ik\big{[}\phi_{s}+V_{1}(p_{t})s^{*}\big{]}+yv_{0}-2F(u_{0}+u_{1})-ru_{0}$
(75) $\displaystyle\frac{\partial v_{0}}{\partial t}$
$\displaystyle=\delta_{x}\Big{[}-\frac{\partial}{\partial
y}\big{[}\phi_{s}+V_{1}(p_{t})s^{*}\big{]}-yu_{0}\Big{]}-F(v_{0}+v_{1})-rv_{0}$
(76) $\displaystyle\frac{\partial u_{1}}{\partial t}$
$\displaystyle=iks^{*}+yv_{1}-2F(u_{0}+u_{1})-ru_{1}$ (77)
$\displaystyle\frac{\partial v_{1}}{\partial t}$
$\displaystyle=\delta_{x}\Big{[}\frac{\partial s^{*}}{\partial
y}-yu_{1}\Big{]}-F(v_{0}+v_{1})-rv_{1}$ (78) $\displaystyle\frac{\partial
s^{*}}{\partial t}$ $\displaystyle=(1+C)s_{m}-w-\alpha(u_{0}+u_{1})-\chi
s^{*}-rs^{*}$ (79) $\displaystyle\gamma\frac{\partial s_{m}}{\partial t}$
$\displaystyle=-Ds^{*}-\alpha(u_{0}+u_{1})-Gw+Cs_{m}-rs_{m}$ (80)
$\displaystyle w$ $\displaystyle=-ik(u_{0}+u_{1})-\frac{\partial}{\partial
y}(v_{0}+v_{1})$ (81) $\displaystyle\frac{\partial u_{s}}{\partial t}$
$\displaystyle=-ik\phi_{s}+yv_{s}-ru_{s}$ (82) $\displaystyle\frac{\partial
v_{s}}{\partial t}$
$\displaystyle=\delta_{x}\Big{[}-\frac{\partial\phi_{s}}{\partial
y}-yu_{s}\Big{]}-rv_{s}$ (83) $\displaystyle\frac{\partial\phi_{s}}{\partial
t}$ $\displaystyle=-\int^{z}_{\infty}w^{*}_{s}S\>dz^{*}-r\phi_{s}$ (84)
$\displaystyle\rho_{s}w^{*}_{s}$
$\displaystyle=-B\Big{(}iku_{0}+\frac{\partial v_{0}}{\partial
y}\Big{)}-\int_{z^{*}=1}^{z}\Big{[}\rho_{s}\Big{(}iku_{s}(y,z^{*})+\frac{\partial}{\partial
y}v_{s}(y,z^{*})\Big{)}\Big{]}dz^{*}$ (85)

where all variables are defined in the main text with the exception of $r$,
which is the sponge coefficient for the sponge layer that is applied at the
boundaries of the domain.

## References

  * Adames and Kim (2016) Adames, Á. F., and D. Kim, 2016: The MJO as a dispersive, convectively coupled moisture wave: Theory and observations. J. Atmos. Sci., 73 (3), 913–941.
  * Adames and Wallace (2014) Adames, Á. F., and J. M. Wallace, 2014: Three-dimensional structure and evolution of the vertical velocity and divergence fields in the MJO. J. Atmos. Sci., 71 (12), 4661–4681.
  * Ahmed et al. (2021) Ahmed, F., J. D. Neelin, and Á. F. Adames, 2021: Quasi-Equilibrium and Weak Temperature Gradient Balances in an Equatorial Beta-plane Model. J. Atmos. Sci., 78 (1), 209–227.
  * Andrews et al. (1987) Andrews, D. G., J. R. Holton, and C. B. Leovy, 1987: Middle atmosphere dynamics. 40, Academic press.
  * Arakawa and Schubert (1974) Arakawa, A., and W. H. Schubert, 1974: Interaction of a cumulus cloud ensemble with the large-scale environment, Part I. J. Atmos. Sci., 31 (3), 674–701.
  * Betts (1982) Betts, A. K., 1982: Saturation point analysis of moist convective overturning. J. Atmos. Sci., 39 (7), 1484–1505.
  * Chao and Chen (2001) Chao, W. C., and B. Chen, 2001: The role of surface friction in tropical intraseasonal oscillation. Mon. Wea. Rev., 129 (4), 896–904.
  * Chumakova et al. (2013) Chumakova, L. G., R. R. Rosales, and E. G. Tabak, 2013: Leaky rigid lid: New dissipative modes in the troposphere. J. Atmos. Sci., 70 (10), 3119–3127.
  * Edman and Romps (2017) Edman, J. P., and D. M. Romps, 2017: Beyond the rigid lid: Baroclinic modes in a structured atmosphere. J. Atmos. Sci., 74 (11), 3551–3566.
  * Emanuel (2020) Emanuel, K., 2020: Slow modes of the equatorial waveguide. J. Atmos. Sci., 77 (5), 1575–1582.
  * Emanuel (1987) Emanuel, K. A., 1987: An air-sea interaction model of intraseasonal oscillations in the tropics. J. Atmos. Sci., 44 (16), 2324–2340.
  * Emanuel et al. (1994) Emanuel, K. A., J. David Neelin, and C. S. Bretherton, 1994: On large-scale circulations in convecting atmospheres. Quart. J. Roy. Meteor. Soc., 120 (519), 1111–1143.
  * Gill (1980) Gill, A. E., 1980: Some simple solutions for heat-induced tropical circulation. Quart. J. Roy. Meteor. Soc., 106 (449), 447–462.
  * Horel and Wallace (1981) Horel, J. D., and J. M. Wallace, 1981: Planetary-scale atmospheric phenomena associated with the Southern Oscillation. Mon. Wea. Rev., 109 (4), 813–829.
  * Hsu and Li (2012) Hsu, P.-c., and T. Li, 2012: Role of the boundary layer moisture asymmetry in causing the eastward propagation of the Madden–Julian oscillation. J. Climate, 25 (14), 4914–4931.
  * Khairoutdinov and Emanuel (2018) Khairoutdinov, M. F., and K. Emanuel, 2018: Intraseasonal variability in a cloud-permitting near-global equatorial aquaplanet model. J. Atmos. Sci., 75 (12), 4337–4355.
  * Kiladis et al. (2005) Kiladis, G. N., K. H. Straub, and P. T. Haertel, 2005: Zonal and vertical structure of the Madden–Julian oscillation. J. Atmos. Sci., 62 (8), 2790–2809.
  * Kuang (2008a) Kuang, Z., 2008a: Modeling the interaction between cumulus convection and linear gravity waves using a limited-domain cloud system–resolving model. J. Atmos. Sci., 65 (2), 576–591.
  * Kuang (2008b) Kuang, Z., 2008b: A moisture-stratiform instability for convectively coupled waves. Journal of the Atmospheric Sciences, 65 (3), 834–854.
  * Mapes (2000) Mapes, B. E., 2000: Convective inhibition, subgrid-scale triggering energy, and stratiform instability in a toy tropical wave model. J. Atmos. Sci., 57 (10), 1515–1535.
  * Mapes and Houze (1995) Mapes, B. E., and R. A. Houze, 1995: Diabatic divergence profiles in western Pacific mesoscale convective systems. J. Atmos. Sci., 52 (10), 1807–1828.
  * Matsuno (1966) Matsuno, T., 1966: Quasi-geostrophic motions in the equatorial area. Journal of the Meteorological Society of Japan. Ser. II, 44 (1), 25–43.
  * Moskowitz and Bretherton (2000) Moskowitz, B. M., and C. S. Bretherton, 2000: An analysis of frictional feedback on a moist equatorial Kelvin mode. Journal of the atmospheric sciences, 57 (13), 2188–2206.
  * Neelin and Held (1987) Neelin, J. D., and I. M. Held, 1987: Modeling tropical convergence based on the moist static energy budget. Mon. Wea. Rev., 115 (1), 3–12.
  * Neelin and Yu (1994) Neelin, J. D., and J.-Y. Yu, 1994: Modes of tropical variability under convective adjustment and the Madden–Julian oscillation. Part I: Analytical theory. J. Atmos. Sci., 51 (13), 1876–1894.
  * Neelin and Zeng (2000) Neelin, J. D., and N. Zeng, 2000: A quasi-equilibrium tropical circulation model—Formulation. J. Atmos. Sci., 57 (11), 1741–1766.
  * Sobel and Maloney (2013) Sobel, A., and E. Maloney, 2013: Moisture modes and the eastward propagation of the MJO. J. Atmos. Sci., 70 (1), 187–192.
  * Sobel and Bretherton (2000) Sobel, A. H., and C. S. Bretherton, 2000: Modeling tropical precipitation in a single column. J. Climate, 13 (24), 4378–4392.
  * Sobel et al. (2001) Sobel, A. H., J. Nilsson, and L. M. Polvani, 2001: The weak temperature gradient approximation and balanced tropical moisture waves. J. Atmos. Sci., 58 (23), 3650–3665.
  * Straub and Kiladis (2002) Straub, K. H., and G. N. Kiladis, 2002: Observations of a convectively coupled Kelvin wave in the eastern Pacific ITCZ. J. Atmos. Sci., 59 (1), 30–53.
  * Wang (1988) Wang, B., 1988: Dynamics of tropical low-frequency waves: An analysis of the moist Kelvin wave. J. Atmos. Sci., 45 (14), 2051–2065.
  * Wang and Rui (1990) Wang, B., and H. Rui, 1990: Dynamics of the coupled moist Kelvin–Rossby wave on an equatorial $\beta$-plane. J. Atmos. Sci., 47 (4), 397–413.
  * Wheeler and Kiladis (1999) Wheeler, M., and G. N. Kiladis, 1999: Convectively coupled equatorial waves: Analysis of clouds and temperature in the wavenumber–frequency domain. J. Atmos. Sci., 56 (3), 374–399.
  * Xu and Emanuel (1989) Xu, K.-m., and K. A. Emanuel, 1989: Is the tropical atmosphere conditionally unstable? Mon. Wea. Rev., 117 (7), 1471–1479.
  * Yano and Emanuel (1991) Yano, J.-I., and K. Emanuel, 1991: An improved model of the equatorial troposphere and its coupling with the stratosphere. J. Atmos. Sci., 48 (3), 377–389.
  * Yoo and Son (2016) Yoo, C., and S.-W. Son, 2016: Modulation of the boreal wintertime Madden-Julian oscillation by the stratospheric quasi-biennial oscillation. Geophys. Res. Lett., 43 (3), 1392–1398.