Temperature Structures Associated with Different Components of the Atmospheric Circulation on Tidally Locked Exoplanets

The most observationally accessible exoplanets tend to be on close-in, short-period orbits. Many are sufficiently close to their host star that the gravitational interaction between star and planet synchronizes their orbital and rotational periods (Dole 1964; Guillot et al. 1996). A planet in this orbital state is referred to as tidally locked. As tidally locked planets always present the same face to their host star, they have a permanent dayside and nightside.

These planets can be characterized by observing their thermal emission (Burrows 2014; Crossfield 2015). Time-resolved observations of thermal emission over an orbital phase are known as thermal phase curves, and show the disk-integrated infrared flux emitted from a planet as it orbits its host star (Cowan & Agol 2008; Demory et al. 2016), revealing the longitudinal structure of its brightness temperature. Eclipse maps show the two-dimensional structure of brightness temperature on a planet's dayside (Majeau et al. 2012).

Atmospheric temperatures and atmospheric circulation are closely related. This is not only due to the ability of an atmosphere to transport heat (Koll & Abbot 2016), and thus maintain a temperature structure against differential heating and cooling (Showman et al. 2013), but also through the pressure gradient forces in the horizontal momentum equations. These pressure gradient forces mean that any distribution of winds must be balanced by a particular temperature distribution. The most familiar example of such a balance is the geostrophic thermal wind relation (Holton 2004, Chapter 3)

$$\begin{eqnarray}&&f\displaystyle \frac{\partial u}{\partial \mathrm{lnp}}=\displaystyle \frac{R}{a}\displaystyle \frac{\partial T}{\partial \vartheta };\quad f\displaystyle \frac{\partial v}{\partial \mathrm{lnp}}=-\displaystyle \frac{R}{a\cos \vartheta }\displaystyle \frac{\partial T}{\partial \lambda },\end{eqnarray} \tag{ 1 }$$

which can be derived by assuming the flow is in hydrostatic and geostrophic balance. Above, p is pressure, a is the planetary radius, R is the specific gas constant, and $f=2{\rm{\Omega }}\sin \vartheta$ is the Coriolis parameter, with Ω the planetary rotation rate. λ is longitude and ϑ is latitude. The thermal wind relation states that vertical variation in the horizontal wind u = (u, v) must be accompanied by horizontal variation in the temperature T. The utility of geostrophic balance itself is limited to scenarios where the Coriolis term is dominant in the momentum equations, but the requirement that accelerations in the horizontal momentum equations are balanced by horizontal pressure contrasts is a generic one. We will show how an understanding of how different circulation features are balanced by atmospheric temperature structure aids interpretation of observations of thermal emission from extrasolar planets.

Tidally locked planets are heated on their dayside and cool from their nightside, and this pattern of heating and cooling strongly influences their atmospheric circulation (see reviews by Showman et al. 2013; Pierrehumbert & Hammond 2019; Zhang 2020). Heating on the dayside drives divergent motion (Showman et al. 2013), which in turn generates stationary Rossby waves (Showman & Polvani 2010; see also Sardeshmukh & Hoskins 1988). Interaction between these two circulation components can lead to momentum convergence toward the equator, and the acceleration of a westerly superrotating jet (Showman & Polvani 2011).

In Hammond & Lewis (2021), we showed that each of these circulation components can be isolated from one another by separating the horizontal velocity u into a rotational (divergence-free) component u _r and a divergent (vorticity-free) component u _d (called a Helmholtz decomposition):

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{u}} & = & {{\boldsymbol{u}}}_{{\rm{r}}}+{{\boldsymbol{u}}}_{{\rm{d}}},\\ & = & {\boldsymbol{k}}\times {{\rm{\nabla }}}_{p}{\psi }+{{\rm{\nabla }}}_{p}{\chi }.\end{array}\end{eqnarray} \tag{ 2 }$$

The zonal-mean zonal jet and stationary Rossby waves are contained within the zonal-mean and eddy components of ${{\boldsymbol{u}}}_{{\rm{r}}}={\overline{{\boldsymbol{u}}}}_{{\rm{r}}}+{{\boldsymbol{u}}}_{{\rm{r}}}^{{\prime} }$, respectively, and the divergent circulation (comprising, e.g., thermally-direct overturning; see Hammond & Lewis 2021, and Kelvin-like waves; see Appendix A) is contained within u _d. Above, χ is the velocity potential and ψ is a streamfunction, defined in terms of the divergence δ = ∇_p · u and vorticity ζ = k · ∇_p × u via the relations ${{\rm{\nabla }}}_{p}^{2}\chi =\delta$ and ${{\rm{\nabla }}}_{p}^{2}\psi =\zeta$. ∇_p is the horizontal gradient operator acting on isobaric surfaces, and k is the unit vector in the vertical direction.

In this paper, our objective is to use the Helmholtz decomposition to understand how different components of atmospheric circulation contribute to the distribution of atmospheric temperature, and by extension thermal emission, on tidally locked planets. This is achieved by decomposing the geopotential ϕ = gz into a component that is balanced by purely rotational winds (the rotational geopotential), and a component that is due to divergent winds and rotational-divergent interactions (the divergent geopotential; Section 3; see also Trenberth & Chen 1988). We then use the hydrostatic relation

$$\begin{eqnarray}&&\displaystyle \frac{\partial \phi }{\partial \,\mathrm{ln}\,{\rm{p}}}=-{RT}\end{eqnarray} \tag{ 3 }$$

to split the global temperature field into rotational and divergent components (Section 4). Finally, we show how each circulation component contributes separately to simulated thermal phase curves by linearly expanding I = σ T⁴ around the horizontally averaged temperature T₀ (Section 5). A summary of our results is included at the end of the paper (Section 6).

The diagnostics introduced in Sections 3–5 are applied to output from general circulation model (GCM) simulations of dry, temperate terrestrial tidally locked planets. The GCM configuration we use is described in Section 2. The atmospheric circulation on tidally locked planets depends strongly on planetary rotation rate, as has been demonstrated by a number of numerical modeling studies (Merlis & Schneider 2010; Carone et al. 2015; Noda et al. 2017; Haqq-Misra et al. 2018; Hammond & Pierrehumbert 2018). To illustrate how partitioning the geopotential and temperature into rotational and divergent components depends on planetary rotation, we present results from four simulations with different rotation periods, P = 4, 8, 16, and ∞ days (i.e., zero rotation). The zero-rotator is included to illustrate the effect of rotation in the other simulations. We also analyzed output from a simulation with P = 32 days. The results from this analysis were nearly identical to those for the P = 16 day simulation and so are not shown. We note that excluding moisture from our simulations will affect the simulated circulation, however we stress that the method we present to partition the geopotential and temperature does not depend on this assumption. To keep the number of simulations we analyze manageable, we only present results from simulations of terrestrial planets. However, there is no reason why the techniques described herein cannot be applied to hot Jupiters, and this is something we aim to do in the near future.

The rotational and divergent geopotential components are defined by assuming different balances of terms in the divergence equation (Equation (17); introduced in Section 3). Appendix A illustrates the relationship between balances of terms in the divergence equations, and balances of terms in decomposed momentum equations for the rotational and divergent winds, using the simple example of the linearized shallow water equations. By deriving momentum equations for the rotational and divergent velocities, we show that a steady linear Kelvin-like wave can exist on a sphere in the absence of drag, in contrast to the equatorial beta plane where drag is a condition for a linear Kelvin wave (Showman & Polvani 2011).

The GCM simulations are run using Isca, which is a framework for building idealized GCMs of varying complexity (Vallis et al. 2018). Isca is built on the GFDL ¹ spectral dynamical core, which is used to integrate the primitive equations of motion forward in time. In the present study we use a dry GCM configuration, i.e., water vapor and moist processes are not included in the model atmosphere. The model is configured to simulate the atmospheric circulation of terrestrial planets with a lower boundary modeled as a slab ocean. The main components of the model are described below.

All planetary parameters aside from the rotation rate and properties of the stellar irradiation (e.g., radius, surface pressure, gravitational acceleration) are set to be those of the Earth. As stated in the introduction, we use results from four simulations with different rotation periods, P = 4, 8, 16, and ∞ days. To facilitate comparison between our simulations and those of previous work (e.g., Koll & Abbot 2016; Haqq-Misra et al. 2018), we quote the external nondimensional equatorial deformation radii ² for each of our simulations, which are ${{ \mathcal L }}_{\mathrm{eq}}/a=1.08$, 1.53, 2.16, and ∞ for the P = 4, 8, 16, and ∞ day simulations, respectively.

2.1. Dynamical Core

2.2. Radiative Transfer

We use a semi-gray radiative transfer code to model radiative heating and cooling. Radiative fluxes are split into two wavelength bands: one covering longwave (thermal emission) wavelengths, and the other covering shortwave (stellar heating) wavelengths.

In the shortwave band, the atmosphere is assumed to be transparent, and all incoming radiation is absorbed at the surface (i.e., the surface albedo has been folded into the solar constant). The instellation at the top-of-atmosphere is imposed with the following distribution, appropriate for a tidally locked planet:

$$\begin{eqnarray}S=\left\{\begin{array}{ll}{S}_{0}\cos \vartheta \cos \left(\lambda -{\lambda }_{0}\right) & \mathrm{on}\ \mathrm{the}\ \mathrm{dayside},\\ 0 & \mathrm{on}\ \mathrm{the}\ \mathrm{nightside},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where λ₀ = 0° is the substellar longitude, and S₀ =1000 W m⁻².

In the longwave band, upward and downward radiative fluxes are computed according to

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{\rm{I}}}^{\uparrow }}{{\rm{d}}{\tau }}={I}^{\uparrow }-{\sigma }{T}^{4},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{\rm{I}}}^{\downarrow }}{{\rm{d}}{\tau }}={\sigma }{T}^{4}-{I}^{\downarrow },\end{eqnarray} \tag{ 6 }$$

where τ = τ₀(p/p₀) is the longwave optical depth, and we set τ₀ = 1 and p₀ = p_s. The longwave radiative heating in the model is then

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}\,=\,\cdots +\displaystyle \frac{g}{{c}_{p}}\displaystyle \frac{\partial \left({I}^{\uparrow }-{I}^{\downarrow }\right)}{\partial p}.\end{eqnarray} \tag{ 7 }$$

The boundary conditions for the longwave fluxes are I^↓(p = 0) = 0 and ${I}^{\uparrow }(p={p}_{{\rm{s}}})=\sigma {T}_{{\rm{s}}}^{4}$, where T_s is the surface temperature.

2.3. Surface Energy Budget, Boundary Layer Sensible Heat Transport, and Convection

The surface is modeled as a static slab ocean that can exchange energy with the atmosphere in the vertical but cannot transport heat horizontally. The surface temperature evolves according to

$$\begin{eqnarray}&&C\displaystyle \frac{\partial {T}_{{\rm{s}}}}{\partial t}=S+{I}_{\mathrm{sfc}}^{\downarrow }-\sigma {T}_{{\rm{s}}}^{4}-{ \mathcal H },\end{eqnarray} \tag{ 8 }$$

where C is a specified heat capacity that corresponds to an ocean mixed-layer depth of 2.5 m, ${I}_{\mathrm{sfc}}^{\downarrow }$ is the downward flux of longwave radiation incident on the surface, and ${ \mathcal H }$ is the surface sensible heat flux. ${ \mathcal H }$ is computed according to

$$\begin{eqnarray}&&{ \mathcal H }={\rho }_{{\rm{a}}}{c}_{p}{ \mathcal C }\left|{{\boldsymbol{u}}}_{{\rm{a}}}\right|\left({T}_{{\rm{s}}}-{T}_{{\rm{a}}}\right),\end{eqnarray} \tag{ 9 }$$

where ρ is the density, the subscript "a" indicates quantities that are evaluated at the lowest model level, and ${ \mathcal C }=0.001$ is a dimensionless bulk transfer coefficient.

Turbulent vertical transport of heat within the boundary layer is parametrized following Thatcher & Jablonowski (2016) as

$$\begin{eqnarray}&&{{ \mathcal D }}_{T}=\overline{w^{\prime} \theta ^{\prime} }=-{ \mathcal K }\displaystyle \frac{\partial \theta }{\partial z},\end{eqnarray} \tag{ 10 }$$

where z is height, w = D z /D t is the vertical velocity, and θ is the potential temperature. The diffusion coefficient ${ \mathcal K }$ is calculated according to

$$\begin{eqnarray}\begin{array}{c}{ \mathcal K }=\left\{\begin{array}{ll}{ \mathcal C }\left|{{\boldsymbol{u}}}_{{\rm{a}}}\right|{z}_{{\rm{a}}} & {\rm{f}}{\rm{o}}{\rm{r}}\,p\gt {p}_{{\rm{p}}{\rm{b}}{\rm{l}}},\\ { \mathcal C }\left|{{\boldsymbol{u}}}_{{\rm{a}}}\right|{z}_{{\rm{a}}}\exp \left[-{\left(\displaystyle \frac{{p}_{{\rm{p}}{\rm{b}}{\rm{l}}}-p}{{p}_{{\rm{s}}{\rm{t}}{\rm{r}}{\rm{a}}{\rm{t}}}}\right)}^{2}\right] & {\rm{f}}{\rm{o}}{\rm{r}}\,p\leqslant {p}_{{\rm{p}}{\rm{b}}{\rm{l}}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 11 }$$

p_pbl = 850 hPa is the top of the boundary layer, and p_strat = 100 hPa controls the rate of decrease of boundary layer diffusion with height. The contribution to the model temperature tendency from the boundary layer diffusion is then

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}\,=\cdots +\,\displaystyle \frac{g}{{c}_{p}}\displaystyle \frac{\partial {{ \mathcal D }}_{T}}{\partial p}.\end{eqnarray} \tag{ 12 }$$

The model includes a dry convection scheme, which instantaneously restores the temperature structure to the dry adiabat whenever it is convectively unstable.

2.4. Surface Friction and Sponge Layer

Surface friction is parametrized following Held & Suarez (1994) and Thatcher & Jablonowski (2016) as a linear Rayleigh drag

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}t}\,=\,\cdots -{k}_{{\rm{f}}}{\boldsymbol{u}},\end{eqnarray} \tag{ 13 }$$

where k_f is defined as

$$\begin{eqnarray}&&{k}_{{\rm{f}}}={k}_{{\rm{f}},{\rm{s}}}max\left(0,\displaystyle \frac{{\sigma }-{{\sigma }}_{{\rm{b}}}}{1-{{\sigma }}_{{\rm{b}}}}\right).\end{eqnarray} \tag{ 14 }$$

k_f,s = 1 day⁻¹ is the drag timescale at the surface, and the top of the frictional boundary layer is located at σ_b = 0.7. We use a linear drag in our simulations as it simplifies the analytic partitioning of the geopotential into components associated with rotational and divergent circulation (see the simplicity of the "k_f δ" term in Equation (18), Section 3). The circulation in our experiments is similar to that obtained from models that use more sophisticated boundary layer parameterizations (e.g., Koll & Abbot 2016; Noda et al. 2017), and so it is unlikely that using the Rayleigh drag strongly affects our results.

A linear drag is also applied in the upper atmosphere above p_sponge = 50 hPa following Polvani & Kushner (2002), which damps the horizontal winds on a timescale of 0.5 day⁻¹. This acts as a sponge layer, and is included to suppress wave reflection at the model top.

3.1. Definition of Rotational Geopotential

The primitive horizontal momentum equations can be written in the form (Holton 2004, Chapter 4)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}t}+{\omega }\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}p}+\left(f+{\zeta }\right){\boldsymbol{k}}\times {\boldsymbol{u}}=-{{\rm{\nabla }}}_{p}\left({\phi }+\displaystyle \frac{1}{2}{\left|{\boldsymbol{u}}\right|}^{2}\right)-{k}_{{\rm{f}}}{\boldsymbol{u}},\end{eqnarray} \tag{ 15 }$$

where ${\left|{\boldsymbol{u}}\right|}^{2}={u}^{2}+{v}^{2}$, ω = Dp/Dt is the pressure vertical velocity, and $f=2{\rm{\Omega }}\sin \vartheta$ is the Coriolis parameter. k_f is a linear drag coefficient that is nonzero close to the surface. Equations for the vorticity ζ and divergence δ can then be obtained by taking k · ∇_p × (15) and ∇_p · (15), respectively,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{D}}\left(f+{\zeta }\right)}{{\rm{D}}t}+\left(f+{\zeta }\right){\delta }+{\boldsymbol{k}}\cdot \left(\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}p}\times {{\rm{\nabla }}}_{p}{\omega }\right)=-{k}_{{\rm{f}}}{\zeta },\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{{\rm{\partial }}{\delta }}{{\rm{\partial }}t}+{{\rm{\nabla }}}_{p}\cdot \left({\omega }\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}p}\right)+{{\rm{\nabla }}}_{p}\cdot \left[\left(f+{\zeta }\right){\boldsymbol{k}}\times {\boldsymbol{u}}\right]\\ \,=\,-{{\rm{\nabla }}}_{p}^{2}\left({\phi }+\displaystyle \frac{{\left|{\boldsymbol{u}}\right|}^{2}}{2}\right)-{k}_{{\rm{f}}}{\delta }.\end{array}\end{eqnarray} \tag{ 17 }$$

The geopotential ϕ does not appear in the vorticity equation (Equation (16)), and thus all of the dynamics that determines ϕ is described by the divergence equation (Equation (17)).

In order to separate the geopotential into components associated with the rotational and divergent winds, we substitute u = u _r + u _d into Equation (17), which yields

$$\begin{eqnarray}\begin{array}{lll}\displaystyle \frac{\partial \delta }{\partial t}+{{\rm{\nabla }}}_{p}\cdot \left[\omega \displaystyle \frac{\partial \left({{\boldsymbol{u}}}_{{\rm{r}}}+{{\boldsymbol{u}}}_{{\rm{d}}}\right)}{\partial p}\right]+{u}_{{\rm{d}}}\beta +J\left(\chi ,{{\rm{\nabla }}}_{p}^{2}\psi \right) & \,-\,{{\rm{\nabla }}}_{p}\cdot \left[\left(f+{{\rm{\nabla }}}_{p}^{2}\psi \right){{\rm{\nabla }}}_{p}\psi \right] & \,=\,-{{\rm{\nabla }}}_{p}^{2}\left[\phi +\displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\psi \right|}^{2}}{2}+\displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\chi \right|}^{2}}{2}+J\left(\psi ,\chi \right)\right]-{k}_{{\rm{f}}}\delta \end{array}\end{eqnarray} \tag{ 18 }$$

where J(ψ, χ) = k · (∇_p ψ × ∇_p χ) is the Jacobian, and β = a⁻¹∂f/∂ϑ, with a the planetary radius. A derivation of Equation (18) is given in Appendix B.

To decompose the geopotential into components associated with the rotational and divergent circulation we expand it as ϕ = ϕ₀ + ϕ_r + ϕ_d, where ϕ₀(t, p) is the horizontal-mean geopotential, which does not contribute to Equation (18) as ${{\rm{\nabla }}}_{p}^{2}{\phi }_{0}=0$. We then define the rotational geopotential as that which is balanced by purely rotational circulation (with no vertical motion ³ ), whence we obtain

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{p}^{2}{\phi }_{{\rm{r}}}\equiv {{\rm{\nabla }}}_{p}\cdot \left[\left(f+{{\rm{\nabla }}}_{p}^{2}\psi \right){{\rm{\nabla }}}_{p}\psi \right]-{{\rm{\nabla }}}_{p}^{2}\left(\displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\psi \right|}^{2}}{2}\right).\end{eqnarray} \tag{ 19 }$$

This equation is called the nonlinear balance equation (Lorenz 1960) and takes the place of the divergence equation for purely non-divergent horizontal circulation (see, e.g., Holton 2004, Chapter 11). For stationary circularly symmetric flow, it is equivalent to the gradient wind approximation. If ψ (i.e., the rotational wind distribution) is known, then Equation (19) can be inverted to obtain ϕ_r (e.g., by expressing Equation (19) in terms of spherical harmonics).

We note that we have used the primitive equations as our starting point, which are simplified with respect to the full, deep-atmosphere Navier–Stokes equations. Previous work has shown that the approximations made to derive the primitive equations are appropriate for the atmospheres of terrestrial planets and hot Jupiters (Mayne et al. 2014), but may not be appropriate for sub-Neptune type planets (Mayne et al. 2019).

3.2. Divergent Geopotential

Given the rotational geopotential as defined by Equation (19), we define the divergent geopotential to be that obtained as a residual from

$$\begin{eqnarray}&&{\phi }_{{\rm{d}}}\equiv \phi -({\phi }_{0}+{\phi }_{{\rm{r}}}).\end{eqnarray} \tag{ 20 }$$

The remaining terms in Equation (18), which determine ϕ_d, involve terms that are purely divergent, and terms due to interaction between the rotational and divergent winds.

3.3. Note on the Relationship between Balances in the Divergence Equation and the Momentum Equations

3.4. Application to Simulations

In this section, we apply the geopotential decomposition defined by Equations (19) and (20) to the height fields obtained in our GCM simulations. ϕ_r and ϕ_d are computed from the time-resolved model output, but time-averaged diagnostics are presented for simplicity. Time averaging means the effect of transient eddies (potentially arising due to barotropic and/or baroclinic instabilities; Read et al. 2020) on the decomposed height fields is not shown in our analysis. However, we would like to stress that the nonlinear balance given by Equation (19) holds locally in time. While beyond the scope of the present work, it is possible to consider the role of transients using the framework described in the preceding subsections.

Figure 1 shows the horizontal wind u = (u, v) and geopotential height z in the mid-troposphere (p = 480 hPa) for each GCM simulation. The full circulation is shown in the top row, and subsequent rows show the rotational circulation—split into contributions from the zonal mean (jet) and eddies (stationary waves)—and the divergent circulation. Figure 2 shows the vertical structure of the full circulation, eddy-rotational circulation, and divergent circulation averaged between ±40° latitude. In the P = 4 day simulation, the circulation is mostly contained within the rotational component, whereas in the P = 8 and P = 16 day simulations it is more evenly divided into rotational and divergent contributions. Finally, the circulation in the P = ∞ day simulation is purely divergent. In this subsection, we describe the important features of each component of the circulation that can be identified in Figures 1 and 2. We refer to Appendix C throughout, which contains additional analysis of the contributions of individual terms in the expanded divergence equation (Equation (18)) to each of the rotational and divergent height fields.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The zonal-mean rotational height field ${\overline{\phi }}_{{\rm{r}}}$ in each of the P = 4, 8, and 16 day simulations is in approximate gradient wind balance with the zonal-mean rotational wind ${\overline{u}}_{{\rm{r}}}$:

$$\begin{eqnarray}&&f{\overline{u}}_{{\rm{r}}}+\displaystyle \frac{{\overline{u}}_{{\rm{r}}}^{2}\tan \vartheta }{a}\approx -\displaystyle \frac{1}{a}\displaystyle \frac{\partial {\overline{\phi }}_{{\rm{r}}}}{\partial \vartheta },\end{eqnarray} \tag{ 21 }$$

which gives rise to a local height maximum centered on the equator (i.e., $\partial {\overline{\phi }}_{{\rm{r}}}/\partial \vartheta {| }_{\vartheta =0}=0;$ ${\partial }^{2}{\overline{\phi }}_{{\rm{r}}}/\partial {\vartheta }^{2}{| }_{\vartheta =0}\lt 0$), associated with the superrotating equatorial jet (Hammond & Pierrehumbert 2018). In the P = 4 day simulation, there are additional retrograde jets poleward of ±50° in each hemisphere, which lead to a secondary maximum in the zonal-mean rotational height field at each pole.

The eddy-rotational circulation in each simulation has a horizontal wind and geopotential structure characteristic of a stationary equatorial Rossby wave. The longitude of the eddy-rotational height maximum is located eastward of the substellar point due to a Doppler shift by the zonal-mean zonal wind (Tsai et al. 2014). For the rapidly rotating P = 4 day case, the structure and phase offset of the Rossby waves is very similar to that obtained in a linearized shallow water model with a prescribed background jet (e.g., see Figure 10 in Hammond & Pierrehumbert 2018). As the rotation period is reduced, the Rossby lobes with a positive geopotential anomaly become elongated with respect to those with a negative geopotential anomaly. Similar phenomenology is also exhibited in the GCM simulations presented in Edson et al. (2011, see their Figure 5). In Appendix C we show that this is due to nonlinear effects by partitioning the rotational height into contributions from the linear and nonlinear terms in the nonlinear balance equation (Equation (19)). Note that the eddy-rotational velocity vectors are not always perfectly parallel to the geopotential contours (see, e.g., λ = 10° for the P = 16 day simulation), indicating slight departures of the eddy-rotational flow from geostrophic balance (this is possible as the nonlinear balance described by Equation (19) includes—but is not limited to—geostrophically balanced flow). In the vertical, the stationary Rossby waves have a wavenumber-1 modulation (Tsai et al. 2014), as can be identified in the second row of Figure 2 (most clearly evident in the P = 4 and 8 day simulations).

The divergent circulation generally displays the most complex structure of all of the circulation components in Figures 1 and 2. The exception to this is the P = ∞ day simulation, for which the circulation is purely divergent and takes the form of a single isotropic, thermally-direct, overturning cell that features air rising on the dayside, near the substellar point, and sinking on the nightside. The geopotential structure associated with the overturning circulation is predominantly associated with the ∇_p · (ω∂ u _d/∂p) term in Equation (18), which is dominant in regions of ascent and descent. In the boundary layer the surface friction term is important, and in the outflow regions aloft the horizontal nonlinear $-{{\rm{\nabla }}}_{p}^{2}[{\left|{{\rm{\nabla }}}_{p}\chi \right|}^{2}/2]$ term is important (see Appendix C).

The balance of terms in Equation (18) for the P = 16 day simulation is similar to that in the zero-rotation simulation. In particular, the ∇_p · (ω∂ u _d/∂p) term is still important in these simulations in the region of strong ascent near the substellar point. However, additional terms in Equation (18) that describe interactions between the rotational (mostly ${\overline{u}}_{{\rm{r}}}$) and divergent circulations are also important (Appendix C), and our interpretation is that in the P = 16 day case the overturning circulation is tipped over by the eastward superrotating jet.

In the P = 4 day simulation, the structure of the divergent height field is very different from that in the P = 16 or P = ∞ day simulations, and is harder to interpret. Both its horizontal and vertical structure are somewhat wave-like, although very different to the structure of a linear Kelvin wave on the equatorial beta plane (see, e.g., Matsuno 1966). In Appendix C we show that the total divergent height structure shown in Figures 1 and 2 for the P = 4 day simulation is mostly due to the linear u_d β term, which takes the form of a Kelvin-like wave in a linearized shallow water model (Appendix A), and terms describing interaction between rotational and divergent horizontal winds. Finally, in the P = 8 day simulation, the divergent height field has features that are characteristic of both the slowly rotating simulations, with thermally-direct overturning, and the P = 4 day simulation with nonlinear and wave-like contributions to the height, suggesting that it resides in an intermediate regime.

In this section, we use the hydrostatic relation to obtain components of the atmospheric temperature that are associated with the rotational and divergent geopotential components. This is achieved by requiring that each component of the geopotential (ϕ₀, ϕ_r, and ϕ_d) is in hydrostatic balance with the corresponding component of the temperature (T₀, T_r, and T_d), i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\phi }_{0}}{\partial \mathrm{lnp}}\equiv -{{RT}}_{0};\quad \displaystyle \frac{\partial {\phi }_{{\rm{r}}}}{\partial \mathrm{lnp}}\equiv -{{RT}}_{{\rm{r}}};\quad \displaystyle \frac{\partial {\phi }_{{\rm{d}}}}{\partial \mathrm{lnp}}\equiv -{{RT}}_{{\rm{d}}}.\end{eqnarray} \tag{ 22 }$$

Understanding how each component of the circulation contributes to atmospheric temperature structure is important for interpreting observations of thermal emission from tidally locked planets, such as eclipse maps (Majeau et al. 2012) and phase curves (Cowan & Agol 2008; Demory et al. 2016; Morris et al. 2022). Previous work has shown that hydrostatic equilibrium is a good approximation for the atmospheres of hot Jupiters (Mayne et al. 2014) and sub-Neptunes (Mayne et al. 2019), in addition to terrestrial atmospheres.

Figure 3 shows the temperature structure in the lower troposphere (p = 750 hPa) for each GCM simulation. As with the geopotential height shown in Figure 1, the temperature is decomposed into contributions from the zonal-mean and eddy-rotational circulations, and the divergent circulation. Figure 4 shows the vertical structure of the full circulation, eddy-rotational circulation, and divergent circulation, averaged meridionally over all latitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Generally, the vertical structure of the temperature components shown in Figure 4 is similar to the structure of the height components shown in Figure 2. The main difference between them is a phase offset in their vertical structure, which appears as the temperature is obtained from the vertical derivative of the geopotential. This motivates our presentation of the horizontal temperature slices shown in Figure 3 at a deeper pressure (p = 750 hPa) than used for the height in Figure 1 (p = 480 hPa). This choice reflects the vertical phase offset between the temperature and geopotential, and means that the horizontal slices show contributions to height and temperature from the same circulation features. Consequently, the horizontal structures of the height components shown in Figure 1 and the temperature components shown in Figure 3 are similar.

In each of the simulations with rotation, there is an eastward hotspot shift from the substellar point. In the P = 8 and 16 day simulations, this is almost exclusively due to the rotational standing Rossby waves, which are Doppler shifted eastward from the substellar point by the zonal-mean jet (Tsai et al. 2014; Hammond & Pierrehumbert 2018). The zonal-mean rotational temperature that balances the zonal-mean rotational wind cannot contribute a longitudinal shift to the temperature structure by definition, although it does contribute toward brightening the equator with respect to the poles. The divergent temperature structure in the P = 16 day simulation is similar to that of the P = ∞ day simulation, which only features overturning circulation. The overturning circulation is a direct response to heating at the substellar point, and it appears that this strong coupling to the surface prevents it from being Doppler shifted much by the jet.

By contrast, in the P = 4 day simulation, the divergent temperature structure has a greater eastward shift than the eddy-rotational circulation (which is now only shifted by a few degrees longitude from the substellar point), although its amplitude is comparatively lower. The eastward shift of the divergent temperature in this simulation is consistent with the wave-like nature of its divergent circulation and height field, discussed in the previous section. The temperature structure associated with the zonal-mean rotational circulation is also different from that in the P = 8 and 16 day simulations, due to the retrograde jets that exist in this simulation in midlatitudes, and now contributes a brightening to both the equator and the poles, relative to the midlatitudes.

Thermal phase curves show the infrared flux emitted from a planet as it orbits its star. As a planet progresses through its orbit we see emission from different angles, so these observations characterize the longitudinal distribution of temperature on extrasolar planets (Cowan & Agol 2008; Koll & Abbot 2015). Thermal phase curves are already available for hot Jupiters and high-temperature rocky planets (Demory et al. 2016; Kreidberg et al. 2019), and should increase in quality in the coming years (Bean et al. 2018) following the successful launch of the James Webb Space Telescope. With JWST it should also be possible to produce thermal phase curves for smaller and cooler terrestrial planets (Deming et al. 2009).

In this section, we compute thermal phase curves from the decomposed temperature structures to show how each component contributes to their observable features. We approximate that the broadband thermal emission I^↑ originates from a specific radiating level p_rad, so that

$$\begin{eqnarray}&&{I}^{\uparrow }\left(p={p}_{\mathrm{rad}}\right)=\sigma {T}^{4}\left(p={p}_{\mathrm{rad}}\right).\end{eqnarray} \tag{ 23 }$$

We do this to simulate observing at a wavelength at which the atmosphere is optically thick, instead of using the model's actual semi-gray outgoing longwave radiation (OLR) because in our simulations the OLR is dominated by the surface emission and so is not an instructive demonstration of the potential contribution to the phase curve from the atmosphere. This also corresponds to the case of gaseous planets where the thermal emission is entirely due to the atmosphere and its circulation.

As Equation (23) is nonlinear in T, we cannot easily separate out contributions to I^↑ from each of T_r and T_d. We therefore linearize Equation (23) about the horizontal-mean temperature T₀,

$$\begin{eqnarray}&&{I}^{\uparrow }\approx \sigma {T}_{0}^{4}+4\sigma {T}_{0}^{3}\left(T-{T}_{0}\right)=\sigma {T}_{0}^{4}+4\sigma {T}_{0}^{3}\left({T}_{{\rm{r}}}+{T}_{{\rm{d}}}\right),\end{eqnarray} \tag{ 24 }$$

which separates contributions to I^↑ from T_r and T_d. Normalized thermal phase curves are then computed according to

$$\begin{eqnarray}&&I\left(\xi \right)=\displaystyle \frac{1}{\pi \sigma {T}_{0}^{4}}{\int }_{-\xi -\pi /2}^{-\xi +\pi /2}{\int }_{-\pi /2}^{\pi /2}{I}^{\uparrow }\cos \left(\lambda +\xi \right){\cos }^{2}\vartheta \,{\rm{d}}\vartheta \,{\rm{d}}\lambda ,\end{eqnarray} \tag{ 25 }$$

where ξ is the phase angle of the planet in its orbit (Cowan & Agol 2008).

Figure 5 shows phase curves computed for each GCM simulation using Equation (25), for a radiating level located at p = 750 hPa, which corresponds to the level for which the horizontal temperature structure was shown in Figure 3. Solid black curves show the phase curve computed using the nonlinear expression given by Equation (23), and dashed gray curves show the phase curve computed using the linear approximation given by Equation (24). There is generally good agreement between the phase curves computed using the nonlinear and linear expressions, and the difference between the two is reduced with increasing rotation period as deviations from the horizontal-mean temperature become smaller. This shows that the decomposition of the phase curve into rotational and divergent components by the linear approximation is useful for understanding the isobaric phase curves computed using the nonlinear expression. We would like to note that the linear approximation may be less accurate for hotter planets (e.g., hot Jupiters and terrestrial lava planets), where T − T₀ is larger. Indeed, we have found this to be the case when conducting a preliminary analysis of GCM simulations of the lava planet K2-141b. However, in this scenario we found that σ T⁴ could be approximated accurately by including second order terms in T in the Taylor expansion, and that in this case the "T_r T_d" term was small, so that the rotational and divergent contributions to σ T⁴ could still be separated.

Zoom In Zoom Out Reset image size

The contributions from the rotational and divergent components of the temperature to the linear approximation are shown in Figure 5 as dashed red and blue lines, respectively. In each simulation with rotation, the amplitude of the rotational contribution to the phase curve is larger than that of the divergent contribution, and leads to a negative offset in the phase curve peak from ξ = 90° (corresponding to an eastward hotspot shift in the temperature maps shown in Figure 3). As discussed in Section 4, this is entirely due to the Doppler-shifted Rossby waves, as the temperature structure that balances the jet cannot itself contribute a direct hotspot shift by definition (although it is the velocity of the jet that Doppler shifts the waves).

The negative phase offset associated with the stationary Rossby waves increases monotonically with increasing period. However, this trend is not reflected in the full phase curves, due to contributions from the divergent circulation. In the P = 4 day simulation, the divergent circulation, which in this simulation resembles a nonlinear wave that is Doppler shifted eastward, also contributes an eastward phase shift, that is in fact larger than that associated with the rotational circulation. This means that the peak of the full phase curve has a greater negative phase offset from ξ = 90° than it would if only the rotational circulation (stationary Rossby waves) contributed. By contrast, in the P = 16 day simulation, the divergent circulation now takes the form of overturning circulation, and is restricted to be close to the substellar point, due to strong coupling between the overturning circulation and the instellation. This means that the divergent contribution to the phase curve is not offset from ξ = 90°, and the full phase curve has a lesser negative phase offset than it would if due to the stationary Rossby waves alone.

Our analysis of the vertical temperature structure in the previous section shows that the phase curves in Figure 5 will depend strongly on our choice of radiating level. To illustrate this we show an additional phase curve in Figure 6 for the P = 4 day simulation, this time for p_rad = 200 hPa. In this figure, the sign of the normalized flux for the full and rotational curves is essentially flipped with respect to that shown in Figure 5, and the peak of the full phase curve has a positive offset from ξ = 0°, which would correspond to a westward hotspot shift in temperature. In our simulation, this is due to the wavenumber-1 structure of the rotational stationary Rossby waves in the vertical (see Figure 4). The baroclinic structure of the Rossby waves means that they contribute an eastward temperature offset at high pressures (lower altitudes) and a westward temperature offset at lower pressures (higher altitudes). Note that the direction of each of these temperature offsets would be reversed in the absence of a background zonal jet phase-shifting the waves eastward (see Figure 9(a) of Hammond & Pierrehumbert 2018; see also Tsai et al. 2014). The Rossby waves' vertical structure means that they can contribute either a negative or a positive shift to the phase curve peak depending on the exact pressure level from which they contribute to thermal emission. This is a key result of understanding thermal phase shifts as due to stationary wave shifts rather than heat advection by a jet, and is a potential explanation for unexpected observations of westward phase shifts. For example, Dang et al. (2018) identify a westward hotspot shift in the atmosphere of the hot Jupiter Corot-2b. Whether or not the mechanism described above can be used to explain this observation depends on whether the photospheric level is coincident with the westward-shifted or eastward-shifted component of the Rossby wave's baroclinic temperature structure. This should be investigated in the future using GCM simulations of atmospheric circulation on hot Jupiters (as opposed to the terrestrial simulations used in this work).

Zoom In Zoom Out Reset image size

The phase curves presented in this section are computed using the assumption that broadband thermal emission originates from a specific radiating level p = p_rad (Equation (23)). However, it is important to note that this will not always be the case, especially for an atmosphere where condensible species and clouds are present (Yang et al. 2013; Roman et al. 2021). Insofar as the linearization of σ T⁴ given by Equation (24) is accurate, it may be possible to account for this effect by introducing this linear decomposition into the two-stream equations for the upwelling and downwelling longwave radiative fluxes, to inform a partitioning of the radiative fluxes into rotational and divergent contributions. This is beyond the scope of the present work, but we aim to pursue it in a future publication.

Hammond & Lewis (2021) showed that the Helmholtz decomposition u = u _r + u _d is a useful tool for separating out different components of the atmospheric circulation on tidally locked planets. The jet and stationary Rossby waves are contained within u _r and divergent circulation is contained within u _d, which takes the form of thermally-direct overturning for slowly rotating planets (e.g., our P = 8 and P = 16 day simulations; see also Hammond & Lewis 2021), but has a complex linear Kelvin-like/nonlinear wave structure for more rapidly rotating planets (see our P = 4 day GCM simulation, and the shallow water results in Appendix A).

The objective of this work was to build upon the work in Hammond & Lewis (2021) by relating the atmospheric temperature structure to each of the circulation components contained within u _r and u _d. This was achieved by defining rotational and divergent components of the geopotential using balance relations in the divergence equation, which were then used to define a temperature structure via the hydrostatic relation. To illustrate the utility of decomposing the temperature in this way, we applied it to output from idealized GCM simulations of atmospheric circulation on terrestrial tidally locked planets (run using Isca; Vallis et al. 2018), considering four rotation periods: P = 4, 8, 16, and ∞ days (i.e., zero rotation).

The temperature maps we obtain (Figure 3) show that both the rotational and divergent circulations make non-negligible contributions to the horizontal temperature structure. Understanding what determines the relative strength of these circulations, and how they depend on properties of a particular exoplanet, will be crucial for correctly interpreting observations of thermal emission such as phase curves (Cowan & Agol 2008) and eclipse maps (Majeau et al. 2012).

Our analysis shows that both circulation components can contribute to the hotspot shifts that have been inferred from observations of thermal emission (e.g., Demory et al. 2016; Figures 3 and 5). In our simulations, both the rotational circulation and the divergent circulation contribute an eastward hotspot shift if thermal emission is assumed to come from the lower troposphere (p = 750 hPa). For our P = 8 and P = 16 day simulations, the hotspot shift is almost entirely due to the rotational stationary Rossby waves (which are Doppler shifted eastward by the superrotating jet). In the P = 8 day case, this is because the amplitude of the rotational temperature component is much larger than the divergent temperature component. In the P = 16 day case, both components have a similar amplitude, but the divergent component does not contribute a hotspot shift as it is dominated by thermally-driven overturning, which means it is strongly coupled to the pattern of instellation that is centered on the substellar point. By contrast, in the P = 4 day simulation, both the rotational component and the divergent component (which is now more wave-like) contribute an eastward hotspot shift if emission is from the lower troposphere. However, we show that the rotational component can also contribute a westward hotspot shift if emission comes from higher up in the atmosphere (p = 200 hPa in our simulations; Figure 6), which is due to its wavenumber-1 structure in the vertical (Figure 4).

From this analysis, we suggest that the temperature structure of the atmosphere of a tidally locked planet can be divided into four physically meaningful components:

• 1.  
  The overturning divergent component, providing a temperature peak at the substellar point and uniform temperature elsewhere.
• 2.  
  A wave-like divergent component, providing a zonal-wavenumber-1 sinusoidal temperature modulation on the equator, which can be Doppler shifted by a zonal jet.
• 3.  
  The Rossby wave-like eddy-rotational component, providing a zonal-wavenumber-1 sinusoidal temperature modulation with maxima off the equator, which can be Doppler shifted by a zonal jet.
• 4.  
  The zonal-mean rotational component, i.e., the zonal jet, providing a zonally uniform temperature field with a meridional gradient determined by the gradient of the jet itself.

We conclude by suggesting some opportunities for future work motivated by our results. First, the four temperature components we highlight above could be represented with simple parameterizations and used to fit observed thermal phase curves or eclipse maps. Dobbs-Dixon & Blecic (2022) and Chubb & Min (2022) use parameterized temperature structures representing GCM simulations and idealized models to fit observations. We aim to use our four-component decomposition of the temperature in this way in the future. Our components are each associated with specific wind fields (Hammond & Lewis 2021), raising the possibility of deriving the two-dimensional wind field from a two-dimensional temperature observation.

However, we have only applied the height and temperature decomposition to output from GCM simulations of atmospheric circulation on terrestrial tidally locked planets. Therefore, an obvious and important follow-up to this study would be to apply the same decomposition to output from hot Jupiter GCMs. While the theoretical framework we employ is not specific to terrestrial planets, it is not obvious a priori that the decomposed temperature components in a hot Jupiter atmosphere would exhibit the same structure as identified for our terrestrial experiments. We note, however, that in Hammond & Lewis (2021), a similar partitioning of the horizontal circulation is obtained from a GCM simulation of a terrestrial planet and a simulation of a hot Jupiter. This indicates that a decomposition of the temperature may also be similar for the two types of planet.

To relate the decomposed temperature maps presented in Figure 3 to thermal emission, we made the simplifying assumption that the broadband outgoing longwave radiation can be approximated by thermal emission from a specific radiating level. However, in general we would not expect this to be an accurate assumption, especially for the case of atmospheres that host radiatively active atmospheric constituents that are spatially inhomogeneous (e.g., clouds and condensible substances; Yang et al. 2013; Roman et al. 2021). Therefore, developing a method to partition the vertically integrated upwelling longwave flux into contributions from rotational and divergent temperature structures is an important area for future work.

Finally, we note that the in this study we have only decomposed the circulation obtained from dry GCM experiments. The effects of latent heating due to condensation have a direct effect on the atmospheric circulation (Merlis & Schneider 2010; Leconte et al. 2013), which might then partition differently into rotational and divergent components. In this scenario, the height and temperature structures that balance each circulation component would also be different to those obtained from the dry simulations here. Thus, it would be interesting to repeat the analysis here for a similar model configuration that included moisture.

We hope that the techniques and analysis presented herein will be of use for interpreting observations of thermal emission on tidally locked planets.

The authors thank Daniel Koll, Peter Read, and an anonymous reviewer for providing constructive feedback that greatly improved the quality of this paper. N.T.L. was supported by Science and Technology Facilities Council grant ST/S505638/1. M.H. was supported by a Junior Research Fellowship at Christ Church, Oxford. The authors are grateful to Dr. Man-Suen Chan for IT support.

In this appendix, we use the example of the linearized shallow water equations on the sphere to show how partitioning the height field into rotational and divergent components can be used to obtain momentum equations for the rotational and divergent wind. We use these equations to interpret the existence of a steady linear Kelvin-like wave in the spherical geometry, in the absence of drag, which does not exist on the equatorial beta plane.

In the absence of friction, the linearized shallow water equations may be written as (Vallis 2017, Chapter 3)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{\boldsymbol{u}}}{{\rm{\partial }}t}+f{\boldsymbol{k}}\times {\boldsymbol{u}}=-g{\rm{\nabla }}{\eta },\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \eta }{\partial t}+H\delta =Q.\end{eqnarray} \tag{ A2 }$$

Above, η is the free surface displacement (or height), H is the mean fluid depth, and Q is a mass source or sink term.

As with the primitive equations (Section 3), equations for the vorticity and divergence can be obtained by taking k · ∇ × (A1) and ∇ · (A1), respectively,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \zeta }{\partial t}+f\delta +\left({v}_{{\rm{r}}}+{v}_{{\rm{d}}}\right)\beta =0,\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \delta }{\partial t}-{\rm{\nabla }}\cdot \left[f{\rm{\nabla }}\psi \right]+{u}_{{\rm{d}}}\beta =-g{{\rm{\nabla }}}^{2}\eta .\end{eqnarray} \tag{ A4 }$$

Assuming a steady state, we can partition the divergence equation into two equations that define the rotational and divergent height, η_r and η_d:

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left[f{\rm{\nabla }}\psi \right]\equiv g{{\rm{\nabla }}}^{2}{\eta }_{{\rm{r}}};\quad -{u}_{{\rm{d}}}\beta \equiv g{{\rm{\nabla }}}^{2}{\eta }_{{\rm{d}}},\end{eqnarray} \tag{ A5 }$$

which can then be integrated to yield momentum equations for the rotational and divergent circulations in a steady state

$$\begin{eqnarray}&&f{\boldsymbol{k}}\times {{\boldsymbol{u}}}_{{\rm{r}}}=-g\left({\rm{\nabla }}{{\eta }}_{{\rm{r}}}+{\boldsymbol{k}}\times {\rm{\nabla }}{\gamma }\right),\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&f{\boldsymbol{k}}\times {{\boldsymbol{u}}}_{{\rm{d}}}=-g\left({\rm{\nabla }}{{\eta }}_{{\rm{d}}}-{\boldsymbol{k}}\times {\rm{\nabla }}{\gamma }\right),\end{eqnarray} \tag{ A7 }$$

where γ is defined by

$$\begin{eqnarray}&&g{{\rm{\nabla }}}^{2}\gamma \equiv \beta {v}_{{\rm{r}}}=-\left(f\delta +\beta {v}_{{\rm{d}}}\right).\end{eqnarray} \tag{ A8 }$$

Equation (A8) can be obtained by taking the curl of either Equation (A6) or (A7), and the second equality comes from the vorticity equation in steady state.

Equations (A6) and (A7) each contain an additional term involving the potential function γ, when compared with the momentum equation for the full horizontal velocity u . These terms arise as an integration constant and are due to the variability of f with latitude, and their proper representation is contingent on working in a spherical geometry (Trenberth & Chen 1988).

Figure 7 shows horizontal velocity and height fields from a linear shallow water simulation using the GFDL spectral shallow water model. We ran the model at T85 resolution for 500 days with a timestep of 300 s, and plot output averaged over the final 100 days of the run. The mass source term is configured for a tidally locked planet as

$$\begin{eqnarray}&&Q=\displaystyle \frac{{\eta }_{\mathrm{eq}}-\eta }{{\tau }_{\mathrm{rad}}},\end{eqnarray} \tag{ A9 }$$

following Perez-Becker & Showman (2013), where η_eq is given by

$$\begin{eqnarray}{\eta }_{\mathrm{eq}}=\left\{\begin{array}{ll}H+{\rm{\Delta }}{\eta }_{\mathrm{eq}}\cos \lambda \cos \vartheta & \mathrm{on}\ \mathrm{the}\ \mathrm{dayside},\\ H & \mathrm{on}\ \mathrm{the}\ \mathrm{nightside}.\end{array}\right.\end{eqnarray} \tag{ A10 }$$

Δη_eq is the difference in the radiative equilibrium height between the substellar point and the terminator. In Equation (A9), τ_rad is the radiative relaxation timescale. The GFDL model solves the full nonlinear shallow water equations, and the solution was kept in a linear regime by using a small forcing amplitude, η_eq/H ≪ 1, following Perez-Becker & Showman (2013). An important feature of the model is that it does not include a linear drag (i.e., τ_drag = ∞ in Perez-Becker & Showman 2013). The planetary parameters we used for this simulation were a = 6.371 × 10⁶ m, P = 8 days, and g = 9.81 m s⁻². The equilibrium height was H = 10,000 m, the magnitude of the forcing was Δη_eq = 10 m, and the radiative timescale was τ_rad = 0.1 day. We chose these parameters for demonstrative purposes as they give rotational and divergent components of similar magnitude, but our conclusions are not sensitive to the chosen parameters as long as the forcing is linear. In Figure 7 the velocity and height is split into eddy-rotational and divergent components. The zonal-mean rotational component is not shown as it is small.

Zoom In Zoom Out Reset image size

The decomposed rotational and divergent circulations shown in Figure 7 have structures that resemble the linear solutions for equatorial Rossby and Kelvin waves, respectively, on the equatorial beta plane (Matsuno 1966), although they are modified by the spherical geometry. This partitioning of the Rossby waves into the rotational component and Kelvin-like/gravity waves into the divergent component is consistent with Vanneste & Vial (1994), who show that in the limit of small Lamb parameter, Λ ≡ 4Ω² a²/(gH) ≪ 1 (in our simulation Λ ≈ 0.1), Rossby waves are purely rotational and gravity waves are purely divergent. The Kelvin-like wave obtained here has nonzero meridional velocity, which derives from the spherical geometry (Yamamoto 2019). The Rossby component is very similar to that obtained by Showman & Polvani (2011) on the equatorial beta plane in the absence of friction (as is the case here; see their Figure 14). Notably, Showman & Polvani (2011, their Appendix C) show that in the absence of drag, an equatorial Kelvin wave cannot exist on the equatorial beta plane.

In a steady state, the zonal momentum equation is

$$\begin{eqnarray}&&{fv}=\displaystyle \frac{g}{a\cos \vartheta }\displaystyle \frac{\partial \eta }{\partial \lambda }.\end{eqnarray} \tag{ A11 }$$

As f = 0 at the equator (and also v = 0 due to the hemispheric symmetry of the problem), Equation (A11) requires that ∂η/∂λ = 0 at the equator. This is satisfied in our linear shallow water simulation (see the left-hand panel in Figure 7). However, Showman & Polvani (2011) additionally show that on the equatorial beta plane, equatorial Rossby wave solutions have no longitudinal η variation at the equator, which by Equation (A11), means that an equatorial Kelvin component cannot exist (as it would introduce a nonzero ∂η/∂λ).

By contrast, both the equatorial Rossby and Kelvin-like components in Figure 7 have nonzero height modulations that cancel one another so that the total ∂η/∂λ at the equator is still zero. This is due to the effect of the spherical geometry on the waves, which now are governed by a zonal momentum equation of the form (e.g., for u_r, from Equation (A6))

$$\begin{eqnarray}&&{fv}_{{\rm{r}}}=g\left(\displaystyle \frac{1}{a\cos {\vartheta }}\displaystyle \frac{{\rm{\partial }}{{\eta }}_{{\rm{r}}}}{{\rm{\partial }}{\lambda }}+\displaystyle \frac{1}{a}\displaystyle \frac{{\rm{\partial }}{\gamma }}{{\rm{\partial }}{\vartheta }}\right).\end{eqnarray} \tag{ A12 }$$

If ∂γ/∂θ is nonzero at the equator (which will be the case if v is antisymmetric about the equator; γ obtained in our simulation is shown in Figure 8 for reference), then neither will ${\left.\partial {\eta }_{{\rm{r}}}/\partial \lambda \right|}_{\vartheta =0}=-{\left.\partial {\eta }_{{\rm{d}}}/\partial \lambda \right|}_{\vartheta =0}$. As a consequence, on the sphere and in the limit of small Lamb parameter, the existence of a rotational stationary Rossby wave requires the existence of a divergent stationary Kelvin-like wave, and vice versa. Showman & Polvani (2011) introduced a linear drag to their beta plane model to generate the stationary Kelvin wave necessary for the acceleration of a zonal jet, but we have shown here that this linear drag is not necessary for the formation of a Kelvin-like wave on a sphere.

Zoom In Zoom Out Reset image size

In this appendix we describe how the divergence equation (Equation (17) in the main text)

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \delta }{\partial t}+{{\rm{\nabla }}}_{p}\cdot \left(\omega \displaystyle \frac{\partial {\boldsymbol{u}}}{\partial p}\right)+{{\rm{\nabla }}}_{p}\cdot \left[\left(f+\zeta \right){\boldsymbol{k}}\times {\boldsymbol{u}}\right]\\ \qquad =\,-{{\rm{\nabla }}}_{p}^{2}\left(\phi +\displaystyle \frac{{\left|{\boldsymbol{u}}\right|}^{2}}{2}\right)-{k}_{{\rm{f}}}\delta \end{array}\end{eqnarray} \tag{ B1 }$$

can be expanded into a form that separates contributions from the rotational and divergent winds (i.e., Equation (18) in the main text).

First, applying the Helmholtz decomposition to ∇_p · [(f + ζ) k × u ], we obtain four terms:

• 1.  
  ∇_p · [f k × u _r] = ∇_p · [f k × k × ∇_p ψ] = − ∇_p · [f∇_p ψ],
• 2.  
  ${{\rm{\nabla }}}_{p}\cdot [{\zeta }{\boldsymbol{k}}\times {{\boldsymbol{u}}}_{{\rm{r}}}]=-{{\rm{\nabla }}}_{p}\cdot [{{\rm{\nabla }}}_{p}^{2}{\psi }{{\rm{\nabla }}}_{p}{\psi }]$,
• 3.  
  ∇_p · [f k × u _d] = k × u _d · ∇_p f = u_d β,
• 4.  
  $\begin{array}{rcl}{{\rm{\nabla }}}_{p}\cdot [{\zeta }{\boldsymbol{k}}\times {{\boldsymbol{u}}}_{{\rm{d}}}]={\boldsymbol{k}}\times {{\rm{\nabla }}}_{p}{\chi }\cdot {{\rm{\nabla }}}_{p}{\zeta } & = & {\boldsymbol{k}}\cdot ({{\rm{\nabla }}}_{p}{\chi }\times {{\rm{\nabla }}}_{p}{\zeta })\\ & = & J({\chi },{\zeta })=J({\chi },{{\rm{\nabla }}}_{p}^{2}{\psi })\end{array}$,

where β = ∂f/∂ϑ, for the third and fourth terms we have used (f + ζ)∇_p · ( k × u _d) = − (f + ζ) k · (∇_p × u _d) = 0, and in the fourth term we have defined the Jacobian

$$\begin{eqnarray}&&J({\chi },{\zeta })\equiv {\boldsymbol{k}}\cdot \left({{\rm{\nabla }}}_{p}{\chi }\times {{\rm{\nabla }}}_{p}{\zeta }\right).\end{eqnarray} \tag{ B2 }$$

We can also expand the ${\left|{\boldsymbol{u}}\right|}^{2}$ term into contributions from u _r and u _d as follows:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\left|{\boldsymbol{u}}\right|}^{2}}{2} & = & \displaystyle \frac{{\left|{{\boldsymbol{u}}}_{{\rm{r}}}\right|}^{2}}{2}+\displaystyle \frac{{\left|{{\boldsymbol{u}}}_{{\rm{d}}}\right|}^{2}}{2}+{{\boldsymbol{u}}}_{{\rm{r}}}\cdot {{\boldsymbol{u}}}_{{\rm{d}}}\\ & = & \displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\psi \right|}^{2}}{2}+\displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\chi \right|}^{2}}{2}+\,J\left(\psi ,\chi \right).\end{array}\end{eqnarray} \tag{ B3 }$$

With these substitutions, the divergence Equation (17) becomes

$$\begin{eqnarray}\begin{array}{lll}\displaystyle \frac{\partial \delta }{\partial t}+{{\rm{\nabla }}}_{p}\cdot \left[\omega \displaystyle \frac{\partial \left({{\boldsymbol{u}}}_{{\rm{r}}}+{{\boldsymbol{u}}}_{{\rm{d}}}\right)}{\partial p}\right]+{u}_{{\rm{d}}}\beta +J\left(\chi ,{{\rm{\nabla }}}_{p}^{2}\psi \right) & -\,{{\rm{\nabla }}}_{p}\cdot \left[\left(f+{{\rm{\nabla }}}_{p}^{2}\psi \right){{\rm{\nabla }}}_{p}\psi \right] & \,=\,-{{\rm{\nabla }}}_{p}^{2}\left[\phi +\displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\psi \right|}^{2}}{2}+\displaystyle \frac{{\left|{{\rm{\nabla }}}_{p}\chi \right|}^{2}}{2}+J\left(\psi ,\chi \right)\right]-{k}_{{\rm{f}}}\delta ,\end{array}\end{eqnarray} \tag{ B4 }$$

which is Equation (18) in the main text.

This appendix includes some additional analysis of the terms in the divergence equation (Equation (18)) for our simulations, to identify those that contribute most strongly to the rotational and divergent height components. In particular, we analyze the divergent circulation by highlighting the term that corresponds to overturning circulation, and comparing it to the terms that correspond to wave-like divergent circulation and nonlinear rotational-divergent interactions

In the P = 4 and P = 8 day simulations, the eddy-rotational height (Figure 1) is dominated by the linear ∇_p · [f∇_p ψ] term (not shown). It has the structure of a stationary equatorial Rossby wave that is Doppler shifted eastward by the zonal-mean jet (Tsai et al. 2014; Hammond & Pierrehumbert 2018). In the P = 16 day simulation, the eddy-rotational height also resembles a stationary Rossby wave, but the lobes with a positive geopotential anomaly are noticeably elongated with respect to those with a negative geopotential anomaly (Figure 1). Figure 9 show the eddy-rotational height field for the P = 16 day simulation, divided into contributions from the linear term and the nonlinear ${{\rm{\nabla }}}_{p}\cdot [{{\rm{\nabla }}}_{p}^{2}{{\rm{\nabla }}}_{p}\psi ]-{{\rm{\nabla }}}_{p}^{2}[{\left|{{\rm{\nabla }}}_{p}\psi \right|}^{2}/2]$ terms. This subdivision of the height field shows that the elongation of the stationary Rossby is mostly due to nonlinear rotational-rotational interactions.

Zoom In Zoom Out Reset image size

The divergent height field has a complex structure in each of the simulations presented in Figures 1 and 2. In order to interpret it, we show the divergent height field decomposed into contributing terms in Figures 10 and 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The P = ∞ day simulation represents a limiting case where the divergent circulation takes the form of an isotropic, overturning cell (see Figures 1 and 2). In this scenario, Figure 11 shows that the divergent height structure is due to the ∇_p · (ω∂ u _d∂p) term in regions of ascent near the substellar point, and descent on the nightside. Meanwhile, the k_f δ (surface friction) contribution is important in the boundary layer, and nonlinear horizontal divergent–divergent interactions are important in the outflow area surrounding the region of ascent near the substellar point. The balance of terms in the P = 16 day simulation is similar to that in the P = ∞ day case, although there is a contribution to the divergent height field from the horizontal nonlinear rotational-divergent interaction terms, $-{{\rm{\nabla }}}_{p}^{2}[J(\psi ,\chi )]-J(\chi ,{{\rm{\nabla }}}_{p}^{2}\psi )$ (mostly through ${\overline{u}}_{{\rm{r}}}$). Our interpretation of the divergent height structure in the P = 16 day simulation is that it is consistent with a thermally-direct overturning circulation that is tipped over by the superrotating jet.

By contrast, the most important contributions to the divergent height in the P = 4 day simulation are from the linear u_d β term, and nonlinear horizontal rotational-divergent interactions (note that surface friction is still important in the boundary layer). Figure 10 shows the horizontal structure of the linear u_d β term and horizontal nonlinear rotational-divergent interaction terms for the P = 4 day simulation. The height field associated with the u_d β term is qualitatively similar to that of a linear Kelvin-like wave obtained in the shallow water simulation presented in Appendix A. In the vertical, the height associated with this term has a wavenumber-1 structure (Figure 11). However, unlike in the linear shallow water simulation, the divergent height field is modified by the contribution from the nonlinear rotational-divergent interaction term (largest in the upper troposphere), which yields the complex structure shown in Figure 1. The divergent height in the P = 8 day simulation has contributions from all of the terms discussed for both the P = 4 and P = 16 day simulations, and it appears that this simulation is an intermediate case between the Kelvin-like/nonlinear wave regime of the P = 4 day simulation, and the thermally-direct overturning regime of the P = 16 and P = ∞ day simulations.