THREE-DIMENSIONAL STRUCTURES OF EQUATORIAL WAVES AND THE RESULTING SUPER-ROTATION IN THE ATMOSPHERE OF A TIDALLY LOCKED HOT JUPITER

As has been commonly adopted in the study of equatorial waves, we work on the thermal-tide problem by utilizing Cartesian geometry on an equatorial beta plane for a hot Jupiter with the radius R_p and spin rate Ω. Any fluid quantity q in the atmosphere is linearized in terms of its zonal-mean q₀ and its perturbation associated with planetary-scale waves q'. In an atmosphere of a hot Jupiter without thermal inversion, it can be assumed without loss of generality that the zonal-mean state (i.e., the basic state) of the atmosphere is isothermal and hydrostatic (e.g., Madhusudhan & Seager 2009), with a uniform zonal flow at a speed of u₀ to resemble the planetary-scale super-rotation flow. As a result of strong vertical stable stratification as well as the fact that phase velocities of planetary waves are much less than the sound speed, we can further assume the perturbations to be vertically hydrostatic (e.g., Showman et al. 2011) and nearly incompressible such that the anelastic approximation is applicable to eliminate acoustic waves (Ogura & Phillips 1962). Hence, considering wave motions in an atmosphere of a tidally locked planet driven by stationary thermal forcing, the linearized equations on an equatorial beta plane are given by (e.g., Holton 1992)

$$\begin{equation} \alpha u^{\prime } + u_{0} \frac{\partial u^{\prime }}{\partial x} - \beta y v^{\prime } + \frac{\partial \phi ^{\prime }}{\partial x} =0, \end{equation} \tag{ 1a }$$

$$\begin{equation} \alpha v^{\prime } + u_{0} \frac{\partial v^{\prime }}{\partial x} + \beta y u^{\prime } + \frac{\partial \phi ^{\prime }}{\partial y} =0, \end{equation} \tag{ 1b }$$

$$\begin{equation} \frac{\partial u^{\prime }}{\partial x} + \frac{\partial v^{\prime }}{\partial y} + \frac{1}{\rho _{0}}\frac{\partial }{\partial z} (\rho _{0} w^{\prime }) =0, \end{equation} \tag{ 1c }$$

$$\begin{equation} \alpha \frac{\partial \phi ^{\prime }}{\partial z} + u_{0} \frac{\partial ^{2} \phi ^{\prime }}{\partial x \partial z} + N^{2} w^{\prime } = \frac{\kappa }{H} J, \end{equation} \tag{ 1d }$$

where we have already assumed that the ∂/∂t terms all vanish for steady-state solutions (which are actually treated as quasi-steady for the evolution of u₀ on a longer timescale, as will be discussed later in Section 3.3). In the above equations, x, y, z are the zonal, meridional, and vertical coordinates, and u', v', w' are the horizontal, meridional, vertical velocities, respectively. β = 2Ω/R_p is the Coriolis parameter in the equatorial beta-plane approximation, which approximates the latitudinal variation of the Coriolis force from the equator. ϕ' = p'/ρ₀ is the perturbed geopotential with p' being the perturbed thermal pressure and ρ₀ being the basic density given by

$$\begin{equation} \rho _{0}(z) = \rho _{0}(z_b)e^{-(z-z_b)/H}, \end{equation} \tag{ 2 }$$

where H is the pressure and density scale height in the isothermal atmosphere and z_b is the prescribed location of the bottom boundary. Equation (2) states that z is equivalent to the logarithmic pressure coordinate, which will be applied in this work. N is the buoyancy frequency. α is the linear momentum damping rate (a.k.a. Rayleigh friction) in the momentum equations and thermal damping rate (a.k.a. Newtonian cooling) in the energy equation. These two types of damping are set equal here for the advantage of the mode analysis shown later in this paper to better elucidate underlying physics. Furthermore, J is the stellar heating rate per unit mass in the planetary atmosphere. κ is the specific gas constant R divided by the specific heat at constant pressure c_p.

We note that u₀ is assumed to be uniform over the entire atmosphere under consideration. Although it is not self-consistent to ignore the Coriolis effect and shear of the zonal-mean zonal flow, we add a uniform u₀ here to study the effect of zonal-mean zonal flow on the global-scale waves and their momentum transport, without making the equations too difficult to solve by separation of variables (Phlips & Gill 1987; Arnold et al. 2012).

We first follow the shallow-water approach employed by Wu et al. (2000) in which the solution is separated into a vertical part and a horizontal part, where the vertical part of the solution is first expressed by the vertical eigenfunctions (the free vertical modes). For each of these vertical modes, the horizontal structure is governed by the stationary forced shallow-water equations but with a different equivalent depth respectively. The wave solutions to the forced shallow-water equations can then be expressed by free horizontal modes (i.e., equatorial waves). The summation of each vertical mode in combination with the horizontal modes as the solutions to the corresponding forced shallow-water equations are the complete solutions in the 3D problem.

In this study, we focus solely on the s = 1 response (s is the longitudinal wavenumber), which is expected to be the gravest response of the atmosphere in a tidally locked heating scenario. Following Wu et al. (2000), we eliminate w' in Equation (1) for easier decomposition in shallow-water equations later. The solutions can then be written as the sum of vertical modes Z_m(z), each of which is denoted by m:

$$\begin{equation} {\left(\begin{array}{@{}c@{}} u^{\prime }(x,y,z)\\ v^{\prime }(x,y,z)\\ \phi ^{\prime } (x,y,z) \end{array}\right)} = \sum _{m} {\left(\begin{array}{@{}c@{}}u_{m}(y)\\ v_{m}(y)\\ \phi _{m}(y) \end{array}\right)} Z_m(z) e^{i k x } e^{z/2H}, \end{equation} \tag{ 3 }$$

where e^z/2H is the density weighted factor resulting from the vertical stratification and k = s/R_p = 1/R_p is the wavenumber in the x direction for the s = 1 modes. By defining the relation describing the vertical eigenfunctions (i.e., free vertical modes) as follows

$$\begin{equation} \left(\frac{d^{2}}{d z^{2}}-\frac{1}{4H^{2}} \right)Z_{m}(z)=-\frac{N^{2}}{gh_{m}}Z_{m}(z), \end{equation} \tag{ 4 }$$

where g is the gravitational acceleration on the planetary surface and h_m is the eigenvalue known as the "equivalent depth," Wu et al. showed that Equation (1) can be transformed to sets of stationary forced shallow-water equations for the coefficients u_m, v_m, and ϕ_m of the solutions in Equation (3). That is, the set of shallow-water equations for each vertical mode denoted by m reads

$$\begin{equation} (\alpha + i k u_{0}) u_{m} - \beta y v_{m} + i k \phi _{m} = 0, \end{equation} \tag{ 5a }$$

$$\begin{equation} \beta yu_{m} + (\alpha + i k u_{0}) v_{m} + \frac{\partial \phi _{m}}{\partial y} = 0, \end{equation} \tag{ 5b }$$

$$\begin{equation} i k u_{m} + \frac{\partial v_{m}}{\partial y} + (\alpha + i k u_{0}) \frac{\phi _{m}}{gh_{m}}= -F_{m}, \end{equation} \tag{ 5c }$$

where F_m satisfies

$$\begin{eqnarray} &&Q \equiv \frac{\kappa }{N^{2}H} \left(\frac{\partial }{\partial z} - \frac{1}{2H}\right)J e^{-z/2H} = \sum _{m} F_{m}(y)Z_{m}(z)e^{{\rm i}kz}.\nonumber\\ \end{eqnarray} \tag{ 6 }$$

Namely, F_m is the projection of a quantity related to the heating function onto a vertical eigenmode m. The summation in Equations (3) and (6) represents a sum over all discrete modes or an integration over all continuous modes.

The set of Equation (5) is identical to the shallow-water equation in Matsuno (1966) by replacing α and the ocean depth with α + iku₀ and h_m, respectively. Thus, it can again be solved straightforwardly by expanding the solutions (u_m, v_m, ϕ_m) and the heating function F_m in terms of the meridional eigenfunctions, i.e., equatorial Kelvin waves, equatorial Rossby waves, and inertia-gravity waves, obtained from the free equatorial wave equations subject to the zero-velocity boundary condition as done in Matsuno (1966). It follows that the characteristic speed, length, and time for each vertical m mode in Equation (3) are given by $C_m=\sqrt{gh_{m}}$, $L_m=(\sqrt{gh_{m}}/ \beta)^{1/2}$, and $T_m=(\sqrt{gh_{m}} \beta)^{-1/2}$, respectively. L_m is the so-called equatorial radius of deformation describing the characteristic meridional scale of equatorial trapped waves and T_m is the characteristic wave travel time over L_m with the characteristic speed (i.e., horizontal phase speed) of gravity wave $\sqrt{gh_m}$. We can see from these expressions that a vertical m mode with a larger h_m has a larger horizontal phase speed ∼C_m and a larger horizontal scale ∼L_m.

The expansion procedure leads to an expression describing how the atmosphere response, designated by a, reacts linearly to the thermal forcing, represented by b, in terms of each vertical mode (m) and each type of the corresponding meridional eigenfunction³ (n, l) denoted together by (m, n, l) in the presence of a linear damping and a zonal-mean flow. In other words, we have (see Equation (34) in Matsuno 1966)

$$\begin{equation} a_{m,n,l} = \frac{1}{\alpha - i (\omega _{m,n,l} - \omega) } b_{m,n,l}, \end{equation} \tag{ 7 }$$

where ω is the stationary forcing frequency Doppler-shifted by the zonal-mean zonal flow (i.e., ω = u₀k) and ω_{m, n, l} is the eigenfrequency of the free-wave mode (m, n, l). For sign convention in this study, we always choose k to be positive. Consequently, positive (negative) ω and ω_{m, n, l} represent the horizontal phase velocity of the thermal forcing and free modes going westward (eastward). Note that the simple relation between the atmospheric response and thermal forcing in Equation (7) results from the same damping rate applied to both the Rayleigh friction and Newtonian cooling (Matsuno 1966). While this procedure limits the parameter space in which the Prandtl number equals just one (see SP), its advantage is that it is easy to understand how an atmosphere responds to thermal forcing in terms of different contributions from eigenmodes. Equation (7) neatly tells us that the atmospheric response of each mode, which is proportional to the heating projection, is mainly determined by the relative magnitude of damping rate and the Doppler-shifted free-wave frequency as viewed by the uniform zonal-mean flow. Specifically, the response amplitude in Equation (7) can be rewritten in the polar form

$$\begin{equation} a_{m,n,l} = r_{m,n,l} \, \exp ({\rm i}\theta _{m,n,l}), \end{equation} \tag{ 8 }$$

where

$$\begin{eqnarray} r_{m,n,l} &=& \frac{b_{m,n,l}}{ \sqrt{\alpha ^2 + (\omega - \omega _{m,n,l})^2}}, \nonumber\\ \theta _{m,n,l} &=& \arctan \left(\frac{\omega _{m,n,l} - \omega }{\alpha } \right). \end{eqnarray} \tag{ 9 }$$

The phase difference in θ_{m, n, l} is crucial in introducing the correlation of horizontal velocities. For example, the horizontal momentum flux contributed from the modes with the same m and n is given by

$$\begin{eqnarray} (\overline{u^{\prime }v^{\prime }})_{m,n} &=& {1\over 2}{\rm Re} \sum _{l,l^{\prime }} [ u_{m,n,l} v^*_{m,n,l^{\prime }}]\nonumber\\ &=& {1\over 2}r_{m,n,l}\,r_{m,n,l^{\prime }} \sum _{l,l^{\prime }}{\rm sin} (\theta _{m,n,l}-\theta _{m,n,l^{\prime }}) |u_{m,n,l}||v_{m,n,l^{\prime }}|,\nonumber\\ && \end{eqnarray} \tag{ 10 }$$

where the overbar means taking a zonal mean and * means taking a complex conjugate. In deriving the above relation, we have used the same convention as Matsuno (1966) for free-wave eigenfunctions in which u_{m, n, l} is real and v_{m, n, l} is pure imaginary. It is evident from the above equation that the term "sin" associated with the phase difference between two modes provides the correlation between u' and v', giving rise to the horizontal momentum flux. For instance, when l = l', there is no correlation from one single horizontal mode and the momentum transfer is zero. Therefore, different horizontal modes are necessary in order to create the phase differences and thus the horizontal momentum flux, as the roles played by the combined equatorial Kelvin and Rossby waves in SP.

Finally, we stress that even though the phase speed of Kelvin waves is not a function of k, it is not constant because h_m changes with m. In this respect, baroclinic Kelvin waves with vertical propagation are "dispersive" in a vertically stratified atmosphere. Likewise, baroclinic Rossby waves are not only dispersive in the horizontal direction but are also "dispersive" in the same way as baroclinic Kelvin waves due to different vertical modes denoted by m. We shall discuss the dispersive effect due to vertical stratification in Section 3.2.3.

To solve Equation (4) for eigenfunctions and eigenvalues, boundary conditions at the top and bottom boundary of the atmosphere of interest need to be specified. Since we shall compare our 3D linear analysis with our numerical results based on the 3D radiative hydrodynamical simulation in DA and the 3D GCM in SP, we employ the free-slip bottom boundary condition as adopted by DA, namely, the derivatives in the z direction of u' and v' at the bottom boundary are equal to zero, which gives one constraint:

$$\begin{equation} \frac{d Z_{m}(z)}{d z} + \frac{Z_{m}(z)}{2H} = 0,\quad {\rm at } z=z_{b}. \end{equation} \tag{ 11 }$$

This bottom boundary condition is similar to the impermeable condition used by SP for their circulation simulation.

Moreover, we use the outgoing radiation condition as z → ∞. It should be noted that the free-slip bottom boundary in our linear model reflects waves rather than allowing waves to travel downward according to Equation (1d). It is not reasonable unless the waves propagating to z = z_b have been strongly damped.

The eigenfunctions satisfying the boundary conditions consist of a single barotropic mode

$$\begin{equation} Z_{m}(z) = e^{-z/2H} \end{equation} \tag{ 12 }$$

and the continuous baroclinic modes (see Wu et al. 2000)

$$\begin{equation} Z_{m}(z) = \frac{A {\rm sin}\, mz + B {\rm cos}\, mz}{ \sqrt{A^{2}+B^{2}}}, \end{equation} \tag{ 13 }$$

where

$$\begin{eqnarray*} &&A = 2 H m \sin (m z_b) - \cos (m z_b), \\ &&B = 2 H m \cos (m z_b) + \sin (m z_b), \end{eqnarray*}$$

with $m \equiv \sqrt{N^2/gh_m - 1/4H^2}$ being any real positive number, representing the vertical wavenumber for the continuous baroclinic mode.

For mathematical simplicity, we assume the heating rate Q in Equation (6) has the separable form p(z)S(y)e^ikx like the heating function used in Wu et al. (2000), and then we rewrite

$$\begin{equation} F_{m} = F(m)S(y). \end{equation} \tag{ 14 }$$

Hence

$$\begin{equation} \sum _{m} F(m)S(y)Z_{m}(z) = p(z)S(y), \end{equation} \tag{ 15 }$$

where

$$\begin{equation} F(m)= \sqrt{\frac{2}{\pi }} \int _0^\infty p(z) Z_{m}(z)\, dz. \end{equation} \tag{ 16 }$$

By calculating F(m), the projection of p(z) onto the vertical eigenfunction Z_m(z), we obtain the contribution from each vertical m mode. As shown in the preceding section, each vertical mode has a horizontal structure defined by a set of shallow-water equations with equivalent depth h_m and corresponding length scale and eigenfrequencies. We sum them up to construct the complete solutions.

The heating function J in Equation (6) describes the strength and distribution of the thermal forcing due to the stellar irradiation. In our linear model, the zonally uniform heating (i.e., s = 0) provides the equilibrium temperature, which is used as the background temperature T₀ of the isothermal atmosphere. As has been stated in the paper, we focus ourselves on the s = 1 thermal forcing in the present work.

Figure 1 shows the s = 1 distribution of the stellar heating rate per unit volume absorbed in the atmosphere, i.e., ρ₀J, adopted in the linear analysis (left column) in comparison with that in the radiative hydrodynamic simulation (right column). We apply a simple analytical function to mimic the vertical profile of s = 1 thermal forcing in Dobbs-Dixon et al. (2010) as shown in the top row of Figure 1; namely,

$$\begin{equation} \rho _{0}(z) J(z) =\rho _0(z) F_{*,s=1} \kappa _e \exp {[\kappa _{e} H(\rho _{0}(z_{t})-\rho _{0}(z))]}, \end{equation} \tag{ 17 }$$

where F_{*, s = 1} is the s = 1 component of the stellar flux, κ_e is a parameter presenting the opacity effect, and z_t is the vertical location of the top boundary. The vertical coordinate is set to zero at 1 bar, while the bottom boundary where the free-slip boundary condition applies lies at the pressure level of 20 bars.⁴ Hereafter, whenever a vertical structure of the atmosphere is presented, we show the results in the vertical domain from the bottom boundary to the altitude at 3 × 10⁻⁴ bars. κ_e is set as 0.0175 cm² g⁻¹ for a close fit to the vertical heating profile in the simulation, as shown in the top panels of Figure 1. The heating-center level of the stellar heating rate per unit volume lies at p ≈ 0.1 bar. Most of the stellar irradiation is absorbed between 0.1 and 1 bar.

Zoom In Zoom Out Reset image size

The meridional heating profiles in the linear theory are assumed to be Gaussian functions (see the middle left panel for ρJ at the substellar point where x = 0), with location and width matching the stellar heating term used in the simulation (see the middle right panel for ρJ at ϕ = 0).⁵ The slight difference between the latitudinal heating profile in the linear theory and in the simulation is due primarily to the fact that the atmosphere in the simulation is slightly colder and thus vertically thinner at higher latitudes, whereas the background atmosphere in the linear theory is assumed independent of latitudes. As illustrated in the figure, the width of the intense stellar heating in the y direction is large scale, comparable to the equatorial Rossby deformation radius, as previously noted by SP. A sinusoidal variation along the longitude, i.e., exp (ikx), is applied in the linear theory to represent the day–night heating contrast; ρJ at the equator shown in the bottom left panel is similar to the s = 1 thermal forcing in the simulation plotted in the bottom right panel.

Since the characteristic length scale L_m of each vertical m mode now differs, it is not possible to have one set of free-wave eigenfunctions (n = ±1 modes only) with the meridional scale (equatorial Rossby deformation radius) matching the heating profile, such that the solutions can be expanded with only a few terms like those in Matsuno (1966), in Gill (1980), or in SP. Theoretically, infinite terms are required in order to expand the meridional Gaussian heating. In practice, higher n modes contribute much less at the equator so the solutions can be truncated. Because the eigenfunctions with even n are anti-symmetric to the equator, n should be odd numbers to provide the eigenfunctions for constructing a stellar heating symmetric to the equator. Modes up to n = 7 are used to get a fairly good approximation in our analysis. The presence of higher n modes because of the vertical structure decomposition does not alter the results significantly in most cases, compared with the shallow-water solutions, but differences arise when the damping term is small, which will be further discussed in Section 3.1.

The power spectrum F²(m) corresponding to the vertical heating in the linear theory is presented in Figure 2. As described previously, the energy spectrum of the baroclinic mode should be continuous. Nevertheless, the 20 discrete baroclinic modes shown in Figure 2 are able to construct a reasonably good profile of the vertical stellar heating. To further confirm that we do not terribly miss any dominant modes, the power spectrum composed of 200 vertical modes are also calculated and illustrated in Figure 2, showing the spectral peaks agree well with that constructed only by 20 modes. In what follows, we shall use the 20 vertical modes for the linear calculation. A baroclinic mode with a smaller m has a longer vertical wavelength, and therefore has a larger equivalent depth h_m and a larger meridional scale ∼L_m. The projection on the barotropic mode is not shown but it is negligible. Generally, only a very deep heating profile that extends to the bottom boundary would give rise to noticeable projection on the barotropic mode.

Zoom In Zoom Out Reset image size

To roughly examine the sensitivity of the wave solutions to the vertical heating structure, we have further tried various values for the opacity parameter κ_e in Equation (17). The results for different κ_e are generally similar to each other except that the vertical location of wave patterns shifts. It is because a larger (smaller) κ_e presents an atmosphere absorbing the stellar heating at a higher (lower) level. We have also tested the linear theory with general heating profiles using Gaussian functions with different vertical widths, and found that a narrower heating profile projects more onto shorter vertical modes as expected. Having performed these tests, it should be kept in mind that the vertical heating profile, the efficiency of re-radiated radiation energy, and the background states are determined by the opacity of the gas. Modifying one without the other is inconsistent. Therefore, for the purpose of the consistency that has been taken care in our radiative hydrodynamical simulations, we shall focus on the particular case of HD 189733b as the fiducial model in this study.

Equation (7) describes the response amplitude of each mode, which reaches a maximum when ω = ω_{m, n, l}. In the limit of no dissipation (α = 0), the response grows without bound in the linear theory, which implies a resonance analogous to the singularity of the topographic Rossby waves in the presence of a westerly wind (Holton 2004). The resonances result from horizontal discrete modes; as ω → ω_{m, n, l}, the free mode becomes stationary relative to the heating in the reference frame of the zonal-mean flow and will be amplified until nonlinearity becomes important. A more detailed analysis on resonances in both linear theory and multi-level GCMs, as well as the positive feedback leading to transition to super-rotation can be found in Arnold et al. (2012) in the study of Earth's atmosphere.

The response amplitude of each mode, mainly equatorial Rossby and Kelvin waves in the range of the mean-flow speed that we investigate, is represented in terms of ∑_m|∑_na_{m, n, l}| and is shown in Figure 3. Given a u₀, the results are calculated by summing over all m and n for each type of equatorial wave. The figure is plotted using the atmospheric parameters of the hot Jupiter HD 189733b as an example (see Section 3 for the detailed information of parameters and damping rates α). In the plot for α = 0.1, we can see that the resonances of Rossby and Kelvin waves occur when the velocity of the zonal-mean flow is close to the horizontal phase velocity (ω_{m, n, l}/k) of each free wave (u₀ ≈ −1000 to −2000 m s⁻¹ for the Kelvin wave and u₀ ≈ 600 m s⁻¹ for the Rossby wave). When α is small (=0.01), the resonance is the strongest since the frequency term outweighs the damping term in Equation (7). Thus the resonance occurs at numerous vertical modes associated with various horizontal phase velocities as explained in Section 2.1. In general, smaller-scale modes have smaller equivalent depths h_m, and thus smaller phase speeds. Therefore, multiple resonance peaks are present in the top left panel. In contrast, stronger damping rates suppress the resonant effects more significantly, i.e., a large α diminishes the frequency term in Equation (7). This is illustrated in the bottom left panel for α = 0.5.

Zoom In Zoom Out Reset image size

Generally, owing to the Rossby-wave resonance, the response amplitude of the Rossby wave is greater than that of the Kelvin wave when u₀ > 0 and when the damping is weak or modest (α = 0.01 or 0.1). This is illustrated in the top two panels in the right column of the figure. In contrast, wave resonances are substantially suppressed by strong damping, leading to the response amplitude of each mode simply proportional to the thermal forcing amplitude of the corresponding mode. The Kelvin mode couples better with the assumed meridional heating distribution because similar to the heating profile, the Kelvin mode also varies in latitude as a Gaussian function centered at the equator. It follows that the response amplitude of the Kelvin wave becomes stronger than that of the Rossby wave when the damping is strong (α = 0.5). This is demonstrated in the bottom panel of the right column in the figure.

Several early theoretical studies have investigated interactions between a mean flow and waves in Earth's atmosphere. One of the well-known examples is the quasi-biennial oscillation (QBO), the alternating zonal wind in a period about 28 months in the tropical stratosphere, driven by the propagation and dissipation of equatorial Kelvin and Rossby-gravity waves. The interaction of the waves and the zonal-mean flow is deemed to be responsible for the QBO (e.g., Holton & Lindzen 1972; Andrews et al. 1987). In the case of hot Jupiters, the formation of atmospheric super-rotation can be dominantly attributed to the Reynolds stress associated with planetary-scale equatorial waves driven by stellar irradiation as explained by SP. Omitting vertical and meridional mean flows and ignoring the friction acting on the mean flow, the acceleration of a zonal-mean zonal flow is governed by (e.g., Andrews & McIntyre 1976)

$$\begin{equation} \frac{ \partial \overline{u} }{\partial t} = - \partial _{y} \overline{u^{\prime } v^{\prime }} - \frac{1}{\rho _{0}} \partial _{z} \rho _{0} \overline{u^{\prime } w^{\prime }}. \end{equation} \tag{ 18 }$$

The first and second terms on the right-hand side are the negative divergence of the horizontal and vertical components of the wave momentum flux, respectively. By definition, $u_0=\overline{u}$. Because the Reynolds stress arises from the linear theory to second order, Equation (18) describes a slow evolution of u₀ on a timescale much longer than the relaxation time of the atmosphere in response to the change of u₀ in a linear wave problem. The different timescales in the linear theory enables us to estimate the acceleration of the zonal flow based on the steady-state wave solutions for a given u₀. The relaxation time of the atmosphere will be defined and discussed in Section 3.3.2 where the implication to the evolution of u₀ is presented.

The divergence of the wave momentum fluxes in Equation (18) can be related to thermal forcing and dissipation. More specifically, combing the momentum, continuity and energy equations in Equation (1), we first obtain the zonal-mean relation between the energy fluxes, thermal forcing, and dissipation for stationary forced waves as follows

$$\begin{eqnarray} &&\frac{\partial }{\partial y}\overline{\phi ^{\prime } v^{\prime }} + \frac{\partial }{\partial z}\overline{\phi ^{\prime } w^{\prime }} - \frac{\overline{\phi ^{\prime } w^{\prime }}}{H} = -{\alpha \over {2}}(u^{\prime 2}+v^{\prime 2})\nonumber\\ &&\quad + \frac{1}{N^{2}}\left[ \overline{Q^{\prime } \frac{\partial \phi ^{\prime }}{\partial z}} - {\alpha \over {2}} \left(\frac{\partial \phi ^{\prime }}{\partial z} \right)^{2} \right], \end{eqnarray} \tag{ 19 }$$

where Q' ≡ (κ/H)J and the square of a complex number denotes the square of the amplitude of the number. The wave energy fluxes on the left-hand side of the above equation can then be replaced by the wave momentum fluxes through the procedure similar to the Eliassen-Palm first theorem (Lindzen 1990), which can be derived from the linearized momentum equation. Namely,

$$\begin{equation} \overline{\phi ^{\prime } v^{\prime }} = - u_{0}(\overline{u^{\prime } v^{\prime }}) - \overline{(\epsilon u^{\prime }) v^{\prime }}, \end{equation} \tag{ 20 }$$

and

$$\begin{equation} \overline{\phi ^{\prime } w^{\prime }} = - u_{0}(\overline{u^{\prime } w^{\prime }}) - \overline{(\epsilon u^{\prime }) w^{\prime }} + \overline{(\beta _{i} v^{\prime }) w^{\prime }}, \end{equation} \tag{ 21 }$$

where ≡ α/(ik) and β_i ≡ βy/(ik). Substituting Equations (19)–(21) into Equation (18) for the divergence of wave momentum fluxes,⁶ we obtain

$$\begin{eqnarray} {\partial u_0 \over \partial t}&=&{1\over u_0}\left[ {-}{\alpha \over {2}} (u^{\prime 2}+v^{\prime 2})+{\overline{Q^{\prime }\partial _z \phi ^{\prime }}-(\alpha / {2}) (\partial _z \phi ^{\prime })^2 \over N^2}\right.\nonumber\\ &&\left. +\partial _y \overline{(\epsilon u^{\prime })v^{\prime }}+\partial _z \overline{(\epsilon u^{\prime })w^{\prime }}-\partial _z\overline{(\beta _iv^{\prime })w^{\prime }}-{\overline{(\epsilon u^{\prime })w^{\prime }} \over H} \right].\nonumber\\ && \end{eqnarray} \tag{ 22 }$$

The right-hand side of the above equation contains the terms involving the thermal forcing ($\overline{Q^{\prime }\partial _z \phi ^{\prime }}$), dissipation (terms associated with α), and Coriolis effect (the term with β_i). Once these second-order terms are balanced out, the net acceleration of the zonal-mean zonal flow is zero and thus there exist equilibrium solutions in Equation (18) for the atmosphere. The stability of the equilibrium solutions is relevant to the evolution of u₀, which will be explored in Section 3.3.

We begin by investigating the case with no zonal-mean zonal flow. We show our solutions at 0.03 bar pressure level to compare with the linear shallow-water results in SP. The three panels in the left column of Figure 4 display the solutions on the x–y plane at the pressure level of 0.03 bar from the 3D linear model in the three damping scenarios. For comparison, the three panels in the right column of the figure are solutions obtained from linear shallow-water equations with no vertical structure, similar to the results with τ = 100, τ = 10 and τ = 1 in Figure 3 of SP (corresponding to the dissipation time constants of 1 month, 3 Earth days, and 6 hr), where τ = 1/α. We merely elaborate the change of the zonal phase of wave responses relative to the thermal forcing as α changes, which is not mathematically analyzed by SP.

Zoom In Zoom Out Reset image size

We inspect the two extreme damping scenarios first. In the strong damping scenario, according to Equation (7), the response amplitude is dictated by the factor 1/α, so the response of each mode almost has the same phase pattern as the perturbed heating profile. In other words, the geopotential response simply resembles the horizontal heating structure in our 3D linear theory (or the radiative equilibrium height in the linear shallow-water model). In the weak damping scenario, the response amplitude is instead dominated by the factor −1/iω_{m, n, l} in Equation (7). In addition, the Rossby-wave response is stronger than other waves in the weak damping scenario (see Figure 3). Hence, the geopotential response is mainly manifested by planetary-scale Rossby gyres, which have a 90° phase difference westward to the horizontal heating due to the imaginary term −1/iω_{m, n, l}.

In the modest damping scenario, α and ω_{m, n, l} have the same order of magnitude. The responses from Kelvin and Rossby waves are comparable in magnitude (see Figure 3). Since ω_Kelvin and ω_Rossby have opposite signs (eastward for Kelvin wave and westward for Rossby waves), the two types of waves with opposite propagating directions have competing magnitude, thus leading to the formation of the Matsuno–Gill pattern, described as a "chevron" pattern of velocity fields (northwest–southeast tilts in the northern hemisphere and southwest-northeast tilts in the southern hemisphere) in SP. This particular pattern, producing $\overline{u^{\prime }v^{\prime }} < 0$ in the north of the equator and $\overline{u^{\prime }v^{\prime }} > 0$ in the south of the equator, is suggested by SP (also referred to Equation (10)) to be crucial for equatorward momentum transport that induces super-rotation.

The above results based on the 3D linear model with u₀ = 0 agree fairly well with the solutions from the linear shallow-water model as shown in the panels in the right column of Figure 4, and are therefore in reasonable agreement with the analytical solutions when the drag and radiative time constants are equal as shown in Figure 3 of SP.

Having shown that the solutions around the heating-center level from the 3D linear equations are generally consistent with the shallow-water solutions, there actually exists a slight difference in the weak damping scenario. Unlike the gyres shown in the top right panel of Figure 4, the gyres from 3D solutions exhibit a minor eastward tilt in the longitudinal direction, as can be seen from the shape of the gyres in the top left panel, guided by the green arrows which roughly point to the local flow directions. In fact, the eastward tilt in the presence of vertical stratification is a result of the combination of Rossby waves summed over all vertical modes, which gives rise to higher orders of n, as discussed in Sections 2.1 and 2.3. The phenomenon is illustrated in Figure 5. Rossby waves of higher n have wider latitudinal extension and slightly slower phase velocity (Matsuno 1966), such that they have smaller westward shifts than the Rossby waves of lower n. Referring to Equation (7) again, with a small damping rate and no mean flow, we can see that Rossby modes with higher n have smaller westward shifts. The eastward tilt is not seen in the linear shallow-water model (e.g., top right panel of Figure 4 or Figure 3 in SP) since there are no n > 1 waves in the shallow-water model forced by the Gaussian heating function with the same latitudinal extent as the equatorial Rossby deformation radius. For the same horizontal heating profile, vertical stratification introduces the eastward tilt at high latitudes.

Zoom In Zoom Out Reset image size

We have demonstrated that the atmospheric responses based on our 3D linear theory around the heating-center level in the absence of a zonal flow qualitatively agree with the 2D linear analysis in SP. Now we extend our linear analysis to study one of the key issues of our paper, namely, to examine the atmospheric responses in the presence of a uniform eastward mean zonal flow.

3.2.1. Horizontal Structures Around the Heating-center Level: Eastward Shift of Rossby Modes as the Eastward Zonal Flow Speed Increases

Figure 6 compares the linear solutions with the steady-state simulation results around the heating-center level (p = 0.03 bar for both of the linear model and simulation). The figure shows the wave velocity fields superimposed on the color-coded geopotential contours in the cases of different damping rates. We adopt the typical zonal-mean flow speed 1000 m s⁻¹ for u₀ in the linear model. The simulation results shown in the figure are Fourier transformed for the s = 1 component.

Zoom In Zoom Out Reset image size

Despite the fact that the nature of viscous dissipation is different from linear damping and that the equatorial super-rotating velocities are not the same in the numerical simulations with different viscosities, we present the simulation results for three different viscosities progressively differing also by a factor of 10 to qualitatively compare with the linear results for the three damping scenarios. Figure 6 shows that the wave velocity at higher latitudes are weaker in the linear theory than those in the numerical simulations. It is because the linear solutions with an a prior global uniform zonal flow on the equatorial β-plane decline with latitude due to the zero-velocity boundary condition (Matsuno 1966; Holton 2004), whereas in the atmospheric simulations in the spherical geometry, there actually exist flows across the polar regions and hence no global zonal flows are present at high latitudes. With all the uncertainties and discrepancies, Figure 6 nonetheless shows that the linear solutions resemble the overall structure around the heating-center level in the numerical simulation.

We further compare Figures 4 and 6 to investigate the changes of atmospheric responses after a eastward zonal flow is developed. As shown in the figures for weak and modest damping scenarios, the solutions with a zonal flow exhibit cyclones and anticyclones like the case for u₀ = 0 but with almost opposite phases. In the weak damping scenario shown in the upper panels of Figure 6, since now a_{m, n, l} is predominately determined by the factor −1/i(ω_{m, n, l} − ω) in response to the thermal forcing, a sufficiently large ω = u₀k (i.e., ω > ω_{m, n, l} due to the large u₀) leads to a sign change and thus the opposite phase relative to the no-zonal-mean-flow case. In addition, the pressure variations increase at the equator compared to the no-zonal-mean-flow case in the weak damping scenario. The reason is that a weak frictional force only requires a small pressure gradient to balance out in the absence of Coriolis forces and a zonal flow at the equator. When u₀ is nonzero, the pressure gradient consequently increases to maintain the force balance against advection flows. In the modest damping scenario, the middle panels of Figure 6 show that the structure not only exhibits a phase shift eastward to the structure for the non-zero-mean flow case (refer to Equation (7)), the Rossby component also exceeds the Kelvin component due to the Rossby resonance (see Figure 3). Therefore the chevron pattern is no longer present once the zonal flow is established. In the strong damping scenario, α is significantly large, and the presence of u₀ only causes a slightly eastward phase shift as shown in the bottom panels of Figure 6 in comparison with the bottom panels of Figure 4. In this strong damping regime, the wave velocity deviates significantly from the geostrophic configuration and generally flows from the day to night side, driven primarily by the pressure gradient with a slight influence of the Coriolis force.

Figure 7 further presents the solutions with increasing speed of the zonal-mean flow from 0, 500, to 1000 m s⁻¹ in the modest damping rate scenario. We decompose the solutions (left column) in terms of Kelvin (middle column) and Rossby (right column) components As indicated by the strength of wave geopotentials, waves in the case of u₀ = 500 m s⁻¹ displayed in the middle row have the strongest response due to the resonant effect. When u₀ = 0, the geopotential maximums of Kelvin component and Rossby component have a phase difference by about 90°, forming the chevron pattern (i.e., $\sin (\theta _{m,n,l}-\theta _{m,n,l^{\prime }})\approx 1$ in Equation (10)). However, the phase of Rossby component progressively shifts eastward with increasing u₀, while that of Kelvin component remains almost stationary. This is because the pattern speeds of the Kelvin component are generally too large to be affected by the change of u₀ smaller than 1000 m s⁻¹. As a result, when the mean flow reaches a certain speed, say u₀=1000 m s⁻¹ in this example, the maxima of Kelvin component and Rossby component have almost the same zonal phase. Thus, the chevron pattern progressively weakens as the speed of the zonal-mean flow increases, superseded by more symmetric Rossby gyres. This evolution can be found in Figure 3 of Arnold et al. (2012) as well. Therefore the structure provides weaker horizontally equatorward momentum transport as a result of the smaller zonally averaged flux $\overline{u^{\prime }v^{\prime }}$, in agreement with Equation (10) with $\sin (\theta _{m,n,l}-\theta _{m,n,l^{\prime }})\approx 0$. The weakening of the zonal acceleration due to a change of the equatorward momentum transport after the super-rotating flow is developed was first mentioned in SP. The authors reported that the Kelvin-wave structure is shifted eastward and the midlatitude velocity structure differs significantly in the final equilibrated state by observing the morphology of the temperature and velocity fields in their simulations. Although it seems to appear in that way, our linear-mode analysis suggests that the same change of the atmospheric structure near the equator is in fact attributed to the eastward shift of the Rossby wave with almost no zonal phase shift of the Kelvin wave. The topic of momentum fluxes will be discussed in more detail in Section 3.3.

Zoom In Zoom Out Reset image size

3.2.2. Vertical Structures on the x–z Plane: Baroclinic Equatorial Waves and Wavefront Tilt

In the section, we study the vertical structures of equatorial waves projected onto the x–z plane along the equator, primarily focusing on the tilt of wavefronts (constant phase lines) of baroclinic equatorial waves, which do not appear in the shallow-water analysis.

Figure 8 displays the vertical wave structures on the x–z plane corresponding to those in Figure 6 at the equator. ϕ' superimposed with velocities u' and w' are plotted for three different damping scenarios in the linear model to compare to our simulation results with three viscosities. Once again, u₀ = 1000 m s⁻¹ is adopted in the linear model to represent the final equilibrium state in which the equatorial jet has fully developed. Inspection of Figure 8 indicates that the linear solutions qualitatively resemble the simulation results. The nature of the equatorial waves may be identified from these plots. Despite the Coriolis effect, the vertical structures of baroclinic equatorial waves appear similar to those of internal gravity waves according to dispersion relations (e.g., Holton 1992). Since equatorial Kelvin modes have eastward horizontal phase velocities and Rossby modes have westward horizontal phase velocities, such properties along with those for internal gravity waves can determine the tilt of wavefronts in their vertical structures. As the baroclinic waves propagate vertically downward from the level of excitation, the wave fronts, which are perpendicular to the wave vectors, exhibit an eastward or westward tilt depending on the directions of their horizontal phase velocities. As indicated by the ϕ' contours in Figure 8, the wavefronts tilt westward with increasing depth from the heating level when α = 0.01 and 0.1 but tilt eastward with increasing depth from the heating level when α = 0.5.

Zoom In Zoom Out Reset image size

The difference in the vertical structures arises because the Rossby-wave response is stronger than the Kelvin-wave response in both the weak and modest damping scenarios but weaker in the strong damping scenario, as has been explained earlier in the paper. To assert the relation between the tilt of wavefronts and the types of waves, we decompose ϕ' into the Rossby and Kelvin components ($\phi ^{\prime }_K$ and $\phi ^{\prime }_R$, respectively) in the moderate and strong damping scenarios and plot them in Figure 9. The same decomposition is performed for the velocity field. The figure confirms that the wave front of $\phi ^{\prime }_R$ below the heating level exhibits a downward–westward tilt and that of $\phi ^{\prime }_K$ exhibits an downward–eastward tilt. $\phi ^{\prime }_R$ is larger (smaller) than $\phi ^{\prime }_K$ in the modest (strong) damping scenario. Thus, the wave front of ϕ' (${=}\phi ^{\prime }_R+\phi ^{\prime }_K$) below the heating level preferentially exhibits a downward–westward tilt in the modest damping scenario and a downward–eastward tilt in the strong damping case, as shown in Figure 8 for α = 0.1 and α = 0.5, respectively.

Zoom In Zoom Out Reset image size

u' and w' also exhibit the same wavefront tilts as ϕ' but with different phase differences depending on the degree of dissipation, albeit difficult to recognize by arrows in Figure 8. The horizontal velocity u' around the heating-center level agrees with that shown in Figure 6. The velocity field is almost horizontal since the vertical velocity is relatively small due to strong vertical stratification and hence is hard to see in the figure. Nevertheless, upwelling (downwelling) flows primarily occur in the regions where the horizontal westward (eastward) flows take place. Therefore, the presence of small w' leads to an overall upward–westward and downward–eastward velocity correlation, which results in downward transfer of momentum and will be discussed in Section 3.3.1.

We further present the linear results in Figure 10 for the vertical structures of waves projected on the vertical plane at equator for u₀ = −1000, 0, 500, 1000 m s⁻¹ in the modest damping scenario, as have been done for the horizontal structures of waves shown in Figure 7. Figure 7 shows that when u₀ = 0, the horizontal chevron pattern leads to more symmetric horizontal flows from the day to night sides near the equator. As u₀ increases and therefore the Rossby-wave component shifts eastward, the day-to-night horizontal flows in the equatorial region become less symmetric. This trend of the flow changing with u₀ is also seen in vertical wave structures plotted in Figure 10. As indicated from the color-coded geopotential contours, the wave response to the stellar heating for u₀ = 500 m s⁻¹ is more pronounced than in other cases, and so is the magnitude of the wavefront tilt from the vertical in the heating region (p ≲ 0.1 bar), which results from the Rossby-wave resonance. The resonance amplifies the wave response so that Rossby waves can propagate farther before being dissipated (see Gill 1980). This may also be roughly identified by the WKBJ dispersion relation for internal gravity waves with the vertical wavenumber k_z,

$$\begin{equation} \frac{(\omega - u_0 k)^2}{k^2} \approx \frac{N^2}{k^2_z}, \end{equation} \tag{ 23 }$$

where the damping rate is omitted for a simple illustration. When a resonance occurs, ω − u₀k is small, leading to large k_z while k and N are fixed. Thus the larger k_z brings the wavefront tilt more from the vertical in the heating region. We also show in the upper left panel that the largest wavefront tilt for u₀ < 0 occurs near the Kelvin-wave resonance (i.e., u₀ ≈ −1000 m s⁻¹) in the heating region, in support of the relation between the wavefront tilt angle $\arctan (k_z/k)$ and u₀.

Zoom In Zoom Out Reset image size

Because the bottom boundary of the linear model reflects waves and its location at p = 20 bars eliminates the modes with longer vertical wavelengths, we examine the issue by extending the location of the bottom boundary down to 200 bars, while keeping the background states the same for a simple test. We find that the results in the original vertical domain (i.e., p ⩾ 20 bars) do not change significantly. This arises because the waves are substantially dissipated beneath the heating region. Equation (19) indicates that the effect of wave dissipation with the altitude in the absence of thermal forcing can be examined by inspecting the profile of the vertical wave energy flux $\overline{p^{\prime }w^{\prime }}$.⁷ We integrate the wave energy fluxes over all latitudes for the three damping scenarios and plot them as a function of z in Figure 11. As expected, the energy fluxes derived based on two different locations of the bottom boundary agree better for stronger damping (i.e., lager α). Nonetheless, the flux profiles approximately coincide in each damping scenario except for the region near the bottom boundaries where the vertical flux profiles decline quickly to meet the free-slip boundary condition. Overall, the vertical wave energy flux decreases by one to two orders of magnitude from the altitude just below the heating region (p ≳ 1 bar) to p ≲ 20 bars, inferring that with the bottom boundary at 20 bars, the amplitude of the reflected waves is small enough to not affect the wave solutions significantly. Specifically, the vertical group velocity for the Rossby modes is given by

$$\begin{equation} \frac{\partial \omega _{\rm Rossby}}{\partial m} \approx - \omega _{\rm Rossby}^{2} \frac{m(2n+1)}{ N k \sqrt{m^2+1/4H^2}}. \end{equation} \tag{ 24 }$$

For n = 1 Rossby waves, the vertical group velocity is about 0.05–1 m s⁻¹. Thus, the estimated vertical wave travel time crossing the vertical domain from the heating-center level to the bottom boundary is roughly 10⁶ s, which is comparable to the small damping timescale, and longer than the moderate and strong damping timescales.

Zoom In Zoom Out Reset image size

3.2.3. Horizontal Structure Below and Above the Heating Level: Dispersive Phenomenon Induced by Baroclinic Modes

Figures 12 and 13 provide another view to look at the tilted wavefronts, showing horizontal wave geopotentials overlaid with projected velocity fields at different depths obtained from the linear theory (left column) and from the numerical simulation (right column). The wave structures in the modest damping scenario are shown in Figure 12 to compare to those in the strong damping scenario displayed in Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected from the preceding section, the Rossby gyres around the equator are gradually shifted westward at greater depths. We can see from Figure 12 that on the wider horizontal scale, the Rossby gyres are also gradually tilted westward as going deeper, but interestingly there are smaller-scale gyres emerging from inside of the outer gyres around 1 bar, right below the bottom edge of the heating. The phenomena are seen in both the linear model and simulation. This change of horizontal patterns at various depths is also present in Figure 6 of Kraucunas & Hartmann (2005) where the authors study zonal winds in the tropical upper troposphere of the Earth.

This particular pattern is a result of dispersion of a group of equatorial waves of different scales. As explained earlier in the paper, baroclinic Rossby modes dominate over baroclinic Kelvin modes in the modest damping scenario, which permits us to consider Rossby waves only to interpret the results in Figure 12 without loss of generality. Baroclinic Rossby waves with smaller vertical wavelengths have smaller equivalent depths h_m, and thus have smaller latitudinal scales (see Section 2.1). Rossby modes of smaller scales would tilt more horizontally compared to those of larger scales, meaning that small-scale Rossby gyres are shifted more westward than their large-scale counterparts at the same depth as the waves propagate downward from the level of excitation. Consequently, smaller-scale Rossby gyres appear to be separated from the larger-scale counterparts over a certain range of depth, forming the "double-gyre" pattern as illustrated in Figure 12.

On the other hand, in the strong damping scenario shown in Figure 13, baroclinic Kelvin modes dominates over baroclinic Rossby modes. Hence we consider predominant Kelvin waves only in the strong damping scenario. By the same explanation for the modest damping case but with the opposite tilt for the wave front of Kelvin waves, small-scale Kelvin waves are shifted more eastward than their large-scale counterparts at the same depth as the waves propagate downward from the level of excitation. As a result, smaller-scale Kelvin waves appear to emerge from the larger-scale counterparts and shift more eastward at greater depths, as illustrated in Figure 13. Note that the Kelvin modes do not possess any meridional velocities, and therefore the meridional velocity fields in the figure arise from the less-dominant Rossby modes in the strong damping scenario.

To be more specific, we further consider predominant vertical modes of the solution and inspect their amplitudes, defined by Re[∑_na_{m, n, l}Z_m(z)], along the depth in the vertical domain. The left (right) panel of Figure 14 shows the predominant vertical modes with m = 0.26 and m = 0.79 related to the Rossby (Kelvin) component in the modest (strong) damping scenario.⁸ These two vertical modes correspond to the equivalent depths about 150 km and 50 km respectively, selected from the power spectrum in Figure 2, referred to as the larger-scale mode (solid) and smaller-scale mode (dashed). In fact, the large excitation of these two modes results not only from a good match with the vertical heating profile but also from a strong coupling with the longitudinal heating profile. The dotted line is the sum of the two modes and so the constructive or destructive phase interference can be seen. The two modes have the same phase and thus the maximum constructive effect around the heating level. This implies that they are simultaneously excited at the heating level. However, they have opposite signs and equal amplitude near 1 bar, implying that after traveling over a certain vertical distance the dispersive equatorial waves spatially split, leading to the dispersive wave patterns along the latitude.

Zoom In Zoom Out Reset image size

Before closing the section, we would like to reiterate that as explained in Section 2.1, the wave dispersion for Rossby waves can arise from baroclinic modes (denoted by m) in addition to horizontal modes (denoted by n, see Figure 5) resulting from the shallow-water model. Generally, we find that modes of higher order n are one order of magnitude smaller than the n = 1 mode in amplitude, and therefore do not contribute to the double gyres. It is the baroclinic modes that are responsible for double gyres, a dispersive phenomenon unique to vertical stratification.

We have studied the 3D structure and physical properties of equatorial waves in the presence of a zonal-mean zonal flow. We now proceed to the feedback of these waves to mean zonal flows. We first compare the distribution of wave momentum fluxes from the linear theory with the results from simulations. We then apply the linear model to explore the possible evolution of super-rotating flows due to the feedback, i.e., the interaction between waves and mean zonal flows.

3.3.1. Spatial Distributions of Momentum Fluxes and Mean-flow Accelerations

Figure 15 shows the spatial distributions of horizontal and vertical momentum fluxes as well as their gradients from our numerical simulation with ν = 10⁹ cm² s⁻¹ in the final equilibrium state, which are compared with those from our linear model for α = 0.1 and u₀ = 1000 m s⁻¹ illustrated in Figure 16. Owing to thermal forcing and dissipation, the gradients of wave momentum fluxes give the force acting on different parts of the atmosphere, as elaborated in Section 2.5. The linear results for the momentum fluxes and their gradients are qualitatively consistent with the results from the simulation, except that in the simulations the distributions of the momentum fluxes and corresponding gradients lie deeper in the atmosphere than those obtained from the linear model, which might be due to different vertical density profiles for the background state of the atmosphere. We can see from the figures that the horizontal momentum flux, $\rho _{0} \overline{u^{\prime }v^{\prime }}$, is negative in the northern hemisphere and positive in the southern hemisphere, creating the equatorward momentum transport to accelerate the equatorial jet around the heating level. On the other hand, the vertical momentum flux $\rho _{0} \overline{u^{\prime }w^{\prime }}$ is negative around the heating level and slightly positive beneath it in the simulation, thereby decelerating the super-rotating flow around the heating level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The equatorward momentum flux has been discussed in Section 3.2.1. Here, we focus on the mechanism for the vertical momentum transport by studying velocity fields in our linear model. As has been described in Section 3.2.2. for u' and w' at the equator, upwelling (downwelling) flows primarily occur in the horizontal westward (eastward) flows in Figure 8. Therefore in the vertical domain of the moderate damping case, the horizontal and vertical wave velocities have negative correlation $\overline{u^{\prime }w^{\prime }}<$ 0, implying an upward transport of westward momentum and a downward transport of eastward momentum, leading to the eastward acceleration at the bottom of heating region and westward acceleration (i.e., deceleration) above and around the heating level, in agreement with the result in the right panels of Figure 16. Moreover, the mechanism of jet acceleration and deceleration in our 3D linear analysis generally agrees with the 3D numerical results performed by SP as well. The vertical (horizontal) momentum transport decelerates (accelerates) the super-rotating jet.

The net force per unit volume acting on the super-rotating jet as well as zonal flows away from the equator are further illustrated in Figure 17, which plots the divergence of the total wave momentum flux for different u₀ in the linear theory and for the final equilibrium state in the simulation (see Equation (18)). The forces are strongest due to the Rossby resonance when u₀ = 500 m s⁻¹ and subside when u₀ = 1000 m s⁻¹. The results hint that the zonal-mean flow would be accelerated from rest, eastward around the equatorial region but westward at mid-latitudes and higher altitudes. The result for u₀ = 1000 m s⁻¹ from the linear theory (bottom left panel) resembles that from the simulation (bottom right panel), strongly suggesting that equatorial waves play a significant role in the formation and evolution of zonal flows observed in the simulation. The implication of our linear wave theory to possible evolutions of the equatorial super-rotation flow will be explored in more detail in the next section.

Zoom In Zoom Out Reset image size

3.3.2. Implication to Evolution of a Superrotating Flow

As has been suggested by the 3D simulation of SP, the evolution of the zonal flow mainly depends on the interaction with waves, and so the final equilibrium state occurs when the averaged accelerations from the waves vanishes. In our 3D linear analysis, we do not evolve linearized wave equations but solve for the solutions in a "quasi-steady state" for any given u₀ as a free parameter, and then calculate the resulting acceleration. Here we relax the "steady state" into "quasi-steady state," to describe the possible evolution of u₀ by a sequence of intermediate quasi-steady states. As long as the acceleration is not too abrupt, u₀ may vary continuously from one quasi-steady state to the next. Figure 18 shows the accelerations of the equatorial jet as a function of u₀ for the three damping scenarios derived from our linear model. The accelerations in the figure are taken averaged across the vertical domain roughly from 20 to 3 × 10⁻⁴ bars, and from 20° south of the equator to 20° north of the equator, representing the latitudinal region of the equatorial jet as suggested by Figure 17. The red (blue) curve indicates the acceleration contributed from the horizontal (vertical) wave transport. The net acceleration is depicted by the black curve.

Zoom In Zoom Out Reset image size

Figure 18 can approximately describe the evolution of the zonal flow if each steady state for a given u₀ behaves as an intermediate steady state such that the entire evolution of u₀ can be considered as a sequence of these intermediate steady states. There is no straightforward way to justify whether this approximation is reasonable unless a time evolution is performed to compare. Indeed, the approximation should be invalid if the evolution timescale of u₀ in the steady state is much shorter than the relaxation timescale for the entire atmosphere reaching an intermediate steady state. Hereafter, we just boldly assume that the approximation is applicable in order to explore the possible evolution of a super-rotating flow when the relaxation time is shorter than the evolution time of u₀, estimated based on the steady-state results. As explained in Section 2.5, the evolution of u₀ resulting from the Reynolds stress is a second-order process in the linear theory and is therefore slower than the relaxation of the entire atmosphere in response to the change of u₀. However as we shall see, the linear theory breaks down when the dissipation is not strong enough to suppress the wave resonance.

The relaxation timescale is given by the longer of the dynamical and thermal timescales when planetary waves are not severely damped.⁹ We still apply the same estimate in our strong damping scenario with α = 0.5, in the sense that the damping rate is not much larger than the dynamical time, which will be defined below. Because in our model the Newtonian cooling rate is set to equal the Rayleigh friction rate, the thermal timescale is just the damping time 1/α. The dynamical timescale of the atmosphere is characterized by the wave travel time across the equatorial Rossby deformation radius ¹⁰:

$$\begin{equation} T \sim \frac{L}{\sqrt{gH}} \sim 10^{4} \ {\rm s}, \end{equation} \tag{ 25 }$$

where L is the equatorial Rossby radius, presenting the horizontal scale of planetary waves (for HD 189733b, $L = (\sqrt{gH}/\beta)^{1/2} \approx 5 \times 10^{7}$ m), and $\sqrt{gH}$ is the speed of the shallow-water gravity wave. Since α is normalized roughly by 1/T, the thermal timescales corresponding to α = 0.5, 0.1, and 0.01 are about 2T, 10T, and 100T, respectively. Obviously, the relaxation time is mainly governed by the thermal timescale in the damping scenarios that we consider. On the other hand, the evolutionary timescale of the zonal flow τ_jet may be estimated from Figure 18 by taking u₀ divided by the corresponding acceleration. On average, for an eastward traveling jet (i.e., u₀ > 0), τ_jet is ≳ 10⁶ s for α = 0.5, ≳ 10⁵ s for α = 0.1, and ≳ 10⁴ s for α = 0.01. As a result, in the strong and moderate damping scenarios, the waves, although damped non-negligibly, are able to propagate and interact promptly to reach an intermediate steady state before u₀ has changed so significantly that the new u₀ would create a different wave state. In the small damping scenario, τ_jet can be comparable to the relaxation time because the Rossby resonance is not severely suppressed, even leading to significant westward accelerations in the regime of u₀ = 200–400 m s⁻¹. Consequently, our linear model breaks down in the small damping scenario. We henceforth focus only on the evolutions in the strong and moderate damping scenarios.

In the moderate and strong damping scenarios, we can see that the net acceleration (black curves) due to the wave momentum flux continuously accelerates the mean flow from an initial rest state, and gradually declines to zero after u₀ exceeds 1000 m s⁻¹. The bump of the net acceleration in the moderate damping scenario results from the resonance. The time evolution process generally resembles that shown in Figure 11 of SP, despite different implementations of the damping processes. Note that the accelerations via horizontal and vertical transports in our model initially increase because of the Rossby resonance and then subside, which is also suggested by Arnold et al. (2012). It would be interesting to know whether the decline in the acceleration curves in SP starting from u₀ ∼ 500 m s⁻¹ is attributable to the Rossby resonance. Also note that we ignore the drag acting on the zonal-mean flow, which would increase the rate that the acceleration curve approaches zero, and so a more realistic value of the final equilibrium velocity u₀ is expected to have a smaller value. Since the momentum flux is the product of wave velocities, the wave flux and the corresponding zonal acceleration scale quadratically with heating amplitude. Thus, our linear wave mechanism suggests that there is no threshold in the heating magnitude for generating super-rotation (similar to SP). However, this is because the drag acting on the zonal-mean flow is ignored here, which may be responsible for inhibiting a super-rotating jet when heating is too weak.

The intersections of the the net acceleration curve with the x-axis in Figure 18 indicate the equilibrium states. The negative slope of the acceleration curve around an intersection point implies a stable equilibrium such as u₀ around 3000 m s⁻¹ for α = 0.1 and around 2800 m s⁻¹ for α = 0.5. We may conclude from our linear model that when the damping rate is not too small, there is one unstable equilibrium state for negative u₀, and only one stable equilibrium state for positive u₀ over the tested range. The final equilibrium state u₀ should lie at the positive stable equilibrium point, or linearly runaway when the westward u₀ exceeds the unstable point, which occurs only when the initial state provides a very strong westward zonal-mean flow (u₀ = −1000 m s⁻¹ for α = 0.5 and u₀ = −2000 m s⁻¹ for α = 0.1). That is to say, the speed of the zonal-mean flow (or the total zonal-mean kinetic energy) in the final equilibrium state may have little dependence on the initial conditions according to our linear wave model with the strong and modest damping.