Phase Shift of Planetary Waves and Wave–Jet Resonance on Tidally Locked Planets

A 2D linear shallow-water model in equatorial β-plane with Matsuno-Gill-type forcing (Matsuno 1966; Gill 1980; Showman & Polvani 2011; Heng & Workman 2014; Penn & Vallis 2017) is used in this study. The model is a 1.5-layer shallow-water system, including an active upper layer of constant density that represents the free troposphere and underlying a quiescent layer representing the lower troposphere. Atmospheric flows are forced by steady heating and cooling (corresponding to mass source and sink in the model) and damped by radiative relaxation and linear friction. The heating and cooling are assumed small enough for linear theory to apply. For easy to obtain analytic solutions, a uniform mean flow (U), representing an atmospheric superrotation, is specified in the model, same as that employed in Phlips & Gill (1987), Arnold et al. (2012) and Tsai et al. (2014). To enable analytic solutions, nonlinear advection terms are not considered, except the zonal advection by the specified mean flow. These lead to a linear system for the flow in the upper layer:

$$\begin{eqnarray}&&\displaystyle \frac{\partial u}{\partial t}+U\displaystyle \frac{\partial u}{\partial x}-\beta {yv}+g\displaystyle \frac{\partial \eta }{\partial x}=-\displaystyle \frac{u}{{\tau }_{\mathrm{drag}}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial v}{\partial t}+U\displaystyle \frac{\partial v}{\partial x}+\beta {yu}+g\displaystyle \frac{\partial \eta }{\partial y}=-\displaystyle \frac{v}{{\tau }_{\mathrm{drag}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \eta }{\partial t}+U\displaystyle \frac{\partial \eta }{\partial x}+H\left(\displaystyle \frac{\partial u}{\partial x}+\displaystyle \frac{\partial v}{\partial y}\right)=S(x,y)-\displaystyle \frac{\eta }{{\tau }_{\mathrm{rad}}},\end{eqnarray} \tag{ 3 }$$

where t is the time, x is the eastward distance, y is the northward distance, u is the zonal velocity, v is the meridional velocity, β = df/dy is the gradient of the Coriolis parameter with the northward distance, η is the height anomaly from the reference fluid height in the absence of forcing (H), and g is the reduced gravity (i.e., gravity times the fractional density difference between the layers Δρ/ρ). Momentum dissipation is represented by Rayleigh friction with a drag timescale of τ_drag. Heating and cooling in the atmosphere are represented by mass source and sink, respectively, writing as S = (h_eq  −  H)/τ_rad, where τ_rad is the radiative relaxation timescale (i.e., the timescale for the atmosphere to reach local radiative equilibrium) and h_eq is the 2D radiative-equilibrium height. These linear equations imply that the height anomaly is much smaller than the average fluid thickness (η  ≪  H). Momentum exchange between the upper layer and the lower layer (Showman & Polvani 2010) is neglected in the model. The momentum exchange term is important in generating the superrotation in a nonlinear shallow-water system (Showman & Polvani 2010), but it is unimportant in a linear shallow-water system when investigating the behavior of waves rather than the superrotation. In math, the momentum exchange term is quadratic and does not appear in the linear shallow-water equations (Showman & Polvani 2011). Similar to Showman & Polvani (2011), we have ignored the term here in order to obtain analytical solution, and moreover the superrotation wind is specified in our model.

We nondimensionalize the equations with a velocity scale of $c=\sqrt{{gH}}$ corresponding to the speed of gravity wave (such as the nondispersive Kelvin wave), a length scale of ${L}_{R}={\left(c/\beta \right)}^{1/2}$ corresponding to the equatorial Rossby deformation radius, and a timescale of $T={\left(c\beta \right)}^{-1/2}$ corresponding to the time for the gravity wave to propagate throughout the distance of L_R. The height scale is equal to H, the dissipation and radiative timescales are nondimensionalized with T, and S with H/T. These characteristic scales are the same as those used in Matsuno (1966) and Showman & Polvani (2011). For steady solutions under forcing and damping, the equations become

$$\begin{eqnarray}&&{U}^{* }\displaystyle \frac{\partial {u}^{* }}{\partial {x}^{* }}+\displaystyle \frac{\partial {\eta }^{* }}{\partial {x}^{* }}-{y}^{* }{v}^{* }=-\displaystyle \frac{{u}^{* }}{{\tau }_{\mathrm{drag}}^{* }},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{U}^{* }\displaystyle \frac{\partial {v}^{* }}{\partial {x}^{* }}+\displaystyle \frac{\partial {\eta }^{* }}{\partial {y}^{* }}+{y}^{* }{u}^{* }=-\displaystyle \frac{{v}^{* }}{{\tau }_{\mathrm{drag}}^{* }},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{U}^{* }\displaystyle \frac{\partial {\eta }^{* }}{\partial {x}^{* }}+\left(\displaystyle \frac{\partial {u}^{* }}{\partial {x}^{* }}+\displaystyle \frac{\partial {v}^{* }}{\partial {y}^{* }}\right)={S}^{* }({x}^{* },{y}^{* })-\displaystyle \frac{{\eta }^{* }}{{\tau }_{\mathrm{rad}}^{* }},\end{eqnarray} \tag{ 6 }$$

within which all variables with stars are nondimensional. For brevity, we omit stars in the following text. The model now contains three main important parameters: U, τ_drag, and τ_rad. Our goal in Section 3 below is to determine the dependence of the phase shift and the resonance on these three parameters.

For 1:1 tidally locked planets, the forcing is a sinusoidal function in longitude and a parabolic cylinder function in latitude, in order to approximately mimic the mass source on the dayside and the mass sink on the nightside, same as that employed in previous studies (such as Showman & Polvani 2010, 2011), written as

$$\begin{eqnarray}&&S(x,y)=\displaystyle \frac{{\rm{\Delta }}{h}_{{eq}}}{{\tau }_{\mathrm{rad}}}\displaystyle \sum _{n}{\psi }_{n}(y){e}^{{ikx}},\end{eqnarray} \tag{ 7 }$$

where ψ_n(y) are parabolic cylinder functions, subscript n is the order of the functions (Abramowitz & Stegun 1965), and Δh_eq represents the dimensionless strength of the radiative-equilibrium force. The substellar point is at kx = 0 and the antistellar point is at kx = π, thus positive e^ikx corresponds to mass source on the dayside and negative corresponds to mass sink on the nightside, as shown in Figure 2(b). Due to the hemisphere-scale mass source/sink on tidally locked planets, we consider ψ₀(y) only in the following solutions. Note that a flat nightside is more realistic according to the spatial pattern of the stellar radiation (Perez-Becker & Showman 2013; Showman et al. 2013; Komacek & Showman 2016; Zhang & Showman 2017), but it is not used in this work; this is due to that obtained analytic solution is not easy under a flat nightside when the radiative timescale and the dissipation timescale are unequal. Furthermore, the waves do not exhibit essential differences when comparing a flat nightside to a sinusoidal function on the nightside, although the detailed spatial pattern and amplitude of the waves do vary (Perez-Becker & Showman 2013).

Zoom In Zoom Out Reset image size

The model used here is similar to that in Showman & Polvani (2011) but a mean flow is included, so that our model is able to consider the effect of the mean flow on the waves. Tsai et al. (2014) also considered the effect of a uniform mean flow on the waves in both 2D and 3D linear models. But, the model used here is relatively easier to uncover the underlying mechanisms. Moreover, in the work of Tsai et al. (2014), the radiative timescale is assumed to be equal to the dissipation timescale (i.e., τ_rad = τ_drag) due to the constraint in their method of obtaining analytic solutions, while, in our model, this assumption is unnecessary, so that we can discuss the separate effects of varying τ_rad and varying τ_drag. Hammond & Pierrehumbert (2018) introduce a nonuniform jet U(y) in their 2D shallow-water model, which is more realistic and can better match their results up with 3D atmospheric circulation simulations in the horizontal height field, but it was inconvenient for obtaining the results analytically; a pseudo-spectral method is required for a nonuniform jet (Hammond & Pierrehumbert 2018).

To solve Equations (4)–(6), we follow the method used in Showman & Polvani (2011) but including the mean jet of U. Due to the longitudinally sinusoidal form of the forcing (see Figure 2(b)), all the variables would have the same form in longitude:

$$\begin{eqnarray}&&\left[u,v,\eta ,S\right]=\left[u(y),v(y),\eta (y),S(y)\right]{e}^{{ikx}}.\end{eqnarray} \tag{ 8 }$$

The goal is to solve for the unknowns of u, v, and η, under given S(y). In order to solve the equations, it is convenient to introduce several variables: $\alpha \equiv {\tau }_{\mathrm{drag}}^{-1}$, $\gamma \equiv {\tau }_{\mathrm{rad}}^{-1}$, $\tilde{\alpha }\equiv \alpha \,+{ikU}$, $\tilde{\gamma }\equiv \gamma +{ikU}$, $q(y)=\sqrt{\tilde{\gamma }}\eta (y)+\sqrt{\tilde{\alpha }}u(y)$, and $r(y)\,=\sqrt{\tilde{\gamma }}\eta (y)-\sqrt{\tilde{\alpha }}u(y)$, where q(y) and r(y) are complex variables due to the introduction of the zonal-mean jet U. So, we can convert Equations (4)–(6) to three equivalent equations for q(y), r(y), and v(y):

$$\begin{eqnarray}&&{ikq}(y)+\sqrt{\tilde{\alpha }}\displaystyle \frac{{dv}(y)}{{dy}}=\sqrt{\tilde{\gamma }}{yv}(y)+\sqrt{\tilde{\alpha }}S(y)-\sqrt{\tilde{\alpha }\tilde{\gamma }}q(y),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{ikr}(y)-\sqrt{\tilde{\alpha }}\displaystyle \frac{{dv}(y)}{{dy}}=\sqrt{\tilde{\gamma }}{yv}(y)-\sqrt{\tilde{\alpha }}S(y)+\sqrt{\tilde{\alpha }\tilde{\gamma }}r(y),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{yq}(y)-{yr}(y)+\sqrt{\displaystyle \frac{\tilde{\alpha }}{\tilde{\gamma }}}\displaystyle \frac{{dq}(y)}{{dy}}+\sqrt{\displaystyle \frac{\tilde{\alpha }}{\tilde{\gamma }}}\displaystyle \frac{{dr}(y)}{{dy}}=-2{\tilde{\alpha }}^{3/2}v(y).\end{eqnarray} \tag{ 11 }$$

The eigenfunctions of q(y), r(y), v(y), and S(y) are parabolic cylinder functions:

$$\begin{eqnarray}&&\left[q(y),r(y),v(y),S(y)\right]=\displaystyle \sum _{n}\left[{\hat{q}}_{n},{\hat{r}}_{n},{\hat{v}}_{n},{\hat{S}}_{n}\right]{\psi }_{n}(y),\end{eqnarray} \tag{ 12 }$$

where ${\hat{S}}_{0}=\gamma {\rm{\Delta }}{h}_{{eq}}$. The parabolic cylinder functions (ψ_n(y)) and their recursion relations are

$$\begin{eqnarray}&&{\psi }_{n}(y)=\exp \left(-\displaystyle \frac{{y}^{2}}{2{{ \mathcal P }}^{2}}\right){H}_{n}\left(\displaystyle \frac{y}{{ \mathcal P }}\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\psi }_{n}}{{dy}}=\displaystyle \frac{2n{\psi }_{n-1}}{{ \mathcal P }}-\displaystyle \frac{y{\psi }_{n}}{{{ \mathcal P }}^{2}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\psi }_{n}}{{dy}}=-\displaystyle \frac{{\psi }_{n+1}}{{ \mathcal P }}+\displaystyle \frac{y{\psi }_{n}}{{{ \mathcal P }}^{2}},\end{eqnarray} \tag{ 15 }$$

where H_n(y) is the Hermitian polynomial and ${ \mathcal P }\equiv {\left(\tilde{\alpha }/\tilde{\gamma }\right)}^{1/4}$ (Abramowitz & Stegun 1965). Substituting Equation (12) into Equations (9)–(11) and using the recursion relations, we obtain three series equations:

$$\begin{eqnarray}&&{ik}{\hat{q}}_{n}+\sqrt{\tilde{\alpha }\tilde{\gamma }}{\hat{q}}_{n}-{\left(\tilde{\alpha }\tilde{\gamma }\right)}^{1/4}{\hat{v}}_{n-1}=\sqrt{\tilde{\alpha }}{\hat{S}}_{n},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{ik}{\hat{r}}_{n}-\sqrt{\tilde{\alpha }\tilde{\gamma }}{\hat{r}}_{n}-2(n+1){\left(\tilde{\alpha }\tilde{\gamma }\right)}^{1/4}{\hat{v}}_{n+1}=-\sqrt{\tilde{\alpha }}{\hat{S}}_{n},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&2(n+1){\hat{q}}_{n+1}-{\hat{r}}_{n-1}=-2{\tilde{\alpha }}^{3/2}{{ \mathcal P }}^{-1}{\hat{v}}_{n},\end{eqnarray} \tag{ 18 }$$

where the subscript n ≥ 0, and the terms will be zero if the subscript n < 0. The given forcing is ${\hat{S}}_{0}$, and the order n is no more than 2, as mentioned in Gill (1980) and Showman & Polvani (2011). For n = 0, Equation (16) implies

$$\begin{eqnarray}&&{\hat{q}}_{0}=\displaystyle \frac{\sqrt{\tilde{\alpha }}\left(\sqrt{\tilde{\alpha }\tilde{\gamma }}-{ik}\right)}{\tilde{\alpha }\tilde{\gamma }+{k}^{2}}{\hat{S}}_{0}.\end{eqnarray} \tag{ 19 }$$

Substituting n = 2, n = 0, and n = 1 into Equations (16)–(18), we obtain:

$$\begin{eqnarray}&&{ik}{\hat{q}}_{2}+\sqrt{\tilde{\alpha }\tilde{\gamma }}{\hat{q}}_{2}-{\left(\tilde{\alpha }\tilde{\gamma }\right)}^{1/4}{\hat{v}}_{1}=0,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{ik}{\hat{r}}_{0}-\sqrt{\tilde{\alpha }\tilde{\gamma }}{\hat{r}}_{0}-2{\left(\tilde{\alpha }\tilde{\gamma }\right)}^{1/4}{\hat{v}}_{1}=-\sqrt{\tilde{\alpha }}{\hat{S}}_{0},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&4{\hat{q}}_{2}-{\hat{r}}_{0}=-2{\tilde{\alpha }}^{3/2}{{ \mathcal P }}^{-1}{\hat{v}}_{1}.\end{eqnarray} \tag{ 22 }$$

Eliminating ${\hat{r}}_{0}$ and ${\hat{v}}_{1}$, we obtain

$$\begin{eqnarray}&&{\hat{q}}_{2}=\displaystyle \frac{{\tilde{\alpha }}^{3/2}{k}^{2}+3\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{5/2}\tilde{\gamma }+{ik}{\tilde{\alpha }}^{1/2}}{2\left[{\left(\tilde{\alpha }{k}^{2}+3{\tilde{\alpha }}^{1/2}{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{2}\tilde{\gamma }\right)}^{2}+{k}^{2}\right]}{\hat{S}}_{0}.\end{eqnarray} \tag{ 23 }$$

From Equation (20) and the analytic solution of ${\hat{q}}_{2}$, we get an analytic solution for ${\hat{v}}_{1}$:

$$\begin{eqnarray}&&{\hat{v}}_{1}=\displaystyle \frac{\sqrt{\tilde{\alpha }\tilde{\gamma }}+{ik}}{{\left(\tilde{\alpha }\tilde{\gamma }\right)}^{1/4}}{\hat{q}}_{2}.\end{eqnarray} \tag{ 24 }$$

Combining Equations (22) and (24), we obtain ${\hat{r}}_{0}\,=\left(4+2{\tilde{\alpha }}^{3/2}\sqrt{\tilde{\gamma }}+2{ik}\tilde{\alpha }\right){\hat{q}}_{2}.$ Note that ${\hat{r}}_{2}$ is equal to zero. The final solution is:

$$\begin{eqnarray}\begin{array}{rcl}\hat{\eta } & = & \displaystyle \frac{{\hat{q}}_{0}{\psi }_{0}+{\hat{q}}_{2}{\psi }_{2}+{\hat{r}}_{0}{\psi }_{0}}{2\sqrt{\tilde{\gamma }}}{e}^{{ikx}}\\ & = & {\hat{S}}_{0}\left[\sqrt{\displaystyle \frac{\tilde{\alpha }}{\tilde{\gamma }}}\displaystyle \frac{\sqrt{\tilde{\alpha }\tilde{\gamma }}-{ik}}{2(\tilde{\alpha }\tilde{\gamma }+{k}^{2})}{e}^{{ikx}}\right]\exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right)\\ & & +{\hat{S}}_{0}\left[\displaystyle \frac{{\tilde{\alpha }}^{3/2}{k}^{2}+3\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{5/2}\tilde{\gamma }+{ik}{\tilde{\alpha }}^{1/2}}{4\sqrt{\tilde{\gamma }}\left[{\left(\tilde{\alpha }{k}^{2}+3{\tilde{\alpha }}^{1/2}{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{2}\tilde{\gamma }\right)}^{2}+{k}^{2}\right]}{e}^{{ikx}}\right]\\ & & \times \left(4\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}{y}^{2}-2\right)\exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right)\\ & & +{\hat{S}}_{0}\left[\left(2{\tilde{\gamma }}^{-1/2}+{\tilde{\alpha }}^{3/2}+{ik}\tilde{\alpha }{\tilde{\gamma }}^{-1/2}\right)\right.\\ & & \times \left.\displaystyle \frac{{\tilde{\alpha }}^{3/2}{k}^{2}+3\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{5/2}\tilde{\gamma }+{ik}{\tilde{\alpha }}^{1/2}}{2\left[{\left(\tilde{\alpha }{k}^{2}+3{\tilde{\alpha }}^{1/2}{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{2}\tilde{\gamma }\right)}^{2}+{k}^{2}\right]}{e}^{{ikx}}\right]\\ & & \times \exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right),\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}\begin{array}{rcl}\hat{u} & = & \displaystyle \frac{{\hat{q}}_{0}{\psi }_{0}+{\hat{q}}_{2}{\psi }_{2}-{\hat{r}}_{0}{\psi }_{0}}{2\sqrt{\tilde{\alpha }}}{e}^{{ikx}}\\ & = & {\hat{S}}_{0}\left[\displaystyle \frac{\sqrt{\tilde{\alpha }\tilde{\gamma }}-{ik}}{2(\tilde{\alpha }\tilde{\gamma }+{k}^{2})}{e}^{{ikx}}\right]\exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right)\\ & & +{\hat{S}}_{0}\left[\displaystyle \frac{{\tilde{\alpha }}^{3/2}{k}^{2}+3\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{5/2}\tilde{\gamma }+{ik}{\tilde{\alpha }}^{1/2}}{4\sqrt{\tilde{\alpha }}\left[{\left(\tilde{\alpha }{k}^{2}+3{\tilde{\alpha }}^{1/2}{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{2}\tilde{\gamma }\right)}^{2}+{k}^{2}\right]}{e}^{{ikx}}\right]\\ & & \times \left(4\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}{y}^{2}-2\right)\exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right)\\ & & -{\hat{S}}_{0}\left[\left(2{\tilde{\alpha }}^{-1/2}+\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{ik}{\tilde{\alpha }}^{1/2}\right)\right.\\ & & \times \left.\displaystyle \frac{{\tilde{\alpha }}^{3/2}{k}^{2}+3\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{5/2}\tilde{\gamma }+{ik}{\tilde{\alpha }}^{1/2}}{2\left[{\left(\tilde{\alpha }{k}^{2}+3{\tilde{\alpha }}^{1/2}{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{2}\tilde{\gamma }\right)}^{2}+{k}^{2}\right]}{e}^{{ikx}}\right]\\ & & \times \exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right),\end{array}\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\hat{v} & = & {\hat{v}}_{1}{\psi }_{1}{e}^{{ikx}}\\ & = & {\hat{S}}_{0}\left[\displaystyle \frac{(\sqrt{\tilde{\alpha }\tilde{\gamma }}+{ik})({\tilde{\alpha }}^{3/2}{k}^{2}+3\tilde{\alpha }{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{5/2}\tilde{\gamma }+{ik}{\tilde{\alpha }}^{1/2})}{{\left(\tilde{\alpha }\tilde{\gamma }\right)}^{1/4}\left[{\left(\tilde{\alpha }{k}^{2}+3{\tilde{\alpha }}^{1/2}{\tilde{\gamma }}^{1/2}+{\tilde{\alpha }}^{2}\tilde{\gamma }\right)}^{2}+{k}^{2}\right]}{e}^{{ikx}}\right]\\ & & \times {\left(\tilde{\gamma }/\tilde{\alpha }\right)}^{1/4}y\exp \left(-\sqrt{\displaystyle \frac{\tilde{\gamma }}{\tilde{\alpha }}}\displaystyle \frac{{y}^{2}}{2}\right),\end{array}\,\end{eqnarray} \tag{ 27 }$$

where $\hat{\eta }$, $\hat{u}$, and $\hat{v}$ are complex quantities and η, u, and v are the real parts of them.

The solution can be decomposed to Kelvin and Rossby components (Gill 1980; Vallis 2006). For the Kelvin wave, its meridional velocity (v) is equal to zero, and its u and η satisfy ${\eta }_{K}=\tfrac{{\hat{q}}_{0}{\psi }_{0}}{2\sqrt{\tilde{\gamma }}}{e}^{{ikx}}$ and ${u}_{K}=\tfrac{{\hat{q}}_{0}{\psi }_{0}}{2\sqrt{\tilde{\alpha }}}{e}^{{ikx}}$, respectively. The Rossby waves are ${\eta }_{R}=\tfrac{{\hat{q}}_{2}{\psi }_{2}+{\hat{r}}_{0}{\psi }_{0}}{2\sqrt{\tilde{\gamma }}}{e}^{{ikx}}$, ${u}_{R}=\tfrac{{\hat{q}}_{2}{\psi }_{2}-{\hat{r}}_{0}{\psi }_{0}}{2\sqrt{\tilde{\alpha }}}{e}^{{ikx}}$, and ${v}_{R}={\hat{v}}_{1}{\psi }_{1}{e}^{{ikx}}$. Below, Figures 3 and 4 show the typical solutions under different values of τ_rad, τ_drag, and U.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The atmospheric circulations under zero and non-zero mean flows are shown in Figure 3. When U = 0, the spatial pattern of the circulation results from the superposition of equatorial Kelvin wave that propagates to the east and off-equatorial Rossby waves that propagate to the west (see the leftmost panels). When the drag or cooling effect is small, an outstanding feature of the circulation pattern is the high-pressure center (i.e., the local maximum of the height anomaly (η); zonal-mean values have been subtracted) and low-pressure center (i.e., the local minimum of η) located off the equator. In the vicinity of the equator, the flow is directing from the dayside to the nightside, caused by the surface inclination associated with the mass sources and sinks. In higher latitudes, due to the effect of the Coriolis force, the flow field is approximately in geostrophical balance with the pressure (surface elevation) field, within which anticyclonic or cyclonic flow fields are established where high- or low-pressure cells are located (Matsuno 1966). The flow off the equator is connected to the flow on the equator through converging or diverging motions toward or from the equator at the end or beginning of each cell. When the drag and cooling effect is strong (top left corner of Figure 3), the spatial pattern is analogous to the source/sink. However, wave pattern still exhibits slight northwest-southeast tilting in the northern hemisphere and southwest-northeast tilting in the southern hemisphere although the zonal propagation of Rossby and Kelvin waves is nearly inhibited. This is due to the fact that the drag and Coriolis forces are comparable, and a three-way horizontal force balance between those two forces and pressure-gradient force leads to the tilting (Showman et al. 2013). Note that inertia gravity waves are not included in the solutions. This is due to the facts that we are looking for a steady solution under stationary forcing and that the frequencies of the Rossby and Kelvin waves are much smaller than those of inertia gravity waves (Matsuno 1966).

When U = 0, the eddy high-pressure center (marked with red dot in each panel) is around the substellar point when the radiation timescale (τ_rad) or the drag timescale (τ_drag) is short, or on the west of the substellar point when τ_rad and τ_drag are intermediate or long. When τ_rad (or τ_drag) is short, the effect of radiation relaxation (or friction damping) is strong, so that the waves are unable to effectively propagate zonally, the height anomaly field is more analogous to the radiative-equilibrium height, and therefore the high-pressure center is close to the substellar point and the low pressure is close to the antistellar point. When τ_rad and τ_drag are long, the eddy height anomaly field is dominated by off-equatorial Rossby waves, which effectively propagate toward the west, so that the high-pressure centers (corresponding to wave crest) are on the west of the substellar point and the low-pressure centers (corresponding to wave trough) are on the east of the substellar point (the leftmost panels in Figure 3, see also Figure 3 in Showman & Polvani (2011)).

As U > 0, a clear response is that both the high and low-pressure centers shift toward the east, meaning an eastward shift in the phase of the waves (Figure 3). The eastward phase shift is mainly from the Rossby waves while the Kelvin wave component remains almost stationary (figure not shown), similar to that found in Figure 7 of Tsai et al. (2014). This is due to the fact that the phase speed of the Kelvin wave ($\sqrt{{gH}}$) is always faster than the mean flow, so it is less affected. The magnitude of the phase shift depends on the magnitude of U and on the values of τ_rad and τ_drag (Figures 5(a) and (b)). For a given U, the magnitude of the phase shift is smaller when τ_rad or τ_drag is shorter (i.e., the radiative relaxation or the drag friction is stronger), and vice versa. Mathematically speaking, Equations (4)–(6) shows that when the radiation relaxation term and the friction term are strong, the effect of the advection terms would be relatively smaller.

Zoom In Zoom Out Reset image size

For U < 0, both the high- and low-pressure centers shift toward the west (Figure 4), meaning a westward phase shift of the waves. The westward phase shift is mainly from the Kelvin wave while the Rossby wave component remains almost stationary (figure not shown). The magnitude of the westward phase shift also depends on the magnitude of U and on the values of τ_rad and τ_drag, and it is relatively small (large) when τ_rad or τ_drag is shorter (longer) (Figure 5(c) and (d)).

Figure 5 shows that, under a given τ_rad and τ_drag, the degree of the phase shift is a monotonic, nonlinear, increasing function of U and more importantly there are limits for the phase shift, +π for U > 0 and -π for U < 0. When ∣U∣ is increased, the transition from a small phase shift to the upper limits of ±π under large τ_rad and τ_drag is sharper than that under small τ_rad and τ_drag. The limits of ±π in the phase shift can be explained based on the energy balance of the whole system. For Equations (1)–(3), the corresponding global energy equation is given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\displaystyle \int \displaystyle \int ({K}_{E}+{P}_{E}){dxdy}=\displaystyle \int \displaystyle \int (\eta S){dxdy}\\ \quad -\displaystyle \int \displaystyle \int \left(\displaystyle \frac{2{K}_{E}}{{\tau }_{\mathrm{drag}}}+\displaystyle \frac{2{P}_{E}}{{\tau }_{\mathrm{rad}}}\right){dxdy},\end{array}\end{eqnarray} \tag{ 28 }$$

where K_E ≡ (u² + v²)/2 and P_E ≡ η²/2 represent nondimensional kinetic energy and potential energy, respectively. This equation is obtained through multiplying Equation (1) by u, multiplying Equation (2) by v, multiplying Equation (3) by η, and then spatially integrating the equations under periodic boundary condition in the x direction and zero normal velocity condition at the north and south side boundaries.

The first term on the right side of Equation (28) is the correlation between the height anomaly (η) and the forcing (S), which represents the exchange between the source energy and the kinetic and potential energy of the waves. The second term is the damping and friction effects in the system. Because the second term is negative everywhere, the first term must be positive at least in global mean and thereby a positive correlation between η and S is required. From Figure 6(a), it is clear to see that under a westward mean flow of such as $U=-2\sqrt{{gH}}$, the phase shift reaches the lower limit of −π; this is because a further westward shift of η will cause a zero or negative integration of ηS and a steady state cannot reach. Similarly, under an eastward mean flow of such as $U=+\sqrt{{gH}}$, the phase shift reaches the upper limit of +π, because a further eastward shift of η will cause a zero or negative integration of ηS (Figure 6(b)). These results mean that the center of the Kelvin wave cannot shift beyond the west terminator and the center of the Rossby waves cannot shift beyond the east terminator.

Zoom In Zoom Out Reset image size

The two limits in the phase shift can also be understood based on a solution at the equator. At y = 0, the shallow-water equations (Equations (4)–(6) and Equation (8)) reduce to

$$\begin{eqnarray}&&\tilde{\alpha }{u}^{0}+{ik}{\eta }^{0}=0,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\tilde{\gamma }{\eta }^{0}+{{iku}}^{0}=\gamma {S}_{0}-{v}_{y}^{0},\end{eqnarray} \tag{ 30 }$$

where the superscript 0 means quantities at y = 0, and the subscript y represents the partial derivative with respect to the latitude. The solution of the height field is

$$\begin{eqnarray}\begin{array}{rcl}{\eta }^{0} & = & \displaystyle \frac{{\alpha }^{2}\gamma +\alpha {k}^{2}+\gamma {k}^{2}{U}^{2}+{ikU}\left({k}^{2}-{\alpha }^{2}-{k}^{2}{U}^{2}\right)}{{\left(\alpha \gamma -{k}^{2}{U}^{2}+{k}^{2}\right)}^{2}+{k}^{2}{U}^{2}{\left(\alpha +\gamma \right)}^{2}}\\ & & \times \left[\gamma {S}_{0}-2{\left(\displaystyle \frac{\gamma +{ikU}}{\alpha +{ikU}}\right)}^{1/4}{\hat{v}}_{1}\right],\end{array}\end{eqnarray} \tag{ 31 }$$

where ${\hat{v}}_{1}$ is the same as that in Equation (24). When α and γ are small (i.e., τ_drag and τ_rad are large) and U is large enough (kU ≫ α, γ), ${\hat{v}}_{1}$ is proportion to ${\tilde{\alpha }}^{-2}$ and hence U⁻² (ref. Equation (27)), so it tends to be close to 0. Thereby, the value of η⁰ in Equation (31) is close to  − iγS₀/kU. This means that the absolute value of the maximum phase deviation from S₀ is π/2 and the upper limit of the phase shift is +π. The same principle applies to the lower limit.

Besides the phase shift of planetary waves, mean flow can also influence the amplitude of the waves. As shown in Figure 7, the amplitude of the waves (compared to the condition under a zero mean flow) exhibits a resonance behavior. For U > 0, the wave amplitude reaches a peak when the value of U approaches the westward phase speed of the Rossby wave, which is $\sqrt{{gH}}/3$. For U < 0, the wave amplitude also reaches a peak when the absolute value of U approaches the eastward phase speed of the Kelvin wave, which is $\sqrt{{gH}}$. When τ_rad or τ_drag is small, the resonance behavior is weak or absent due to the strong relaxing or damping effect, while when τ_rad and τ_drag are large, the resonance is more significant. This resonance phenomena has also been found in the previous studies of Arnold et al. (2012), Tsai et al. (2014), and Herbert et al. (2020) but only for identical τ_rad and τ_drag.

Zoom In Zoom Out Reset image size

The physical mechanism for the resonance can be interpreted based on the energy balance equation of Equation (28). When the mean flow speed is equal to the phase speed of the Kelvin wave (left panel in Figure 8) or that of the Rossby wave but in opposite sign (right panel in Figure 8), the waves are right trapped in the mass source and sink regions, so that the correlation between the forcing (S) and the wave response (η) reaches a peak. This also means that the energy conversation from the source to the kinetic and potential energy of the waves reaches a maximum, as shown in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The resonance behavior can also be explained based on a simplified analytic solution. When τ_rad is equal to τ_drag, the solution of η reduces to:

$$\begin{eqnarray}\begin{array}{rcl}\eta & = & \displaystyle \frac{{S}_{0}}{2\left[{\alpha }^{2}+{k}^{2}{\left(U+1\right)}^{2}\right]}(\alpha \cos ({kx})\\ & & +k(U+1)\sin ({kx}))\exp (-{y}^{2}/2),\end{array}\end{eqnarray} \tag{ 32 }$$

for the Kelvin component, and

$$\begin{eqnarray}\begin{array}{rcl}\eta & = & \displaystyle \frac{{S}_{0}}{6[{\alpha }^{2}+{k}^{2}{\left(U-1/3\right)}^{2}]}\left(\alpha \cos ({kx})\right.\\ & & +\left.k(U-1/3)\sin ({kx}\right)\left(2{y}^{2}+1\right)\exp (-{y}^{2}/2),\end{array}\end{eqnarray} \tag{ 33 }$$

for the Rossby component, similar to that found in Arnold et al. (2012). From these two equations, it is clear to see that the amplitude of the Kelvin wave reaches a maximum when the speed of the mean flow is approximately equal to $-\sqrt{{gH}}$, and the amplitude of the Rossby wave reaches a maximum when the speed of the mean flow is approximately equal to $\sqrt{{gH}}/3$, if all else is equal. The same explanation can be found in Equation (7) of Tsai et al. (2014): the wave responses are maximum when the Doppler-shifted frequency is equal to the frequency for free modes.

Note that, because of the interaction between the Rossby and Kelvin waves, the mean flow speed for the occurrence of the resonance is not exactly equal to the phase speed ($U/\sqrt{{gH}}=1/3$) but has a small deviation ($U/\sqrt{{gH}}\approx \,0.31$). Equation (33) shows that the amplitude of the Rossby wave is maximum at $U/\sqrt{{gH}}=1/3$, while Equation (32) shows that the amplitude of the Kelvin wave decreases as the westerly wind speed increases. Therefore, the coupled Rossby–Kelvin wave pattern resonates when $U/\sqrt{{gH}}$ is somewhat less than 1/3.

The resonance behavior as well as the phase shift of the waves can influence horizontal eddy momentum flux convergence, as shown in Figure 10. In this figure, a positive value means that the waves accelerate the mean flow and a negative value decelerate the mean flow. On the equator, the momentum flux convergence is positive as long as the mean flow is not too small and the corresponding value is negative in higher latitudes, indicating that zonal momentum is transported from the higher latitudes to the equatorial region. This acts to maintain the equatorial superrotation against friction and dissipation (Herbert et al. 2020). Markedly, there is a peak in the equatorial acceleration when the resonance between the Rossby waves and the mean flow occurs (panels (e), (f), (h), and (i)). However, not all the acceleration peaks are occurring synchronously with the resonance. For example, when τ_rad = 1 and τ_drag = 1, the acceleration peak happens when the strength of the mean flow is equal to about 0.8 of the gravity wave speed (panel (a)), and when τ_rad = 1 and τ_drag = 100 (panel (c)) and τ_rad = 100 and τ_drag = 1 (panel (g)), the maximum acceleration happens when the strength of the mean flow is close to zero. This is due to the fact that the horizontal eddy momentum flux convergence is determined by the combined condition of the amplitude and the phase of the coupled Rossby–Kelvin waves; when the wave amplitude reaches a maximum, the tilt of the coupled waves may not be in an optimal condition for the equator-ward momentum transport. Moreover, the resonance between the Kelvin wave and the mean flow does not imply a momentum convergence peak, as shown in Figure 10; this is due to that the tilt of the waves is very small or close to zero under this resonance (see the 3rd column of Figure 4). These results suggest that a coupled system with both active waves and active mean flow is required in future work.

Zoom In Zoom Out Reset image size

The phase shift of the planetary waves and the wave–jet resonance are also found in the spin-up period of a 3D atmospheric general circulation model (AGCM) experiment, as shown in Figure 11. The experiment was performed using the global climate model CAM3, same as that used in Yang et al. (2013). The stellar flux was set to 1200 W m⁻², the star temperature is 3400 K, the rotation period (=orbital period) is 37 Earth days, the surface air pressure is 1.0 bar N₂, and atmospheric CO₂ concentration is 300 ppmv. The surface is covered by a 50 m slab ocean with no any continent. This experiment was initialized from a climate state similar to the present-day Earth.

Zoom In Zoom Out Reset image size

As shown in Figure 11, the wave pattern and atmospheric superrotation are established within about 100 Earth days. Initially, the geopotential height field at 200 hPa is roughly related to the surface land-sea distributions with high values over the oceans and low values over the continents (Figure 11(a)). This is because the experiment was started from the summer atmospheric state of modern Earth, although the surface is set to be an aqua-planet in the simulation. On the sixth Earth day, a zonal-number one wave pattern is formed with high-pressure centers in the west of the substellar point and with low-pressure centers in the east of the substellar point. This wave pattern moves eastward gradually following the zonal flow, which exhibits an equatorial superrotation behavior after about 10 Earth days. The amplitude of the waves also evolves with time and reaches a peak in about 50 Earth days (Figure 11(h)) when the wave crest is right on the longitude line of the substellar point (Figure 11(e)). After that, the waves move further east and the wave amplitude decreases. The wave crest never goes across the east terminator, likely implying the existence of an upper limit (π) for the phase shift. These results confirm that the phase shift of planetary waves and the wave–jet resonance on tidally locked planets do occur in a more realistic 3D AGCM simulation. In the 3D AGCM simulations of Arnold et al. (2012), Tsai et al. (2014), and Hammond & Pierrehumbert (2018), they showed similar phase shifts, but the resonance behavior in 3D tidally locked AGCM simulations is the first time to be uncovered here.