EFFECTS OF INITIAL FLOW ON CLOSE-IN PLANET ATMOSPHERIC CIRCULATION

The global dynamics of a shallow, 3D atmospheric layer is governed by the primitive equations (e.g., Pedlosky 1987; Holton 2004). Here, by "shallow" we mean that the thickness of the atmosphere under consideration is small compared to the planetary radius R_p. In atmospheric studies, pressure p is commonly used as the vertical coordinate. In the pressure coordinate system, these equations read

$$\begin{eqnarray} \frac{{D}{\bf v}}{{D}t} + \left(\frac{u}{R_p}\tan \phi \right){\bf k}\times {\bf v}\ & = & {-}\bnabla _p \Phi - f{\bf k}\times {\bf v} + {\cal D}_{\bf v},\quad \ \ \end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray} \frac{{\partial }\Phi }{{\partial }p}\ & = &{-}\frac{1}{\rho }, \end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray} \frac{{\partial }\omega }{{\partial }p}\ & = &{-}\bnabla _p\cdot {\bf v}, \end{eqnarray} \tag{ 1c }$$

$$\begin{eqnarray} \frac{{D}T}{{D} t} -\frac{\omega }{\rho \, c_p}\ & = &\ \frac{\dot{q}_{\rm {net}}}{c_p} + {\cal D}_T\,, \end{eqnarray} \tag{ 1d }$$

where

$$\begin{eqnarray*} &&\frac{{D}}{{D}t}\ =\ \frac{{\partial }}{{\partial }t} + {\bf v}\cdot \bnabla _p + \omega \!\frac{{\partial }}{{\partial }p}\,. \end{eqnarray*}$$

In the above equations, v(x, t) = (u, v) is the eastward, northward velocity in a frame rotating with Ω, the planetary rotation rate; Φ = gz is the geopotential, where g is the gravitational acceleration and z is the distance above the planetary radius R_p; k is the unit vector in the local vertical direction; f = 2 Ωsin ϕ is the Coriolis parameter, the projection of the planetary vorticity vector 2Ω onto k, with ϕ the latitude; $\bnabla _p$ is the horizontal gradient on a constant p-surface; ω ≡ Dp/Dt is the vertical velocity; ρ is the density; ${\cal D}_{\bf v} = - \nu \bnabla^4{\bf v}$ represents the momentum dissipation, with ν the constant viscosity coefficient; T is the temperature; c_p is the specific heat at constant pressure; $\dot{q}_{\rm {net}}$ is the net diabatic heating rate; and ${\cal D}_T = - \nu \bnabla^4T$ represents the temperature dissipation.

Equations (1) are closed with the ideal gas law, p = ρRT, as the equation of state, with R the specific gas constant. A suitable set of boundary conditions, used in this work, is Dp/Dt = 0 at the top and bottom p-surfaces. Hence, the boundaries are material surfaces and no mass flow is allowed to cross the boundaries. With these boundary conditions, the equations admit the full range of motions for a stably stratified atmosphere—except for sound waves. For a discussion of the various aspects of the primitive equations (including superviscosity) and their use for extrasolar planet application, the reader is referred to Cho et al. (2003, 2008) and J. Y-K. Cho & E. Gulsen (2010, in preparation). In this work, as described below, Equations (1) are actually solved in a more general coordinate system.³

To solve Equations (1) in the spherical geometry, we use the Community Atmosphere Model (CAM 3.0). CAM is a well tested, highly accurate pseudospectral hydrodynamics model developed by the National Center for Atmospheric Research (NCAR) for the atmospheric research community (Collins et al. 2004). For hydrodynamics problems not involving sharp discontinuities (e.g., shocks) and irregular geometry, the pseudospectral method is superior to the standard grid and particle methods (e.g., Canuto et al. 1988).

As in many pseudospectral formulations of the algorithm, CAM solves the equations in the vorticity divergence form in the horizontal direction, where $\bzeta= \bnabla\times {\bf v}$ is the vorticity and $\delta = \bnabla \,{\bm\cdot}\, {\bf v}$ is the divergence. In the vertical direction, CAM uses the generalized p-coordinate:

$$\begin{equation} p(\lambda,\phi,\eta,t)\ =\ A(\eta) p_r\, +\, B(\eta) p_s(\lambda,\phi,t)\,, \end{equation} \tag{ 2 }$$

where λ is the longitude, ϕ is the latitude, η is the generalized vertical coordinate, p_r is a constant reference pressure, p_s(λ, ϕ, t) is a deformable pressure surface at the bottom boundary, and A, B ∈ [0, 1]. In the vertical direction, CAM uses the finite differencing method. Superviscosity (∇⁴ operators), as well as a small Robert–Asselin time filter (Robert 1966; Asselin 1972), is applied at every time step in each layer to stabilize the integration. The time stepping is done using a semi-implicit, second-order leapfrog scheme. Note that effects of various numerical dissipations are often subtle and can be significant on the integration, particularly over long times (e.g., Dritschel et al. 2007). Further details of the model and the effects of numerical viscosity on the flow evolution will be described elsewhere.

In all the simulations discussed in this paper, the physical parameters chosen are based on the close-in extrasolar planet, HD209458b. The basic result presented—that the evolution depends on the initial flow state—does not change for a different close-in planet. The physical parameters for the model HD209458b planet are listed in Table 1.

CAM is able to include radiatively active species and their coupling to the dynamics. However, we do not include them in the present work so that the effects discussed are not obfuscated by complications unrelated to the essential result. Our principle motivation is to study the dependence on the initial flow in the most unambiguous way possible. To this end, the flow is forced using the simple Newtonian drag formalism, as in many previous studies of extrasolar planet atmospheres (e.g., Cooper & Showman 2005; Langton & Laughlin 2007; Showman et al. 2008; Menou & Rauscher 2009). This drag is a simple representation of the net heating term in Equation (1d):

$$\begin{eqnarray} &&\frac{\dot{q}_{\rm {net}}}{c_p} = {-}\frac{1}{\tau _{\rm th}}\left(T - T_e\right), \end{eqnarray} \tag{ 3 }$$

where T_e = T_e(λ, ϕ, η, t) is the "equilibrium" temperature distribution and τ_th is the thermal drag time constant.

In this work, both T_e and τ_th are prescribed as barotropic (∂/∂η = 0) and steady (∂/∂t = 0), although simulations relaxing these restrictions have been run to verify robustness of our results. In general, both T_e and τ_th are (as are R and c_p) complicated functions of space and time (Cho 2008). Here,

$$\begin{equation} T_e\ =\ T_m + \Delta T_e \cos \phi \cos \lambda \,, \end{equation} \tag{ 4 }$$

where T_m = (T_D + T_N)/2 and ΔT_e = (T_D − T_N)/2 and T_D and T_N are the maximum and minimum temperatures at the day and night sides, respectively. Most of the simulations described in this paper have T_D = 1900 K, T_N = 900 K, and τ_th = 3 HD 209458 b planet days (where τ_p ≡ 2π/Ω is 1 planet day). Note that we have varied the timescale of the forcing by using a τ_th value in the range from 0.01 to 10 planet days, as well as letting the timescale decrease with height. The main result does not change for values ≳0.1 day, which nearly covers the entire spectrum of τ_th in all past studies using the Newtonian drag formalism.

The spectral resolution in the horizontal direction for most of the runs described in this paper is T42, which corresponds to 128 × 64 grid points in physical space.⁴ We have performed runs with resolutions varying from T21 (64 × 32) to T85 (256 × 128), in order to check convergence of the solutions. The vertical direction is resolved by 26 coupled layers, with the top level of the model located at 3 mbar.⁵ The pressure at the bottom η boundary is initially 1 bar, but the value of the pressure changes in time. This range of pressure is chosen because it encompasses the region where current observations are likely to be probing and where most of the circulation modeling studies have thus far directed their attention. We have also performed simulations in which the domain extends down to 100 bars and again verified that the basic behavior described in this paper is not affected. The entire domain is initialized with an isothermal temperature distribution, T_m = 1400 K.

To examine the robustness of evolved flow states to organized initial flow configurations, we have performed simulations with a wide range of initial conditions. The conditions from four of those runs (labeled RUN1–RUN4) are shown in Figure 1. In all the runs presented, the setup is identical—except for the initial flow configuration. The physical and numerical parameters/conditions are given in Tables 1 and 2, respectively.

Zoom In Zoom Out Reset image size

RUN1 is initialized with a small, random perturbation introduced in the flow. Specifically, values of u and v are drawn from a Gaussian random distribution centered on zero with a standard deviation of 0.05 m s⁻¹. RUN2 is initialized with a zonally symmetric, eastward equatorial jet of the following form:

$$\begin{equation} u_0(\phi)\ =\ U \exp \left\lbrace \frac{(\phi - \phi _0)^2}{2 \sigma ^2}\right\rbrace, \end{equation} \tag{ 5 }$$

where u(t = 0) =  u₀, U = 1000 m s⁻¹, ϕ₀ = 0, and σ = π/12. RUN3 is initialized with a westward equatorial jet described by Equation (5), with U = −1000 m s⁻¹, ϕ₀ = 0, and σ = π/12. RUN4 is initialized with a flow containing three jets. Note that the condition for RUN4 is very similar to the zonal average of the wind field of RUN1 at 50 planetary rotations. The jet profiles presented in Figure 1 are independent of height as well as longitude.

Figure 2 shows the temperature and flow⁶ fields of the four runs at t/τ_p = 40 (or t/τ_th ≈ 14). The fields near the 900 mbar pressure level are shown. (Recall that the η level surfaces of our model are functions of pressure, as described by Equation (2).) The figure illustrates the major point of this paper: given different initial states, there are clear, qualitative (as well as quantitative) differences between the different runs. Qualitatively, there are some common features. For example, most of the runs exhibit a coherent quadrupole flow structure—two large cyclonic and anti-cyclonic vortex pairs straddling the equator.⁷ However, the location of an individual vortex is different in the runs—as is the temperature pattern. In RUN3, a distinct quadrupole pattern is not present but there are more vortices in this run compared to the other runs. The temperature distributions are different because they are strongly linked to the flow. Consequently, the minimum-to-maximum temperature ranges vary from a moderately large 550 K (RUN4) to only about 200 K (RUN3) in the figure.

Zoom In Zoom Out Reset image size

The behavior just described is not restricted to a single altitude. Figure 3 shows the fields corresponding to those presented in Figure 2, but at a higher altitude (p ≈ 85 mbar pressure level). Comparison of Figures 2 and 3 illustrates the structural differences in 3D (vertical) as well as in 2D (horizontal). In RUN2 and RUN4, the large-scale vortices are strongly aligned, forming columns through most of the height extent of the modeled atmosphere; that is, the flow is strongly barotropic. In the other two runs, the flow is not vertically aligned throughout in large parts of the modeled atmosphere—and, therefore, the flow is baroclinic. The two figures also point to the corresponding strong difference in 3D temperature distributions, associated with the flow structures.

Zoom In Zoom Out Reset image size

This is more clearly seen in Figure 4. The figure shows the vertical (height–latitude) cross section of the temperature at 0°(sub-stellar) longitude from the runs presented in Figures 2 and 3. In Figure 4, the hottest and coldest regions are at different locations in all the runs. Near the equator, RUN1 exhibits a strong temperature inversion,⁸ while RUN2 and RUN4 do not. In addition, RUN2 and RUN4 exhibit generally strong decreases in temperature with height, while the others do not. As can be seen, the degree of temperature mixing varies strongly in the vertical direction among the runs—from ∼200 K contrast (RUN3) to ∼500 K contrast (RUN4). In general, the vertical structure is of low order, containing usually a single inversion. In our study, some form of inversion appears to be a generic feature.

Zoom In Zoom Out Reset image size

Preliminary steady state analysis of the primitive equations suggests that the basic behavior described above is due to the way in which the applied, (s, n) = (1, 1), forcing projects onto the normal modes of the planetary atmosphere; here, s is the zonal wavenumber and n is the total (sectoral) wavenumber of the spectral harmonics. In particular, a normal mode decomposition of the atmosphere into the vertical structure and Hough functions (e.g., Chapman & Lindzen 1970; Longuet-Higgins 1967) indicates that the forcing projects mostly onto low-order baroclinic modes, when the initial state is at rest. In contrast, when the initial state contains large-scale jets, the forcing projects more strongly on the barotropic mode, compared to the runs started from rest. Similar behavior has been observed in studies of the Earth's troposphere under tropical forcing (e.g., Geisler & Stevens 1982; Lim & Chang 1983). A more detailed study of coupling between forcing and normal modes is currently being performed and will be described elsewhere.

Furthermore, it is important to note that all of the above features, both dynamical and thermal, can vary in time. All of these features are important for observations (e.g., Knutson et al. 2007) and spectral modeling (e.g., Tinetti et al. 2007). A thorough study of the long-time evolution (over 1000 planetary rotations or more than 330τ_th) of the runs reveals a fundamental difference in their temporal behavior as well. For example, the flow pattern in RUN2 is characterized by two vertically aligned vortex columns in each hemisphere that translate longitudinally around the poles. The temperature in the upper altitude region is more strongly coupled to the flow than it is in the lower altitude regions. The flow pattern in RUN4 is also a set of vertically aligned vortex columns, but the columns oscillate in the east–west direction. The patterns in RUN1 and RUN3 are more complex, exhibiting a mixture of vortex splitting and merger and stationary states at different altitude levels. Figure 5, which shows a time series of the total kinetic energy, gives a quantitative measure of the different temporal behavior of the simulations.

Zoom In Zoom Out Reset image size

When the flow field is time averaged over a long period, the time-mean state is also significantly sensitive to the initial flow. This can be seen in Figure 6, which shows temperature cross sections at an arbitrary longitude (λ = 135°), averaged over 450 planetary rotations (planet day 300–750). The figure clearly shows that the variability observed is not simply a result of a phase shift in a quasi-periodic evolution. The flow and temperature structures in each run are fundamentally different from one another.

Zoom In Zoom Out Reset image size

Interestingly, the full range of flow and temperature behavior described above has been previously captured qualitatively, using the one-layer equivalent-barotropic model (Cho et al. 2008). The equivalent-barotropic equations are a reduced, vertically integrated version of the full primitive equations used in this study (Salby 1989). In many situations, the reduced model can be fruitfully used to study the dynamics of the full model by varying the Rossby deformation radius to represent different heights (or temperatures) of the full multi-layer model (e.g., Cho et al. 2008; Scott & Polvani 2008), and this also appears to be the case for hot exoplanets.

As might be expected from general nonlinear dynamics theory, in fact, the evolution can be strongly sensitive to small differences in the initial flow state. This is illustrated in Figure 7. There, two simulations are presented (panels (a) and (b)), which are identical in all respects except for a minute difference in the initial flow. The simulation in the top panel (RUN5) is started from rest. In contrast, the simulation in the bottom panel (RUN1) is started with a small perturbation: the initial values of u and v at each grid point are set to a Gaussian random distribution, centered on zero with a standard deviation of 0.05 m s⁻¹. Note that the maximum initial wind perturbation magnitude is only about 0.02% of the typical root mean square flow speed in the frames shown (∼500 m s⁻¹).

Zoom In Zoom Out Reset image size

The two panels in Figure 7 show temperature and flow distributions after t/τ_p = 1000 (or t/τ_th ≈ 333) at the p ∼ 420 mbar level. The distributions in the two panels are clearly different. At times they may look more similar than shown here, but in general that is not the case. The two runs generally show a different temporal behavior. Note that in Figure 7(a), there is a high degree of hemispheric symmetry, particularly in the north–south direction. In contrast, Figure 7(b) shows a clear asymmetry in both the north–south and east–west directions. The small asymmetry in Figure 7(a) is entirely due to machine precision and is not physical, since there is no way to break the symmetry in the setup of the run. Therefore, not surprisingly, some mechanism for inducing a noticeable symmetry breaking is necessary. The salient point here is, however, that even a tiny perturbation can lead to a marked difference in the flow and temperature distributions, even at relatively early integration times.

It is important to understand that the dependence on the initial flow state is robust and the behavior is not limited only to the parameters and the ranges, discussed thus far. The dependence has been verified for numerous model parameters and ranges. For example, Figure 8 shows that the strong dependence on the initial wind exists for much shorter τ_th, despite the strong forcing such drag times entail. In this case, τ_th/τ_p = 0.5. The upper panel is from RUN6 and the lower panel is from RUN7. The former run is initialized with only small stirring and no organized jet. In contrast, the latter run is initialized with a westward equatorial jet, identical to the setup of RUN3. At the shown time and height, there are clear differences between the flow and temperature patterns of the two simulations. The coldest area is advected east of the anti-stellar point in RUN6, but west of the anti-stellar point in RUN7. Furthermore, all the vortices have different locations. In RUN7, there is a fairly zonally symmetric jet at high latitudes, leading to a much more homogenized temperature distribution above the mid-latitudes than in RUN6.

Zoom In Zoom Out Reset image size

While a weaker difference might be expected in this case, based on other studies (e.g., Cooper & Showman 2005), the sensitivity is unabated in our simulations, and this holds for the long time average behavior as well. For a time mean over 300 rotations (e.g., planet days 1200–1500), at the level shown in Figure 8, the location of the coldest region differs by 40°in longitude between the two simulations shown in the figure. Figure 9 shows a time series of the total kinetic energy in the two simulations, revealing their different evolution. Even after integrating for a very long time (15,000 rotations), we have checked that the differences between RUN6 and RUN7 are still present. Note that this is a much longer integration time than that reported in any published studies of close-in planet circulation thus far. However, the flow at such long times is inexorably affected by cumulative numerical dissipation and phase errors and the result obtained should not be taken too literally (e.g., Canuto et al. 1988).

Zoom In Zoom Out Reset image size

In addition, we have studied the dependence when the initial jet contains a vertical shear—i.e., the jet distribution is baroclinic. Simulations have been initialized with an eastward jet, which has zero magnitude at the bottom and increases linearly with height so that the lateral flow distribution in the top layer of the model is identical to that in RUN2.⁹ In the baroclinic case, vortex columns evolve that extend throughout most of the atmosphere, as in RUN2. However, unlike in RUN2, the temperature structure here is also strongly barotropic. This appears to be related to the way the forcing projects onto the free modes of the system, as described in Section 3.1.

As alluded to in Section 2.1, the timescale of the fastest motions admitted by Equations (1) with the given boundary condition is ∼40 minutes, for the planetary physical parameters used in this work. Hence, a forcing with timescales of the same order or smaller is not entirely physically self-consistent with the use of Equations (1). Notably, such forcing introduces numerical difficulties. For the physically unrealistic drag time of τ_th = 1 hr, for which excessive numerical dissipation is required to prevent the code from blowing up or being inundated with numerical noise, the temperature evolution is only mildly sensitive to the initial condition, at long-time integrations. However, when such short drag times are only applied after a limited range in pressure—i.e., τ_th is allowed to vary with η—the sensitivity to the initial condition is present.