HTML

Similar to CS, the atmospheric model is a twolayer quasigeostrophic flow confined to a periodic
$\beta $ plane channel with zonal walls at$y = 0$ and$\pi L$ . The equations in pressure coordinates are:$$\frac{\partial }{{\partial t}}({\nabla ^2}{\psi ^1}) + J({\psi ^1}, {\nabla ^2}{\psi ^1}) + \beta \frac{{\partial {\psi ^1}}}{{\partial x}} =  {k'_d}{\nabla ^2}({\psi ^1}  {\psi ^3}) + \frac{{{f_0}}}{{\Delta p}}\omega, $$ (1) $$\begin{split} & \frac{\partial }{{\partial t}}({\nabla ^2}{\psi ^3}) + J\left({\psi ^3}, {\nabla ^2}{\psi ^3} + {f_0}\frac{h}{H}\right) + \beta \frac{{\partial {\psi ^3}}}{{\partial x}} \\ &\quad \quad = {k'_d}{\nabla ^2}({\psi ^1}  {\psi ^3})  \frac{{{f_0}}}{{\Delta p}}\omega  {k_d}{\nabla ^2}{\psi ^3}, \quad\quad\quad\quad\quad\quad \end{split}$$ (2) where
$x$ and$y$ are eastward and northward coordinates, respectively;$t$ is time,${\nabla ^2}$ is the horizontal Laplace operator,$J$ is the Jacobian operator;${\psi ^1}$ and${\psi ^3}$ are the geostrophic streamfunction fields at${p_1} = 250$ and${p_3} = 750$ hPa, respectively;$\omega = {{{\rm d}p}/{{\rm d}t}}$ is the vertical velocity;${f_0}$ is the Coriolis parameter at a central latitude${\phi _0} = 45^\circ N$ , with$\beta = {{{\rm d}f}/{{\rm d}y}}$ as its meridional gradient;$\Delta p = 500$ hPa is the pressure difference between the two layers;$H$ is mean depth of each layer;$h(x, y)$ is the lower boundary topographic height, and we assume that$h \ll H$ . The constants${k_d}$ and${k'_d}$ multiply the surface friction term and the internal friction between layers, respectively.We define
$$\psi = ({\psi ^1} + {\psi ^3})/2,\quad \theta = ({\psi ^1}  {\psi ^3})/2, $$ (3) then the atmospheric motion equations become the following:
$$ \begin{split} \!\!\! \frac{\partial }{{\partial t}}({\nabla ^2}\psi) & \!+\! J(\psi, {\nabla ^2}\psi) \!+\! J(\theta, {\nabla ^2}\theta) \!+\! \beta \frac{{\partial \psi }}{{\partial x}}\!=\!  0.5\, J\left(\psi, {f_0}\frac{h}{H}\right)\\ & \!+\! 0.5\, J\left(\theta, {f_0}\frac{h}{H}\right)  0.5\, {k_d}{\nabla ^2}(\psi  \theta), \end{split}$$ (4) $$\begin{split} \frac{\partial }{{\partial t}}({\nabla ^2}\theta) & + J(\psi, {\nabla ^2}\theta) + J(\theta, {\nabla ^2}\psi) + \beta \frac{{\partial \theta }}{{\partial x}}= 0.5\, J\left(\psi, {f_0}\frac{h}{H}\right)\\ &  0.5\, J\left(\theta, {f_0}\frac{h}{H}\right) + 0.5\, {k_d}{\nabla ^2}(\psi  \theta)  2\, {{k\,'}_{\!\!\!\!\!d}}{\nabla ^2}\theta + \frac{{{f_0}}}{{\Delta p}}\omega . \end{split} $$ (5) In the equation of temperature of the baroclinic atmosphere, a radiative and heat flux scheme is incorporated reflecting the exchanges in energy among the land, atmosphere, and space (Barsugli et al., 1998; Vannitsem et al., 2015; De Cruz et al., 2016):
$$\begin{split}{\gamma _{\rm{a}}} & \left(\frac{{\partial {T_{\rm{a}}}}}{{\partial t}} + J(\psi, {T_{\rm{a}}})  \sigma \omega \frac{p}{R}\right) =  \lambda ({T_{\rm{a}}}  {T_{\rm{g}}}) \\ & + {\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm{g}}^4  2{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm{a}}^4 + {R_{\rm{a}}},\end{split}$$ (6) where
${T_{\rm{a}}}$ and${T_{\rm{g}}}$ are atmospheric and land temperature, respectively;$\sigma $ is the static stability with$p$ as the pressure;$R$ is the gas constant for dry air;$\omega $ is the vertical velocity in pressure coordinates;${\gamma _{\rm{a}}}$ is the heat capacity of the atmosphere for a 1000hPa deep column;$\lambda $ is the heat transfer coefficient between the land and atmosphere;${\sigma _{\rm{B}}}$ is the StefanBoltzmann constant and${\varepsilon _{\rm{a}}}$ is the longwave emissivity of the atmosphere.${\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm{g}}^4$ is the longwave radiation emitted from the land that is absorbed by the atmosphere;$  2{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm{a}}^4$ is the longwave radiation emitted from the atmosphere to the land and space;${R_{\rm{a}}}$ is the shortwave solar radiation directly absorbed by the atmosphere.The land temperature equation is similar to the atmospheric temperature equation as
$${\gamma _{\rm{g}}}\frac{{\partial {T_{\rm{g}}}}}{{\partial t}} =  \lambda ({T_{\rm{g}}}  {T_{\rm{a}}})  {\sigma _{\rm{B}}}T_{\rm{g}}^4 + {\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm{a}}^4 + {R_{\rm{g}}}, $$ (7) where
${\gamma _{\rm{g}}}$ is the heat capacity of the active layer of the land for a mean thickness of 10 m (Monin, 1986);$  {\sigma _{\rm{B}}}T_{\rm{g}}^4$ is the longwave radiation emitted from the land;${\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm{a}}^4$ is the longwave radiation emitted from the atmosphere absorbed by the land;${R_{\rm{g}}}$ is the shortwave solar radiation absorbed by the land.Similar to Vannitsem et al. (2015), the quartic terms in the radiative fluxes are linearized. The details of this linearization are described in Appendix A.
The system of equations is closed by the thermal wind relation:
$$\theta = \frac{R}{{2{f_0}}}\ln \left(\frac{{{p_3}}}{{{p_1}}}\right){T_{\rm{a}}} \approx \frac{R}{{2{f_0}}}{T_{\rm{a}}}.$$ (8) Let
$x$ and$y$ be scaled by$L$ ,$t$ by$f_0^{  1}$ ,$\psi $ and$\theta $ by${L^2}{f_0}$ ,${T_{\rm{a}}}$ and${T_{\rm{g}}}$ by${{{L^2}f_0^2}/R}$ , and let$h' = {h/H}$ ,$\omega ' = {\omega /{({f_0}\Delta p}})$ ,$\beta ' = {{\beta L}/{{f_0}}}$ ,$2k = {{{k_d}}/{{f_0}}}$ ,${k^\prime } = {{{k_d}^{\!\!\!\!\prime} }/{{f_0}}}$ ,$\sigma ' = {{\sigma {{(\Delta p)}^2}}/{(2{L^2}}}f_0^2)$ .Other nondimensional coefficients are
$$\left. \begin{aligned} & {{R\,'}_{\rm{\!\!\!\!g}}} = {{{R_{\rm{g}}}R}/{({\gamma _{\rm{g}}}f_0^3{L^2})}}, \quad {R_{\rm{a}}}^{\!\!\!\!\prime} = {{{R_{\rm{a}}}R}/{(2{\gamma _{\rm{a}}}f_0^3{L^2})}} \\ & {{\lambda \,'}_{\rm{\!\!\!\!g}}} = {\lambda /{({\gamma _{\rm{g}}}}}{f_0}), \quad {{\lambda\, '}_{\!\!\!\!\rm{a}}} = {\lambda /{({\gamma _{\rm{a}}}}}{f_0}) \\ &{\sigma _{\rm B, a}} = {{8{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm a, 0}^3}/{({\gamma _{\rm{g}}}{f_0})}}, \quad {\sigma _{\rm B, g}} = {{4{\sigma _{\rm{B}}}T_{\rm g, 0}^3}/{({\gamma _{\rm{g}}}{f_0})}} \\ &{S_{\rm B, a}} = 8{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}{{T_{\rm a, 0}^3}/{({\gamma _{\rm{a}}}{f_0})}}, \quad {S_{\rm B, g}} = {{4{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{\rm g, 0}^3}/{(2{\gamma _{\rm{a}}}{f_0})}} \end{aligned} \right\}.$$ (9) We obtain the nondimensional equations of the model as
$$\begin{split} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \frac{\partial }{{\partial t'}}({\nabla ^2}\psi ') + J(\psi ', {\nabla ^2}\psi ') + J(\theta ', {\nabla ^2}\theta ') + \beta '\frac{{\partial \psi '}}{{\partial x'}} \\ =  0.5J(\psi ', h') + 0.5J(\theta ', h')  k{\nabla ^2}(\psi '  \theta '), \end{split} \quad\quad\quad\;\;\; $$ (10) $$\begin{split} & \frac{\partial }{{\partial t'}}({\nabla ^2}\theta ') + J(\psi ', {\nabla ^2}\theta ') + J(\theta ', {\nabla ^2}\psi ') + \beta '\frac{{\partial \theta '}}{{\partial x'}} \\ &\quad\quad = 0.5J(\psi ', h') \!\! 0.5J(\theta ', h')\; \!+\! k{\nabla ^2}(\psi ' \!\! \theta ') \!\! 2k'{\nabla ^2}\theta ' \!+\! \omega ' \end{split}, $$ (11) $$ \begin{split}\!\!\!\!\!\!\!\!\!\!\!\!\! & \frac{{\partial \theta '}}{{\partial t'}} + J(\psi ', \theta ')  \sigma '\omega ' =  {\lambda '_{\rm{a}}}(\theta '  0.5\delta {T_{\rm{g}}}^{\!\!\!\!\!\prime}\;) \\ &\quad\quad + {S_{{\rm B}, {\rm g}}}\delta {T_{\rm{g}}}^{\!\!\!\!\!\prime}\;  {S_{\rm B, a}}\theta ' + \delta {R'_{\rm{a}}}, \;\;\;\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \end{split}$$ (12) $$\frac{{\partial \delta {{T'}_{\!\!\!\!\rm{g}}}}}{{\partial t'}} =  {\lambda _{\rm{g}}}^{\!\!\!\!\!\prime} (\delta {T'_{\!\!\rm{g}}}  2\theta ')  {\sigma _{\rm B, g}}\delta {T'_{\rm{g}}} + {\sigma _{\rm B, a}}\theta ' + \delta {R'_{\rm{g}}}.$$ (13) Note that Eqs. (12) and (13) have been linearized. The nondimensional atmospheric temperature,
${T'_{\rm{a}}}$ , is already replaced by nondimensional$\theta '$ in Eq. (12) according to Eq. (8);$\delta {T'_{\rm{g}}}$ is the nondimensional land temperature anomaly;$\delta {R'_{\rm{a}}}$ and$\delta {R'_{\rm{g}}}$ are the nondimensional meridional differential shortwave solar radiation absorbed by the atmosphere and the land, respectively. All of the variables are now dimensionless unless otherwise specified. Hereafter, we omit the primes of the nondimensional variables$\psi '$ ,$\theta '$ ,$\delta {T'_{\rm{g}}}$ and$x'$ ,$y'$ ,$t'$ for simplicity, but others are retained to avoid confusion.We follow the work of CS, and truncate the expansions for
$\psi $ ,$\theta $ ,$\delta {T_{\rm{g}}}$ ,$\omega '$ , and$h'$ as:$$\left. \begin{aligned} & \psi = \sum\limits_{i = 1}^3 {{\psi _i}} {F_i}, \quad \theta = \sum\limits_{i = 1}^3 {{\theta _i}} {F_i}, \quad \delta {T_{\rm{g}}} = \sum\limits_{i = 1}^3 {{T_{{\rm{g}}, i}}} {F_i} \\ & \omega ' = \sum\limits_{i = 1}^3 {{\omega _i}} {F_i}, \quad h' = {h_2}{F_2} \end{aligned} \right\}.$$ (14) We choose
$$\left. \begin{aligned} & {F_1} = \sqrt 2 \cos y \\ & {F_2} = 2\cos (nx)\sin y \\ & {F_3} = 2\sin (nx)\sin y \end{aligned} \right\}.$$ (15) Here, the zonal wavenumber
$n$ may be chosen freely, and it is related to the planetary zonal wavenumber$m = {{na\cos ({\phi _0})}/{L = 2.83n}}$ , where$a$ is the radius of the earth. Note that the channel is periodic in$x$ direction over the scale${{2\pi }/n}$ .The dimensional boundary topography is given by
$$h = H{h_2}{F_2} = 2H{h_2}\cos ({{nx}/L})\sin ({y/L}).$$ (16) In our model, we set
${h_2} = 0.1$ , and thus, the dimensional amplitude of the topography is fixed at$0.2H = $ 1.46 km.The nondimensional meridional differential shortwave solar radiation absorbed by the land and the atmosphere are given by
$$\delta {R'_{\rm{g}}} = {C_{\rm{g}}}^{\!\!\!\! \prime} {F_1}, \quad \delta {R'_{\rm{a}}} = {C'_{\rm{a}}}{F_1}, $$ (17) where
${C'_{\rm{g}}} = {{{C_{\rm{g}}}R}/{({\gamma _{\rm{g}}}f_0^3{L^2})}}$ and${C_{\rm{a}}}^{\!\!\!\! \prime} = {{C_{\rm{a}}}R}/$ ${(2{\gamma _{\rm{a}}}f_0^3{L^2})}$ . The dimensional forms are$$\delta {R_{\rm{g}}} = \sqrt 2 {C_{\rm{g}}}\cos ({y/L}), \quad \delta {R_{\rm{a}}} = \sqrt 2 {C_{\rm{a}}}\cos ({y/L}), $$ (18) and we set
${C_{\rm{a}}} = 0.4{C_{\rm{g}}}$ . Thus, the dimensional meridional differences in solar heating absorbed by the land and atmosphere between the southern wall$y = 0$ and the northern wall$y = \pi L$ of the channel are$2\sqrt 2 {C_{\rm{g}}}$ and$2\sqrt 2 {C_{\rm{a}}}$ , respectively. The variable C_{g} is a dimensional parameter, which is an indicator of the meridional difference in solar heating absorbed by the land between the walls, and it is the most important varying parameter in our coupled model. Based on the observations, a typical value of C_{g} for boreal summer is 15 W m^{–2}, that for boreal winter is 55 W m^{–2}, and that for spring or autumn is 35 W m^{–2}.We can obtain 12 spectral equations and by eliminating
$\omega '$ , the number of equations can be reduced to 9. The final spectral equations are given as follows:$${\dot \psi _1} =  k({\psi _1}  {\theta _1})  c\tilde h({\theta _3}  {\psi _3}), \hspace{88pt}$$ (19) $$({n^2} + 1){\dot \psi _2} =  c{n^2}({\psi _1}{\psi _3} + {\theta _1}{\theta _3}) + \beta n{\psi _3}  {B_1}({\psi _2}  {\theta _2}), $$ (20) $$\begin{split} & ({n^2} + 1){\dot \psi _3} = c[{n^2}({\psi _1}{\psi _2} + {\theta _1}{\theta _2}) + \tilde h({\theta _1}  {\psi _1})] \\ &\quad\quad  \beta n{\psi _2}  {B_1}({\psi _3}  {\theta _3}), \end{split}\quad\quad\quad\;\;\;$$ (21) $$\begin{split} {C_1}{\dot \theta _1} =& c[{\psi _2}{\theta _3}  {\psi _3}{\theta _2}  \sigma '\tilde h({\psi _3}  {\theta _3})]  {B_3}{\theta _1}\\ & + k\sigma '{\psi _1}  {d_1}{\theta _1} + {d_2}{T_{{\rm g}, 1}} + {C_{\rm{a}}}^{\!\!\!\! \prime} \; ,\end{split}\quad\quad\quad\;\;\;$$ (22) $$\begin{split}({n^2} + 1){C_2}{\dot \theta _2} = & c({A_1}{\psi _3}{\theta _1}  {A_2}{\psi _1}{\theta _3}) + \beta n\sigma '{\theta _3}  {B_2}{\theta _2} \\ & + {B_1}\sigma '{\psi _2}  {d_1}{\theta _2} + {d_2}{T_{{\rm {\rm g}}, 2}},\end{split}\;\;\;$$ (23) $$\begin{split}({n^2} + 1){C_2}{\dot \theta _3} = & c[{A_2}{\psi _1}{\theta _2}  {A_1}{\psi _2}{\theta _1} + \sigma '\tilde h({\psi _1}  {\theta _1})] \\ & \beta n\sigma '{\theta _2}  {B_2}{\theta _3} + {B_1}\sigma '{\psi _3}  {d_1}{\theta _3} + {d_2}{T_{{\rm g}, 3}},\end{split}$$ (24) $${\dot T_{{\rm g}, 1}} =  {d_3}{T_{{\rm g}, 1}} + {d_4}{\theta _1} + {C'_{\rm{g}}}, \hspace{102pt}$$ (25) $${\dot T_{{\rm g}, 2}} =  {d_3}{T_{{\rm g}, 2}} + {d_4}{\theta _2}, \hspace{122pt}$$ (26) $$ {\dot T_{{\rm g}, 3}} =  {d_3}{T_{{\rm g}, 3}} + {d_4}{\theta _3}, \hspace{122pt}$$ (27) where the coefficients used here are
$$\left. \begin{aligned} & c = \frac{{8\sqrt 2 n}}{{3\pi }}, \quad \quad \, \tilde h = \frac{{{h_2}}}{2} \\ & {d_1} = {{\lambda '}_{\!\!\!\! \rm{a}}} + {S_{\rm B, a}}, \quad \, {d_2} = {{{{\lambda '}_{\!\!\!\! \rm{a}}}}/2} + {S_{\rm B, g}} \\ & {d_3} = {{\lambda '}_{\!\!\!\! \rm{g}}} + {\sigma _{\rm B, g}}, \quad {d_4} = 2{{\lambda '}_{\!\!\!\! \rm{g}}} + {\sigma _{\rm B, a}} \\ & {A_1} = 1  \sigma '{n^2}, \quad \;\, {A_2} = 1 + \sigma '{n^2} \\ & {B_1} = ({n^2} + 1)k, \quad {B_3} = (2k' + k)\sigma ' \\ & {B_2} = ({n^2} + 1)(2k' + k)\sigma ' \\ & {C_1} = \sigma ' + 1, \quad \;\;\;\;\;{C_2} = \sigma ' + \frac{1}{{{n^2} + 1}} \\ \end{aligned} \right\}.$$ (28)

In this section, we will show the equilibrium solutions (i.e., stationary solutions) and their stabilities of the simple coupled land–atmosphere system Eqs. (19)–(27).
We set all of the time derivatives and wave components to zero in Eqs. (19)–(27), and then obtain a specific equilibrium state:
$$\left. \begin{aligned} & {{\bar \theta }_1} = \frac{{{D_2}}}{{2k'\sigma '  {D_1}}} \\ &{{\bar \psi }_1} = {{\bar \theta }_1} \\ & {{\bar T}_{{\rm{g}}, 1}} = \frac{{{d_4}{{\bar \theta }_1} + {{C'}_{\rm{\!\!\!\!g}}}}}{{{d_3}}} \end{aligned} \right\}, $$ (29) where
$${D_1} = \frac{{{d_2}{d_4}}}{{{d_3}}}  {d_1}, \quad {D_2} = \frac{{{d_2}{{C'}_{\rm{\!\!\!\!g}}}}}{{{d_3}}} + {C'_{\rm{a}}}.$$ (30) This equilibrium state is referred to as “Hadley circulation” in CS. Note that
${\bar \psi _1} = {\bar \theta _1}$ indicates that there is no lower layer zonal flow, i.e., horizontally motionless in the lower layer, while strong westerlies without any meridional perturbations in the upper layer. Note also that the Hadley solution does not interact with the topography.The method to obtain the general equilibrium solutions of Eqs. (19)–(27) and to determine the stabilities of the equilibrium solutions are shown in Appendix B.
Next, we show the results of calculations of the equilibrium solutions and their stabilities. Similar to the “demonstration case” in RP, we preferentially choose planetary zonal wavenumber
$m = 3.7$ ($n = 1.3$ ) in this study. Wavenumber$m = 6$ is also used for comparison purposes. We set$2k = 0.02$ and$k' = 0.005$ (same as Yoden, 1983), and thus, the dimensional surface and internal frictional dissipation times are 5.6 and 22.4 days, respectively. The values of other dimensional parameters are listed in Table 1.Parameter Value Parameter Value $a$ 6371 km ${\sigma _{\rm{B}}}$ 5.6 × 10^{–8} W m^{–2} K^{–4} $\pi L$ 5000 km ${\gamma_{\rm{a}}}$ 1.0 × 10^{7} J m^{–2} K^{–1} $H$ 7.3 km ${\gamma _{\rm{g}}}$ 1.6 × 10^{7} J m^{–2} K^{–1} ${\phi _0}$ 45°N $\lambda $ 10 W m^{–2} K^{–1} ${g_0}$ 9.8 m^{–1} s^{–2} $\sigma $ 2.16 × 10^{–6} m^{2} s^{–2} Pa^{–2} ${f_0}$ 0.0001032 s^{–1} ${R_{{\rm a}, 0}}$ 89 W m^{–2} $\beta $ 1.62 × 10^{–11} m^{–1} s^{–1} ${R_{{\rm g}, 0}}$ 221 W m^{–2} $R$ 287 J kg^{–1} K^{–1} ${T_{{\rm a}, 0}}$ 270.22 K ${\varepsilon _{\rm{a}}}$ 0.76 ${T_{{\rm g}, 0}}$ 280.40 K Table 1. Dimensional parameter values used in our model

The results of the equilibrium solutions and their stabilities for
$m = 3.7$ are shown in Table 2. We will focus on the stable equilibrium states, and the unstable equilibrium states are rarely described.C_{g} (W m^{–2}) ${\psi _1}$ ${\psi _2}$ ${\psi _3}$ ${\theta _1}$ ${\theta _2}$ ${\theta _3}$ ${T_{{\rm g}, 1}}$ ${T_{{\rm g}, 2}}$ ${T_{{\rm g}, 3}}$ $S$ Character 20 0.0258 0 0 0.0258 0 0 0.0603 0 0 － Hadley 30 0.0387 0 0 0.0387 0 0 0.0905 0 0 － Hadley 40 0.0516 0 0 0.0516 0 0 0.1207 0 0 － Hadley 45 0.0580 0 0 0.0580 0 0 0.1358 0 0 － Hadley 50 0.0644 0 0 0.0644 0 0 0.1509 0 0 N 0.0620 –0.0025 –0.0071 0.0633 –0.0030 –0.0069 0.1487 –0.0053 –0.0123 － High 2 0.0633 –0.0035 0.0111 0.0618 –0.0026 0.0109 0.1461 –0.0047 0.0194 － High 1 55 0.0709 0 0 0.0709 0 0 0.1659 0 0 N 0.0620 –0.0115 –0.0092 0.0661 –0.0126 –0.0087 0.1575 –0.0225 –0.0155 － High 2 0.0643 –0.0124 0.0187 0.0612 –0.0104 0.0183 0.1487 –0.0186 0.0328 － Low 1 60 0.0773 0 0 0.0773 0 0 0.1810 0 0 N 0.0804 0.0227 0.0012 0.0711 0.0177 0 0.1699 0.0316 0 N 0.0650 –0.0201 0.0219 0.0610 –0.0171 0.0214 0.1517 –0.0305 0.0383 － Low 1 70 0.0902 0 0 0.0902 0 0 0.2112 0 0 N 0.0816 0.0437 0 0.0685 0.0324 –0.0017 0.1724 0.0579 –0.0031 N 0.0664 –0.0329 0.0246 0.0607 –0.0279 0.0238 0.1583 –0.0498 0.0426 － Low 1 80 0.1031 0 0 0.1031 0 0 0.2414 0 0 N 0.0816 0.0563 –0.0013 0.0674 0.0412 –0.0031 0.1776 0.0736 –0.0055 N 0.0676 –0.0433 0.0257 0.0605 –0.0362 0.0248 0.1652 –0.0648 0.0443 － Low 1 The column $S$ describes the stability of the equilibria. Sign “－” denotes a stable equilibrium, and “N” denotes an unstable equilibrium. Table 2. Nondimensional equilibrium solutions for m = 3.7
In Table 2, for a given realistic value of C_{g}, there may be one (C_{g}
$ \leqslant $ 45 W m^{–2}) or three (C_{g}$ > $ 45 W m^{–2}) equilibrium states, and some of them are stable. Some of the stable states are highindex equilibria, others are lowindex equilibria. The criteria of lowindex equilibria are that there is at least one closed streamline for low or high pressure center in the upper layer and in addition that the magnitude of streamfunction must be no less than 10^{7} m^{2} s^{–1} in the upper layer and 10^{6} m^{2} s^{–1} in the lower layer. Those that do not meet the above criteria belong to highindex equilibria. The Hadley solution is a specific highindex equilibria.For C_{g} = 20 W m^{–2}, the only one equilibrium state is stable (see Table 2). All of the wave components of this equilibrium state are zero, so this is Hadley solution. For C_{g} = 30, 40, and
$45$ W m^{–2}, the results are the same as for C_{g} = 20 W m^{–2}.For C_{g} = 50 W m^{–2}, there are three equilibrium states and only the last two are stable. Note that the first equilibrium state is the Hadley solution, and now it becomes unstable. The streamfunction and temperature fields of the second and third equilibrium states are illustrated in Fig. 1.
Figure 1. The second one (left panels) and third one (right panels) of the three equilibrium states for m = 3.7 at C_{g} = 50 W m^{–2}. They belong to “High 2” and “High 1” equilibria, respectively. The streamfunction fields of the (a, e) upper and (b, f) lower layers, respectively. The temperature fields of (c, g) the atmosphere and (d, h) the land, respectively. The contour intervals are (a, e) 2.0 × 10^{7} m^{2} s^{–1}, (b) 2.0 × 10^{5} m^{2} s^{–1}, (f) 3.0 × 10^{5} m^{2} s^{–1}, and (c, d, g, h) 10 K. The background dotted lines show the topographic heights in the model, with negative regions shaded.
The second and third equilibrium states are both highindex, due to both strong zonal westerlies with weak meridional perturbations in the upper layer (Figs. 1a, e). However, there are wavy easterlies in the lower layer for the second equilibrium state (Fig. 1b) and wavy westerlies in the lower layer for the third equilibrium state (Fig. 1f). The isotherms in both atmospheric and land temperature fields of the two equilibrium states are all quite flat (Figs. 1c, d, g, h) and almost in phase with each upper layer streamfunction field. Both of the two equilibrium states have a characteristic baroclinic structure, i.e., the waves of streamfunction fields displayed westward phase shifts with height. However, they have different wave phases relative to the topography. For the second equilibrium state, its lower layer streamfunction is nearly in phase with the topography, the lower layer converse ridges (anticyclonic flow) are over the mountains (positive topographic heights), lying slightly west of the mountain crests (Fig. 1b), and the upper layer ridges are located to the west side of the mountains (Fig. 1a). We call this a “ridgetype” equilibrium. By contrast, for the third equilibrium state, the lower layer streamfunction is nearly out of phase with the topography, the lower layer lowpressure centers and troughs are over the mountains, also lying slightly west of the mountain crests (Fig. 1f), and the upper layer troughs are located to the west side of the mountains (Fig. 1e). We call this a “troughtype” equilibrium. For simplicity, we refer to the characters of the two equilibrium states as “High 2” and “High 1”, respectively. Here, “High” denotes “highindex”, “2” denotes “ridgetype”, and “1” denotes “troughtype”.
For C_{g} = 55 W m^{–2}, the last two of the three equilibrium states are also stable. The first one is still “High 2” highindex equilibrium, and the second one becomes lowindex equilibrium (Table 2). The streamfunction and temperature fields of this lowindex equilibrium are illustrated in Fig. 2 (left panel). There are relatively weak westerlies with strong meridional flow in both the upper and lower layer streamfunction fields (Figs. 2a, b), particularly closed streamlines in the former. Note that the magnitude of the streamfunction in Fig. 2b is 10^{6} m^{2} s^{–1} and larger than that in Fig. 1b (10^{5} m^{2} s^{–1}), indicating that the amplitude of meridional perturbations in Fig. 2b are larger than those in Fig. 1b. There are also relatively large meridional perturbations in both the atmospheric and land temperature fields (Figs. 2c, d) and even closed isotherms in the latter. In this lowindex equilibrium, the lower layer streamfunction is nearly out of phase with the topography, the lower layer troughs are over the mountains and the upper layer troughs are located on the west side of the mountains, so this is a troughtype equilibrium. We refer to the character of this equilibrium state as “Low 1”, where “Low” denotes “lowindex”.
Figure 2. As in Fig. 1, but for the third one of the three equilibrium states for m = 3.7 at C_{g} = 55 W m^{–2} (left panels) and for the third one of the five equilibrium states for m = 6 at C_{g} = 30 W m^{–2} (right panels). They belong to “Low 1” and “Low 2” equilibria, respectively. The “Low 2” equilibria has nondimensional solutions with
$({\psi _1}, \;{\psi _2}, \;{\psi _3}, \;{\theta _1}, \;{\theta _2}, \;{\theta _3}, \;{T_{{\rm g}, 1}}, \;{T_{{\rm g}, 2}}, \;{T_{{\rm g}, 3}})$ = (0.0216, –0.0079, –0.0051, 0.0311, –0.0098, –0.0043, 0.0771, –0.0175, –0.0077). The contour intervals are (a) 2.0 × 10^{7} m^{2} s^{–1}, (e) 1.0 × 10^{7} m^{2} s^{–1}, (b, f) 1.0 × 10^{6} m^{2} s^{–1}, (c, d) 10 K, and (g, h) 4 K.At C_{g} = 60 W m^{–2}, only the third one of the three equilibrium states is stable, and it is “Low 1” equilibrium. The second one becomes unstable. For C_{g} = 70 and 80 W m^{–2}, the results are the same as for C_{g} = 60 W m^{–2}.
For comparison purposes, we have calculated the equilibrium solutions for wavenumber 6. We find there may be one, three, or five equilibrium states for a given value of C_{g} (figure omitted). Some of them are stable. Besides the stable “High 1”, “High 2”, and “Low 1” equilibrium, a new stable lowindex equilibrium may exist. As illustrated in Fig. 2 (right panel), it has strong meridional perturbations in both the upper layer streamfunction field (Fig. 2e) and temperature fields (Figs. 2g, h). However, there are wavy easterlies in the lower layer streamfunction field (Fig. 2f). Note that the lower layer streamfunction is nearly in phase with the topography, the lower layer converse ridges are over the mountains, and the upper layer ridges are located to the west side of the mountains, so this is a ridgetype equilibrium. We refer to the character of this equilibrium state as “Low 2”.

To further demonstrate the multiple equilibrium states and their stabilities for wavenumbers 3.7 and 6, simple bifurcation diagrams are shown in Fig. 3. The zonal component
$\psi _1^1$ , the wave component$\psi _2^1$ , and$\psi _3^1$ of the upper layer streamfunction${\psi ^1}$ are given by$\psi _1^1 = {\psi _1} + {\theta _1}$ ,$\psi _2^1 = {\psi _2} + {\theta _2}$ , and$\psi _3^1 = {\psi _3} + {\theta _3}$ , respectively. The equilibrium solutions are shown by the 2W m^{–2} interval of the parameter C_{g}.Figure 3. The equilibrium bifurcation associated with the change in meridional differential solar heating parameter C_{g} for m = 3.7 (left panels) and m = 6 (right panels), respectively. The ordinate shows the nondimensional equilibrium values of (a, d) the zonal component
$\psi _1^1$ and the wave components (b, e)$\psi _2^1$ and (c, f)$\psi _3^1$ , respectively. Different branches of equilibrium solutions have different colors. The crosses denote unstable equilibria, the circles denote stable highindex equilibria, and the asterisks denote stable lowindex equilibria.For wavenumber 3.7, there are four equilibrium branches (Fig. 3, left panel). For small values of C_{g}, the Hadley circulation (black) is the only equilibrium and it is stable. As C_{g} is gradually increased, around C_{g} = 50 W m^{–2}, the Hadley circulation loses its stability, and two new equilibria (blue and red) appear. The blue branch represents a troughtype equilibrium, and it is always stable. It includes “High 1” (50
$ \leqslant {C_{\rm{g}}} \leqslant $ 52 W m^{–2}) and “Low 1” equilibrium (${C_{\rm{g}}} \geqslant $ 54 W m^{–2}). The red branch represents a ridgetype equilibrium, and it includes stable “High 2” equilibrium (50$ \leqslant {C_{\rm{g}}} \leqslant $ 54 W m^{–2}). It becomes unstable around C_{g} = 56 W m^{–2}, then it disappears and a new equilibria (green) appears when${C_{\rm{g}}} > $ 56 W m^{–2}. This green branch equilibrium is always unstable.For wavenumber 6, there are five equilibrium branches (Fig. 3, right panel). For small values of C_{g}, the stable Hadley circulation (black) is still the only equilibrium. As C_{g} is increased to around C_{g} = 20 W m^{–2}, the Hadley circulation becomes unstable, and two new equilibria (blue and red) appear. The blue branch represents a troughtype equilibrium and it is always stable. It includes “High 1” (C_{g} = 20 W m^{–2}) and “Low 1” equilibrium (
${C_{\rm{g}}} \geqslant $ 22 W m^{–2}). The red branch represents a ridgetype equilibrium, and it includes stable “High 2” equilibrium (20$ \leqslant {C_{\rm{g}}} \leqslant $ 24 W m^{–2}) and stable “Low 2” equilibrium (26$ \leqslant {C_{\rm{g}}} \leqslant $ 30 W m^{–2}). The ridgetype equilibrium is unstable within 32$ \leqslant {C_{\rm{g}}} \leqslant $ 36 W m^{–2} and it disappears when${C_{\rm{g}}} > $ 36 W m^{–2}. At around C_{g} = 26 W m^{–2}, two more equilibria (green and magenta) appear, and they are both always unstable. The magenta branch disappears when${C_{\rm{g}}} > $ 36 W m^{–2}.The above results indicate that there are multiple equilibrium states with different wave phases and wave amplitudes in the coupled model. For a considerable range of C_{g} values (50
$ \leqslant {C_{\rm{g}}} \leqslant $ 54 W m^{–2} for$m = 3.7$ and 20$ \leqslant {C_{\rm{g}}} \leqslant $ 30 W m^{–2} for$m = 6$ ), two stable equilibria with distinct wave phase relative to the topography, i.e., ridge and troughtype equilibria, may simultaneously exist (Fig. 3). However, only for a small range of C_{g} values (C_{g} = 54 W m^{–2} for$m = 3.7$ and 22$ \leqslant {C_{\rm{g}}} \leqslant $ 24 W m^{–2} for$m = 6$ ), two stable equilibria with distinct wave amplitude, i.e., high and lowindex equilibria, may coexist (Fig. 3). Therefore, the multiple wave phase equilibria associated with ridge and troughtypes are more prominent than the multiple wave amplitude equilibria associated with high and lowindex types. 
The multiple wavelike stationary equilibrium states exist in the model when the topography is present. This is proved in Appendix B.
Figure 5a shows the stability curves of the Hadley circulation in the coupled model. The blue lines enclose the orographically unstable region. In crossing the blue lines from the stable to unstable sides, the variable
$\alpha $ (see Appendix B) changes from a negative real value to a positive real value (pitchfork bifurcation). The red (black dashed) lines separate the baroclinically stable and unstable regions in the presence (absence) of topography. In crossing these lines from the stable to unstable sides, the real part of the complex$\alpha $ changes from negative to positive while the imaginary part is not zero (Hopf bifurcation). The orographic instability of the Hadley circulation is only present when the topography is present. Moreover, there is no overlap between the orographic instability and baroclinic instability in the presence of topography. Besides, just compared the baroclinic stability curves with/without topography, it is seen that the presence of topography stabilizes the Hadley circulation for most wavenumbers.Figure 5. (a) Stability curves of the Hadley circulation in the coupled land–atmosphere model (Case 1). The blue solid lines enclose the region of orographic instability. The red solid (black dashed) lines and the top xaxis and the right yaxis enclose the region of baroclinic instability in the presence (absence) of topography. Comparison of the regions of (b) orographic instability, (c) baroclinic instability in the presence of topography, and (d) baroclinic instability in the absence of topography for the four experiments (Cases 1–4).
In the absence of topography, there is only the Hadley circulation (see Appendix B) or the traveling wave due to the baroclinic instability of the Hadley circulation (which is a Hopf bifurcation). For example, for
$m = 3.7$ at C_{g} = 50 W m^{–2} without the topography, numerical integration starting at arbitrary initial conditions converges to a periodic solution of period 18 days (Fig. 4). Here, the zonal component$\psi _1^3$ and the wave component$\psi _3^3$ of the lower layer streamfunction${\psi ^3}$ are given by$\psi _1^3 = {\psi _1}  {\theta _1}$ and$\psi _3^3 = {\psi _3}  {\theta _3}$ , respectively. In this periodic solution, there is blockinglike flow in the upper layer streamfunction (Fig. 4a); however, there is no zonal flow (the zonal component$\psi _1^3$ remains zero in Fig. 4c) but wave train in the lower layer streamfunction (Fig. 4b). It is seen that the zonal components of upper and lower layer streamfunction ($\psi _1^1$ ,$\psi _1^3$ ) remain constant (Fig. 4c), whereas the wave components of upper and lower layer streamfunction ($\psi _3^1$ ,$\psi _3^3$ ) evolve periodically with time (Fig. 4d). Particularly, this traveling wave moves westward.Figure 4. A traveling wave solution in the absence of topography for m = 3.7 at C_{g} = 50 W m^{–2}, with
$({\psi _1}, \;{\psi _2}, \;{\psi _3}, \;{\theta _1}, \;{\theta _2}, \;{\theta _3}, \;{T_{{\rm g}, 1}}, \;{T_{{\rm g}, 2}}, \;{T_{{\rm g}, 3}})$ = (0.0534, 0.0214, 0.0016, 0.0534, 0.0125, –0.0008, 0.1311, 0.0014, 0.0048) at this moment. (a, b) The streamfunction fields of the upper and lower layers, respectively. The contour intervals are (a) 2.0 × 10^{7} m^{2} s^{–1} and (b) 2.0 × 10^{6} m^{2} s^{–1}. (c) Temporal evolution for the nondimensional equilibrium values of the zonal component$\psi _1^1$ (solid line) and$\psi _1^3$ (dashed line). (d) Temporal evolution for the nondimensional equilibrium values of the wave component$\psi _3^1$ (solid line) and$\psi _3^3$ (dashed line).In the presence of topography, there may exist multiple equilibrium states. Compared Fig. 1 with Fig. 4, it seems that due to the presence of topography, the traveling wave becomes two types of stationary waves. In fact, the first bifurcation (Fig. 3, around C_{g} = 50 W m^{–2} for
$m = 3.7$ and around C_{g} = 20 W m^{–2} for$m = 6$ ) results from the orographic instability of the Hadley circulation (Fig. 5a, around C_{g} = 50 W m^{–2} for$m = 3.7$ and around C_{g} = 20 W m^{–2} for$m = 6$ ), and it is a (supercritical) pitchfork bifurcation. This bifurcation is important, because it determines the occurrence and coexistence of the trough and ridgetype equilibria. Note that the disappearance of the ridgetype equilibria (Fig. 3, around C_{g} = 56 W m^{–2} for$m = 3.7$ and around C_{g} = 36 W m^{–2} for$m = 6$ ) is not related to the occurrence of baroclinic instability of the Hadley circulation (Fig. 5a, around C_{g} = 64 W m^{–2} for$m = 3.7$ and around C_{g} = 28 W m^{–2} for$m = 6$ ).Therefore, the multiple wave phase equilibria associated with the ridge and troughtypes originate from the orographic instability of the Hadley circulation, which is a pitchfork bifurcation.
3.1. Equilibrium solutions
3.2. Bifurcation diagrams
3.3. The origin of the multiple equilibria

In this section, we explore the role of the land–atmosphere coupling in the existence and properties of the equilibrium states.
Four experiments are designed with different diabatic heating terms (see Table 3). For simplicity, here we refer to the coupled land–atmosphere model as Case 1. Equations (6) and (7) indicate that the diabatic heating terms in Case 1 include three terms: heat flux, longwave radiation, and shortwave radiation. For Case 2, we replace all the terms on the right side of Eq. (6) by specified heating
${Q^*}$ (${Q^*} = \sqrt 2 Q\cos ({y/L})$ , where$Q$ is a heating parameter that is similar to C_{g}) and delete the Eq. (7). Thus, Case 2 is just the classic uncoupled model. For Case 3, we delete the heat flux terms (the first terms on the right side) of both Eqs. (6) and (7) (or set$\lambda = 0$ W m^{–2} technically). For Case 4, we delete the longwave radiation terms (the second and third terms on the right side) of both Eqs. (6) and (7) (or set${\sigma _{\rm{B}}} = 0$ W m^{–2} K^{–4} technically).Experiment Heat flux Longwave radiation Shortwave radiation Specified heating Case 1 √ √ √ Case 2 √ Case 3 √ √ Case 4 √ √ The sign (√) indicates that the corresponding element has been included in the related experiment. Table 3. Experiment design for examination of the role of land–atmosphere coupling

Figure 5b compares the orographic instability of the Hadley circulation in the four experiments. Clearly, compared with Case 1, the thresholds of orographic instability in Cases 2 and 4 are both greatly reduced for wavenumbers 1–6. Moreover, the orographically unstable regions in Cases 2 and 4 are both very narrow. Unexpectedly, the orographically unstable regions in Cases 3 and 1 almost completely overlap. Figures 5c and 5d compare the baroclinic instability of the Hadley circulation with and without topography, respectively. Similarly, compared with Case 1, the thresholds of baroclinic instability in Cases 2 and 4 are both greatly reduced for wavenumbers 1–8. However, the baroclinic stability curves in Cases 3 and 1 roughly overlap. The results indicate that compared with the uncoupled model (Case 2), the land–atmosphere coupling (Case 1) greatly stabilizes the Hadley circulation, and this stabilizing effect is primarily attributed to the presence of longwave radiation fluxes, but not the heat fluxes.
In addition, compared with Case 2, Case 4 has lower thresholds for all of the orographic instability and baroclinic instability with and without topography (Figs. 5b–d). It suggests that the presence of heat fluxes extremely destabilizes the Hadley circulation, no matter with or without topography. Nevertheless, in Case 1, which presents both the heat fluxes and the longwave radiation fluxes, the destabilizing effect of the heat fluxes on Hadley circulation is nearly entirely suppressed.

Next, we compare the equilibrium bifurcation in the four experiments. Figures 6a and 6b show the bifurcation diagrams in Case 3 for wavenumber
$m = 3.7$ and$m = 6$ , respectively. The equilibrium solutions are shown by the 2W m^{–2} interval of C_{g}. Even though Cases 3 and 1 have almost overlapped orographic and baroclinic instability curves (Fig. 5), they still have nonnegligible difference in the equilibrium bifurcation. Compared with Fig. 3c, it is seen that almost all the ridgetype equilibrium states become unstable for$m = 3.7$ in Case 3 in which the heat flux is absent (Fig. 6a). However, the equilibrium bifurcation has little change for$m = 6$ in Case 3 (Fig. 6b). These suggest that the presence (absence) of heat fluxes more or less stabilizes (destabilizes) the ridgetype equilibrium.Figure 6. As in Fig. 3, but for Case 3 (without heat flux) with (a) m = 3.7 and (b) m = 6, (c) for Case 2 (without coupling) with m = 3.7, and (d) for Case 4 (without longwave radiation) with m = 3.7. Each ordinate shows the nondimensional equilibrium solution of the wave component
$\psi _3^1$ .Figures 6c and 6d show the bifurcation diagrams in Cases 2 and 4 for wavenumber
$m = 3.7$ , respectively. The equilibrium solutions are shown by the 0.5W m^{–2} interval of$Q$ or C_{g}. Due to the low thresholds of orographic instability of the Hadley circulation in Cases 2 and 4 (Fig. 5b), the first bifurcation of the Hadley circulation occurs at considerably small heating parameter values (around$Q = 12.5$ W m^{–2} in Fig. 6c and around C_{g} = 9 W m^{–2} in Fig. 6d). However, the ridgetype equilibrium subsequently disappears at small heating parameter values (around$Q = 14.0$ W m^{–2} in Fig. 6c and around C_{g} = 10 W m^{–2} in Fig. 6d), mainly because the orographically unstable regions in Cases 2 and 4 are very narrow (Fig. 5b). In this case, only for a very small range of heating parameter values, the ridge and troughtype equilibria may coexist (Figs. 6c, d). By contrast, for a considerable range of heating parameter values, the ridge and troughtype equilibria may coexist in Case 1 (Fig. 3), mainly due to the fairly wide orographically unstable region (Figs. 5a, b). Therefore, compared with the uncoupled model (Case 2), the multiple wave phase equilibria associated with the ridge and troughtypes in the coupled model (Case 1) is more remarkable. 
Here, we compare the streamfunction and temperature fields of equilibrium states in the four experiments for
$m = 3.7$ at the same heating parameter values:$Q = 50$ or C_{g} = 50 W m^{–2}. In Case 2, there is only one stable equilibrium state (Fig. 6c), with blockinglike large amplitude perturbations in both streamfunction and temperature fields (Fig. 7, left panel). Compared with the two equilibrium states in Case 1 (Fig. 1), it is obvious that the meridional perturbations in streamfunction and temperature fields of the equilibrium state in Case 2 are much stronger. To some extent, this result is attributed to the very low threshold of orographic instability of the Hadley circulation in Case 2 (Fig. 5b). In Case 3, the streamfunction and temperature fields of the two stable equilibrium states (Fig. 8) are very similar to those in Case 1 (Fig. 1), while the meridional perturbations of the lower layer streamfunction (Figs. 8b, f) are apparently weaker than those in Case 1 (Figs. 1b, f), whereas the meridional gradients of land temperature (Figs. 8d, h) are moderately greater than those in Case 1 (Figs. 1d, h). In Case 4, the result is similar to Case 2, but the meridional perturbations in streamfunction and temperature fields (Fig. 7, right panel) are stronger than those in Case 2 (Fig. 7, left panel), probably due to the lower threshold of orographic instability of the Hadley circulation in Case 4 than that in Case 2 (Fig. 5b)Figure 7. As in Fig. 1, but for the only stable equilibrium state for m = 3.7 at Q = 50 W m^{–2} in Case 2 (left panels) and for the third one of the three equilibrium states for m = 3.7 at C_{g} = 50 W m^{–2} in Case 4 (right panels). The former has nondimensional solutions with
$({\psi _1}, \;{\psi _2}, \;{\psi _3}, \;{\theta _1}, \;{\theta _2}, \;{\theta _3})$ = (0.0656, –0.0679, 0.0259, 0.0613, –0.0616, 0.0253). The latter has nondimensional solutions with$({\psi _1}, \;{\psi _2}, \;{\psi _3}, \;{\theta _1}, \;{\theta _2}, \;{\theta _3}, \;{T_{{\rm g}, 1}}, \;{T_{{\rm g}, 2}}, \;{T_{{\rm g}, 3}})$ = (0.0662, –0.0863, 0.0278, 0.0615, –0.0776, 0.0272, 0.1762, –0.1551, 0.0544). They both belong to “Low 1” equilibria. The contour intervals are (a, e) 4.0 × 10^{7} m^{2} s^{–1}, (b, f) 2.0 × 10^{6} m^{2} s^{–1}, and (c, g, h) 10 K. Note that there is no land temperature field in Case 2 (left panel).Figure 8. As in Fig. 1, but for the second one (left panels) and third one (right panels) of the three equilibrium states for m = 3.7 at C_{g} = 50 W m^{–2} in Case 3. They have nondimensional solutions with
$({\psi _1}, \;{\psi _2}, \;{\psi _3}, \;{\theta _1}, \;{\theta _2}, \;{\theta _3}, \;{T_{{\rm g}, 1}}, \;{T_{{\rm g}, 2}}, \;{T_{{\rm g}, 3}})$ = (0.0623, –0.0002, –0.0020, 0.0626, –0.0004, –0.0019, 0.1596, –0.0004, –0.0021) and (0.0628, –0.0006, 0.0052, 0.0621, –0.0003, 0.0052, 0.1591, –0.0003, 0.0056), respectively. They belong to “High 2” and “High 1” equilibria, respectively. The contour intervals are (a, e) 2.0 × 10^{7} m^{2} s^{–1}, (b, f) 1.0 × 10^{5} m^{2} s^{–1}, and (c, d, g, h) 10 K.These results indicate that compared with the uncoupled model (Case 2), the land–atmosphere coupling may weaken the atmospheric response to the thermal and topographic forcing, and this weakening effect is mainly contributed by the presence of longwave radiation fluxes. The presence of heat fluxes greatly strengthens the atmospheric response to the thermal and topographic forcing, but in the coupled model which combined the heat fluxes and longwave radiation fluxes, the heat fluxes just strengthen the response of the lower layer flow, and moderately reduce the meridional gradient of the land temperature.

To further understand the reason of the different results in the four experiments, we should compare the heating fields in the four experiments.
Figure 9 demonstrates the heating fields of the “High 2” and “High 1” equilibrium states shown in Fig. 1 (left and right panels), respectively. The zonally symmetric shortwave radiation fields for the two equilibrium states are identical (Figs. 9a, e). The isolines in all of the longwave radiation fields, the heat flux fields, and the net diabatic heating fields are wavelike, while the wave phases relative to the topography are different. For the “High 2” equilibrium state (Fig. 1, left panel), the “heating ridges” are located on the east side of the mountains (Fig. 9b–d); By contrary, for the “High 1” equilibrium state (Fig. 1, right panel), the “heating ridges” are located on the west side of the mountains (Figs. 9f–h). Note that for the lower layer streamfunction of the two equilibrium states, the ridges (high pressure) are always generated on west side of the “heating ridge”, and the troughs (low pressure) are always generated on east side of the “heat ridges”. Furthermore, it is noteworthy that the longwave radiation fluxes increase from low to high latitudes (Figs. 9b, f); thus, the presence of longwave radiation fluxes reduce the meridional gradient of the net diabatic heating field, resulting in a more stable atmosphere flow. On the contrary, the heat fluxes decrease from low to high latitudes (Figs. 9c, g); thus, the presence of heat fluxes increase the meridional gradient of the net diabatic heating field, resulting in a less stable atmosphere flow.
Figure 9. The heating fields of the “High 2” (left panels) and “High 1” (right panels) equilibrium states shown in Fig. 1, respectively. (a, e) The shortwave radiation, (b, f) the longwave radiation, (c, g) the heat flux, and (d, h) the net diabatic heating absorbed by the atmosphere. All of the contour intervals are 10 W m^{–2}. The background dotted lines show the topographic heights in the model, with negative regions shaded.
The net diabatic heating field in Case 2 is zonally symmetric (Fig. 10a). Particularly, the meridional gradient of the net diabatic heating is much greater than that in Case 1 (Figs. 9d, h). The net diabatic heating fields for the “High 2” and “High 1” equilibrium states in Case 3 (Figs. 10c, d) are similar to those in Case 1 (Figs. 9d, h), while the meridional gradients of the net diabatic heating are smaller than those in Case 1. The net diabatic heating field in Case 4 (Fig. 10b) is almost the same as that in Case 2 (Fig. 10a); however, the meridional gradient of the net diabatic heating is greater than that in Case 2. It suggests that compared with the uncoupled model (Case 2), the land–atmosphere coupling reduces the meridional gradient of the net diabatic heating, and this effect is mainly attributed to the presence of longwave radiation fluxes. The presence of heat fluxes greatly increase the meridional gradient of the net diabatic heating. However, in the coupled model that combines the heat fluxes and longwave radiation fluxes, the heat fluxes just moderately increase the meridional gradient of the net diabatic heating.
Figure 10. (a)–(d) The net diabatic heating absorbed by the atmosphere for the equilibrium states shown in Figs. 7 and 8, respectively. The contour intervals are (a, b) 30 W m^{–2} and (c, d) 5 W m^{–2}. The background dotted lines show the topographic heights in the model, with negative regions shaded.
To sum up, compared with the uncoupled model, the multiple wave phase equilibria associated with the ridge and troughtypes in the coupled model is more remarkable, mainly because the land–atmosphere coupling expands the region of orographic instability of the Hadley circulation. Besides, the land–atmosphere coupling greatly stabilizes the Hadley circulation and weakens the atmospheric response to the thermal and topographic forcing. Particularly, these effects of the land–atmosphere coupling are primarily attributed to the presence of longwave radiation fluxes, which increase from low to high latitudes, reducing the meridional gradient of the net diabatic heating. The presence of heat fluxes more or less modify the effects of longwave radiation fluxes.
4.1. Comparing the stability of the Hadley circulation
4.2. Comparing the bifurcation
4.3. Comparing the streamfunction and temperature fields
4.4. Comparing the heating fields

Next, we investigate the wave phases of ridge and troughtype equilibria relative to the topography. It is clear that the wave components
$\psi _2^1$ of the ridge and troughtype equilibria are both negative (Figs. 3b, e). The negative sign of$\psi _2^1$ denotes that this wave component of upper layer streamfunction of the two types of equilibrium is out of phase with the topography. The wave components$\psi _3^1$ of the ridge and troughtype equilibria are negative and positive, respectively (Figs. 3c, f). The negative (positive) sign of$\psi _3^1$ denotes that this wave component of upper layer streamfunction of the ridgetype (troughtype) equilibria has a lag (lead) in phase by 90° relative to the mountain crests. Therefore, the upper layer ridges (troughs) of ridgetype (troughtype) equilibria are located to the west side of the mountains.We have calculated the wave phase of the streamfunction relative to the mountains for wavenumbers 3.7 and 6, and the results are shown in Table 5. The ridgetype (High 2 and Low 2) equilibrium states have lower layer ridges over the mountains, their upper layer ridges are located to the west side of the mountain crests, and they have lower layer easterlies (Mean_U^{3} is negative, see Table 4 for definition of Mean_U^{3}). The troughtype (High 1 and Low 1) equilibrium states have lower layer troughs over the mountains, their upper layer troughs are located to the west side of the mountain crests, and they have lower layer westerlies (Mean_U^{3} is positive). Figures 11b and 12b also show that the ridgetype (troughtype) equilibria has lower layer easterlies (westerlies). Therefore, the distinct characters of ridge and troughtype equilibria are robust.
Notation Variable U^{1} (m s^{–1}) Zonal mean upperlayer ucomponent (east–west) wind U^{3} (m s^{–1}) Zonal mean lowerlayer ucomponent wind Mean_U^{1} (m s^{–1}) Channelaverage upperlayer ucomponent wind, i.e., mean of U^{1} Mean_U^{3} (m s^{–1}) Channelaverage lowerlayer ucomponent wind, i.e., mean of U^{3} Mean_U^{2} (m s^{–1}) Middlelevel ucomponent wind, mean of Mean_U^{1} and Mean_U^{3} AH (gpm) The amplitude of wave components of upperlayer geopotential height field ΔT_{a} (K) Meridional gradient of atmospheric temperature, mean atmospheric temperature at the southern wall minus that at
the northern wallΔT_{g} (K) Meridional gradient of land temperature, mean land temperature at the southern wall minus that at the northern wall AT_{a} (K) The amplitude of wave components of atmospheric temperature field AT_{g} (K) The amplitude of wave components of land temperature field Table 4. List of variables and notations used in this study
m C_{g} (W m^{–2}) Character Phase ΔPhase (°) Mean_U^{3} (m s^{–1}) g_{1} ($ \times {10^{  11}}$m^{–2}) g_{2} ($ \times {10^{  11}}$m^{–2}) Lower Upper 50 High 2 Ridge –12 –84 –0.18 9.10 –1.49 50 High 1 Trough –9 –57 0.23 –6.93 1.17 3.7 55 High 2 Ridge –18 –108 –0.61 2.76 –0.44 55 Low 1 Trough –9 –45 0.44 –3.57 0.61 60 Low 1 Trough –6 –36 0.60 –2.59 0.45 80 Low 1 Trough –6 –24 1.03 –1.46 0.26 20 High 2 Ridge 0 –42 –0.20 8.31 –1.68 20 High 1 Trough 0 –24 0.26 –6.01 1.29 26 Low 2 Ridge –6 –60 –0.77 2.32 –0.44 6 26 Low 1 Trough 0 –12 0.49 –3.09 0.69 28 Low 2 Ridge –6 –66 –1.06 1.74 –0.32 28 Low 1 Trough 0 –12 0.52 –2.90 0.65 32 Low 1 Trough 0 –6 0.55 –2.73 0.61 40 Low 1 Trough 0 –6 0.56 –2.67 0.60 “Ridge (trough)” in the “Phase” column indicates that the ridges (troughs) of lower layer streamfunction are over the mountains. The subsequent column “ΔPhase” gives the phase of the ridges (troughs) of the lower layer and upper layer streamfunction relative to the mountain crests, respectively, and negative values indicate that the ridges or troughs are located to the west side of the mountain crests. Table 5. Wave phase of the equilibrium states relative to the mountains
Figure 11. As in Fig. 3, but for dimensional variables for m = 3.7. Only stable equilibrium states are shown here. The blue (red) branch represents troughtype (ridgetype) equilibria, and the black branch represents the Hadley equilibria.
Figure 12. As in Fig. 11, but for m = 6.
The above phenomena can roughly be explained by the forced topographic Rossby wave theory. The forced topographic Rossby wave solution based on the barotropic potential vorticity equation (Smith, 1979; Nigam and DeWeaver, 2003; Holton and Hakim, 2012) is given by
$$\Psi (x, y) = {\mathop{\rm Re}\nolimits} \left[\frac{{{f_0}h/{H_0}}}{{{{\tilde k}^2} + {{\tilde l}^2}  \beta /\bar u  i\varepsilon ({{\tilde k}^2} + {{\tilde l}^2})/(\bar u\tilde k)}}\right], $$ (31) where
$\operatorname{Re} []$ denotes the real part;$\tilde k$ and$\tilde l$ are zonal and meridional wavenumbers, respectively;$h$ is the boundary topography;${H_0}$ is the height of the homogeneous atmosphere;$\bar u$ is the mean zonal wind speed;$\varepsilon $ is the dissipation factor. In our model,$\tilde k = {n/L}$ ,$\tilde l = {1/L}$ ,$h = 2H{h_2}\cos (nx/L)\sin (y/L)$ ,$\varepsilon = {k_d}$ and we might set${H_0} = H$ . The boundary topography has little effect on the upper layer flow; therefore, we choose the mean lower layer zonal wind speed, i.e., Mean_U^{3}, as$\bar u$ .We might write the boundary topography as
$$h = {\mathop{\rm Re}\nolimits} \{ 2H{h_2}[\cos (nx/L) + i\sin (nx/L)]\sin (y/L)\}, $$ (32) and set
$$\,\,{g_1} = {\tilde k^2} + {\tilde l^2}  {\beta /{\bar u}}, \quad$$ (33) $${{{g_2} = \varepsilon ({{\tilde k}^2} + {{\tilde l}^2})}/{(\bar u\tilde k)}}, $$ (34) then Eq. (31) becomes the following:
$$\begin{split} &\Psi (x, y) = \operatorname{Re} \{ \frac{{2{f_0}{h_2}}}{{{g_1}  i{g_2}}}[\cos ({{nx}/L}) + i\sin ({{nx}/L})]\sin ({y/L})\} \\ &= \operatorname{Re} \{ \frac{{2{f_0}{h_2}({g_1} + i{g_2})}}{{({g_1}  i{g_2})({g_1} + i{g_2})}}[\cos ({{nx}/L}) + i\sin ({{nx}/L})]\sin ({y/L})\}\\ & = \frac{{2{f_0}{h_2}}}{{g_1^2 + g_2^2}}[{g_1}\cos ({{nx}/L})  {g_2}\sin ({{nx}/L})]\sin ({y/L}). \end{split}$$ (35) Examples of calculated values of g_{1} and g_{2} are shown in Table 5. The absolute values of g_{1} are always much greater than that of g_{2}. Thus, the wave phase of the streamfunction
$\Psi $ relative to the topography mainly depends on the sign of g_{1}. If g_{1} is a positive (negative) value, the streamfunction$\Psi $ should be nearly in (out of) phase with the topography, in other words, ridges (troughs) should be over the mountains. Obviously, the wave phase of the lower layer streamfunction of these equilibrium states is exactly consistent with the wave phase predicted by this rough theory.In fact, the wave phase of equilibrium states depends on the direction of zonal wind and horizontal scale of the topography. Due to the conservation of potential vorticity, the absolute vorticity is decreased over the mountains. For westerly flow (
$\bar u > 0$ ), in the case${g_1} > 0$ , i.e.,${\tilde k^2} + {\tilde l^2} > {\beta /{\bar u}}$ [also called “long waves” case (Smith, 1979)], the decrease in absolute vorticity is primarily caused by the generation of negative relative vorticity, then ridges are generated over the mountains; by contrast, in the case${g_1} < 0$ , i.e.,${\tilde k^2} + {\tilde l^2} < {\beta /{\bar u}}$ [also called “ultralong waves” case (Smith, 1979)], the decrease in absolute vorticity is primarily caused by the decrease in planetary vorticity, which is associated with the southward movement of air parcels, and then troughs are generated over the mountains. However, for easterly flow ($\bar u < 0$ ), there is always${\tilde k^2} + {\tilde l^2} > {\beta /{\bar u}}$ , but the decrease in absolute vorticity arises both from the development of negative relative vorticity and from the decrease in planetary vorticity due to the southward motion (Holton and Hakim, 2012); then converse ridges are generated over the mountains. In our model, the occurrence of “long wave” case associated with ridgetype equilibria is purely due to the lower layer easterlies (Table 5 and Figs. 11b, 12b).In a word, the ridgetype (troughtype) equilibrium states have lower layer ridges (trough) over the mountains and have lower layer easterlies (westerlies). The wave phases of equilibrium states relative to the topography depends on the direction of lower layer zonal wind and horizontal scale of the topography. Further discussion is presented in Section 7.

Next, we investigate the high and lowindex equilibria and wave amplitude. To further examine the differences between these two types of equilibrium, we define some dimensional physical variables, as shown in the Table 4. All of the variables are defined over the domain (
$0 \leqslant x \leqslant 2\pi L$ ,$0 \leqslant y \leqslant \pi L$ ) of the channel. The amplitude of wave component of the upper layer geopotential height field is defined as$${\rm AH} = \frac{{{L^2}f_0^2}}{{{g_0}}}\sqrt {{{(\psi _2^1)}^2} + {{(\psi _3^1)}^2}} .$$ (36) The amplitude of wave components of the atmospheric and the land temperature fields are defined as
$$\!\!\!\!\!\!\!{\rm A}{{\rm T}_{\rm{a}}} = \frac{{2{L^2}f_0^2}}{R}\sqrt {{{({\theta _2})}^2} + {{({\theta _3})}^2}}, $$ (37) $${\rm A}{{\rm T}_{\rm{g}}} = \frac{{{L^2}f_0^2}}{R}\sqrt {{{({T_{{\rm g}, 2}})}^2} + {{({T_{{\rm g}, 3}})}^2}}, $$ (38) respectively. These two variables represent the zonal asymmetry of atmospheric and land temperature fields.
As expected, the wave amplitude AH of lowindex equilibrium state is always greater than that of highindex equilibrium state at the same value of C_{g} (Figs. 11d, 12d). This phenomenon also occurs in wave amplitude of the atmospheric temperature field AT_{a} (Figs. 11f, 12f). However, the meridional atmospheric temperature gradient ΔT_{a} of lowindex equilibrium state is always smaller than that of highindex equilibrium state at the same value of C_{g} (Figs. 11e, 12e).
These two types of equilibrium have no robust differences in the mean upper layer zonal wind speed: the Mean_U^{1} of the lowindex equilibrium state is smaller than that of the highindex equilibrium state at C_{g} = 54 W m^{–2 } for
$m = 3.7$ (Fig. 11a). By contrast, the former is greater than the latter at C_{g} = 22, 24 W m^{–2 } for$m = 6$ (Fig. 12a). Besides, the differences in value of Mean_U^{1} between high and lowindex equilibrium states are no more than 0.5 m s^{–1}. In addition, these two types of equilibrium also show no marked differences in nondimensional zonal component$\psi _1^1$ (see the overlap of the red circles and blue asterisks in Figs. 3a, d), which implies that the high and lowindex equilibria have no marked differences in upper layer zonal wind speed. Focusing on the middlelevel zonal wind speed (Figs. 11c, 12c), the lowindex equilibria always has a greater Mean_U^{2} than the highindex equilibria at the same value of C_{g}; however, their differences are also no more than 1.0 m s^{–1}.It is notable that our results regarding the differences between the high and lowindex equilibria in zonal wind differ from previous studies based on the barotropic models, in which the zonal component
${\psi _A}$ (i.e., zonal wind) of highindex equilibria was much greater than that of lowindex equilibria (see Fig. 1 in CD; Figs. 13a, 14a in Huang et al., 2017a). One may argue that the upper layer zonal component$\psi _1^1$ of the magenta branch equilibrium states, which are characterized by small wave amplitude (Figs. 3e, f), are greater than that of the lowindex equilibrium states of troughtype (Fig. 3d) in the baroclinic model, which can be an analogue of the results based on the barotropic model. However, the magenta branch equilibrium states are always unstable in the baroclinic model in this paper as well as in CS and RP. It should be noted that the results based on barotropic models disagree with the observations, e.g., some studies have shown that the probability density distribution of the zonal wind is unimodal (Benzi et al., 1986; Sutera, 1986). They are also inconsistent with the results of numerical experiments based on the general circulation model (Lindzen, 1986).Figure 13. Phase diagrams of dimensional variables for m = 3.7 (left panels) and m = 6 (right panels), respectively. Each ordinate shows the variable AH, and the abscissa gives (a, e) ΔT_{a}, (b, f) ΔT_{g}, (c, g) AT_{a}, and (d, h) AT_{g}. The meaning of colors and symbols are same as that in Fig. 3.
In fact, as emphasized in CS, the wavelike equilibrium is maintained not by the conversion of mean flow kinetic energy, but by the mean flow potential energy in the baroclinic atmosphere. Therefore, in our baroclinic model, as the lowindex equilibria has larger wave amplitudes (Figs. 11d, 12d), there is indeed a reduction in meridional atmospheric temperature gradient of lowindex equilibria (Figs. 11e, 12e) due to the consumption of mean flow potential energy. The low and highindex equilibria certainly have no marked differences in zonal wind speed in our baroclinic model (Figs. 11a, 12a). However, in the barotropic model, the wavelike equilibria could only obtain energy from the mean flow kinetic energy. Therefore, in the barotropic model, the lowindex equilibria with large wave amplitude undoubtedly has a lower zonal wind speed than the highindex equilibria with small wave amplitude.
The relationships between wave amplitude AH and meridional temperature gradient ΔT_{a} and ΔT_{g} are directly shown in Figs. 13a, b, e, f. As the wave amplitude of troughtype equilibria rapidly increases, the meridional atmospheric temperature gradient ΔT_{a} remarkably decreases for
$m = 3.7$ (Fig. 13a, blue branch) and slowly increases for$m = 6$ (Fig. 13e, blue branch). This is because more and more potential energy is consumed to maintain the troughtype equilibria with rapidly increasing wave amplitude. The wave amplitude of ridgetype equilibria increases much slowly, so the meridional atmospheric temperature gradient ΔT_{a} increases rapidly for both$m = 3.7$ and$m = 6$ (Figs. 13a, e, red branch). Of course, the wavelike equilibria hardly draws energy from the land directly, so there is no marked reduction of meridional land temperature gradient ΔT_{g} for the equilibria with large wave amplitude (Figs. 13b, f).In addition, regardless of ridge or troughtype as well as high or lowindex equilibria, the wave amplitudes of both atmospheric and land temperature fields (AT_{a} and AT_{g}) are highly positively correlated with wave amplitude AH (Figs. 13c, d, g, h). In fact, the atmospheric and land temperature fields of equilibrium states are always nearly in phase with each upper layer streamfunction field (Figs. 1, 2). If the atmospheric temperature field showed a lag or lead to the streamfunction field in phase, the meridional perturbations of streamfunction field would continue to grow or decay due to the temperature advection. Thus, there would be no stationary waves, i.e., equilibrium states. As we have only obtained equilibrium solutions from Eqs. (19)–(27), the atmospheric temperature field is surely in phase with the streamfunction field. The formation of the zonal asymmetric structure of the land temperature field should be attributed to the interactions between the land and atmospheric temperature fields through radiative and heat exchange. Therefore, the changes in wave amplitude of both atmospheric and land temperature fields are highly consistent with that of the upper layer streamfunction field (Figs. 11d, f; 12d, f; 13c, d g, h). This result also suggests that the wavelike equilibrium is maintained by the conversion of the mean flow potential energy.
The results in this section show that the lowindex (highindex) equilibrium states have a larger (smaller) wave amplitude and smaller (larger) meridional atmospheric temperature gradient; however, the two types equilibrium states have no marked differences in zonal wind speed. These results are attributed to the wavelike equilibrium that is maintained by the conversion of the mean flow potential energy in the baroclinic atmosphere.

A. Linearization of the quartic terms in the radiative fluxes
We assume that
$${T_{\rm{a}}}(x, y, t) = {T_{{\rm a}, 0}}(t) + \delta {T_{\rm{a}}}(x, y, t),\tag{A1}$$ (A1) $${T_{\rm{g}}}(x, y, t) = {T_{{\rm g}, 0}}(t) + \delta {T_{\rm{g}}}(x, y, t),\tag{A2}$$ (A2) where
${T_{{\rm a}, 0}}$ and${T_{{\rm g}, 0}}$ are spatially uniform averaged temperatures,$\delta {T_{\rm{a}}}$ and$\delta {T_{\rm{g}}}$ are temperature anomalies of the atmosphere and the land, respectively.We assume that the shortwave solar radiation absorbed by the atmosphere and the land are just the function of latitude and time, i.e.,
${R_{\rm{a}}}(y, t)$ and${R_{\rm{g}}}(y, t)$ , and we set$${R_{\rm{a}}}(y, t) = {R_{{\rm a}, 0}}(t) + \delta {R_{\rm{a}}}(y, t),\tag{A3}$$ (A3) $${R_{\rm{g}}}(y, t) = {R_{{\rm g}, 0}}(t) + \delta {R_{\rm{g}}}(y, t),\tag{A4}$$ (A4) where
${R_{{\rm a}, 0}}$ and${R_{{\rm g}, 0}}$ are time dependent spatially uniform shortwave solar radiation, and$\delta {R_{\rm{a}}}$ and$\delta {R_{\rm{g}}}$ are spatially varying counterparts.Neglecting the highorder terms in
$\delta T$ and separating the averaged temperatures and perturbations, i.e., the zerothorder terms in the expansion from the firstorder ones, the atmospheric temperature equation [Eq. (6)] becomes$${\gamma _{\rm{a}}}\frac{{\partial {T_{{\rm a}, 0}}}}{{\partial t}} =  \lambda ({T_{{\rm a}, 0}}  {T_{{\rm g}, 0}}) + {\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm g}, 0}^4  2{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm a}, 0}^4 + {R_{{\rm a}, 0}},\tag{A5}$$ (A5) $$\!\!\begin{split} & {\gamma _{\rm{a}}}(\frac{{\partial \delta {T_{\rm{a}}}}}{{\partial t}} + J(\psi, \delta {T_{\rm{a}}})  \sigma \omega \frac{p}{R}) =  \lambda (\delta {T_{\rm{a}}}  \delta {T_{\rm{g}}}) \\ &\quad\quad + 4{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm g}, 0}^3\delta {T_{\rm{g}}}  8{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm a}, 0}^3\delta {T_{\rm{a}}} + \delta {R_{\rm{a}}}, \end{split}\tag{A6}\qquad$$ (A6) and the land temperature equation [Eq. (7)] becomes
$${\gamma _{\rm{g}}}\frac{{\partial {T_{{\rm g}, 0}}}}{{\partial t}} =  \lambda ({T_{{\rm g}, 0}}  {T_{{\rm a}, 0}})  {\sigma _{\rm{B}}}T_{{\rm g}, 0}^4 + {\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm a}, 0}^4 + {R_{{\rm g}, 0}},\tag{A7}$$ (A7) $${\gamma _{\rm{g}}}\frac{{\partial \delta {T_{\rm{g}}}}}{{\partial t}} =  \lambda (\delta {T_{\rm{g}}}  \delta {T_{\rm{a}}})  4{\sigma _{\rm{B}}}T_{{\rm g}, 0}^3\delta {T_{\rm{g}}} + 4{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm a}, 0}^3\delta {T_{\rm{a}}} + \delta {R_{\rm{g}}}.\tag{A8}$$ (A8) Note that Eqs. (A5) and (A7) for the averaged temperatures are independent of the perturbations, and thus, stationary solutions can be obtained by solving
$$  \lambda ({T_{{\rm a}, 0}}  {T_{{\rm g}, 0}}) + {\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm g}, 0}^4  2{\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm a}, 0}^4 + {R_{{\rm a}, 0}} = 0,\tag{A9}$$ (A9) $$  \lambda ({T_{{\rm g}, 0}}  {T_{{\rm a}, 0}})  {\sigma _{\rm{B}}}T_{{\rm g}, 0}^4 + {\varepsilon _{\rm{a}}}{\sigma _{\rm{B}}}T_{{\rm a}, 0}^4 + {R_{{\rm g}, 0}} = 0.\tag{A10}\quad$$ (A10) According to the parameter values listed in Table 1, particularly,
$\lambda = 10$ W m^{–2} K^{–1}, we get${T_{{\rm a}, 0}} = 270.22$ K,${T_{{\rm g}, 0}} = 280.40$ K, and they are the default values in the main body. For$\lambda = 0$ W m^{–2} K^{–1} in the experiment Case 3, we get${T_{{\rm a}, 0}} = 264.16$ K,${T_{{\rm g}, 0}} = 295.71$ K, and they are only used in Section 4 of the main body. Since stationary solutions are obtained, Eqs. (A5) and (A7) need not be considered any more, and we just focus on Eqs. (A6), (A8), (4), (5), and (8).B. Equilibrium solutions and their stabilities
To obtain the general equilibrium solutions of Eqs. (19)–(27), we set all of the time derivatives to zero. We obtain
$$ k({\psi _1}  {\theta _1})  c\tilde h({\theta _3}  {\psi _3}) = 0,\tag{A11} \qquad\qquad\qquad\quad\quad $$ (A11) $$  c{n^2}({\psi _1}{\psi _3} + {\theta _1}{\theta _3}) + \beta n{\psi _3}  {B_1}({\psi _2}  {\theta _2}) = 0,\tag{A12}\qquad$$ (A12) $$c[{n^2}({\psi _1}{\psi _2} + {\theta _1}{\theta _2}) + \tilde h({\theta _1}  {\psi _1})]  \beta n{\psi _2}  {B_1}({\psi _3}  {\theta _3}) = 0, \tag{A13}$$ (A13) $$c[{\psi _2}{\theta _3}  {\psi _3}{\theta _2}  \sigma '\tilde h({\psi _3}  {\theta _3})]  {B_3}{\theta _1} + k\sigma '{\psi _1} + {D_1}{\theta _1} + {D_2} = 0,\tag{A14}$$ (A14) $$c({A_1}{\psi _3}{\theta _1}  {A_2}{\psi _1}{\theta _3}) + \beta n\sigma '{\theta _3}  {B_2}{\theta _2} + {B_1}\sigma '{\psi _2} + {D_1}{\theta _2} = 0,\tag{A15}\quad$$ (A15) $$\begin{split} & c[{A_2}{\psi _1}{\theta _2}  {A_1}{\psi _2}{\theta _1} + \sigma '\tilde h({\psi _1}  {\theta _1})]  \beta n\sigma '{\theta _2}\\ &\quad\quad  {B_2}{\theta _3} + {B_1}\sigma '{\psi _3} + {D_1}{\theta _3} = 0,\end{split}\tag{A16}\qquad\qquad\qquad\quad$$ (A16) and
$$\left. \begin{aligned} & {T_{{\rm g}, 1}} = \frac{{{d_4}}}{{{d_3}}}{\theta _1} + \frac{{{{C'}_{\rm{\!\!\!\!g}}}}}{{{d_3}}} \\ & {T_{{\rm g}, 2}} = \frac{{{d_4}}}{{{d_3}}}{\theta _2} \\ & {T_{{\rm g}, 3}} = \frac{{{d_4}}}{{{d_3}}}{\theta _3} \end{aligned} \right\}.\tag{A17}$$ (A17) In view of Eq. (A11), Eq. (A14) may be written as
$$c[{\psi _2}{\theta _3}  {\psi _3}{\theta _2}] + ({D_1}  2k'\sigma '){\theta _1} + {D_2} = 0,\tag{A18}$$ (A18) and Eqs. (A12), (A13), (A15), and (A16) become
$$\left[ {\begin{array}{*{20}{c}} {  {B_1}}&{{B_1}}&{  c{n^2}{\psi _1} + \beta n}&{  c{n^2}{\theta _1}} \\ {c{n^2}{\psi _1}  \beta n}&{c{n^2}{\theta _1}}&{  {B_1}}&{{B_1}} \\ {{B_1}\sigma '}&{{D_1}  {B_2}}&{c{A_1}{\theta _1}}&{  c{A_2}{\psi _1} + \beta n\sigma '} \\ {  c{A_1}{\theta _1}}&{c{A_2}{\psi _1}  \beta n\sigma '}&{{B_1}\sigma '}&{{D_1}  {B_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\psi _2}} \\[7.5pt] {{\theta _2}} \\[7.5pt] {{\psi _3}} \\[7.5pt] {{\theta _3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0 \\ {c\tilde h({\psi _1}  {\theta _1})} \\ 0 \\ {  c\sigma '\tilde h({\psi _1}  {\theta _1})} \end{array}} \right].\tag{A19}$$ (A19) Equation (A19) constitutes a linear system for the variables
$({\psi _2},\;{\theta _2}, \;{\psi _3}, \;{\theta _3})$ if the zonal variables$({\psi _1}, \;{\theta _1})$ are specified. The general solution of Eq. (A19) is described in detail in Supplementary Section 1. If we set$\tilde h = 0$ , then the right side of Eq. (A19) becomes a zero matrix. In general, the value of the coefficient determinant on the left side is not equal to zero. Thus, in this case, the solutions of wave components$({\psi _2}, \;{\theta _2}, \;{\psi _3}, \;{\theta _3})$ are all equal to zero. Therefore, the wavelike equilibria (wave components are not zero, i.e., stationary wave solution) cannot exist without the topography in this system.The stability of the equilibrium solution obtained from Eqs. (19)–(27) is determined from the characteristic values of the linear perturbation equations coefficient matrix, which is a nine homogeneous linear equations governing
$({\psi _1}, \;{\psi _2}, \;{\psi _3}, \;{\theta _1}, \;{\theta _2}, \;{\theta _3}, \;{T_{g, 1}}, \;{T_{g, 2}}, \;{T_{g, 3}})$ (see Supplementary Section 1). Set all perturbation quantities be proportional to${e^{ + \alpha t}}$ , we obtain a nineorder equation in the variable$\alpha $ . If the maximum real part of$\alpha $ is greater than zero, the equilibrium is unstable, otherwise it is stable.