CIRCULATION AND DISSIPATION ON HOT JUPITERS

Strongly irradiated extrasolar giant planets—"hot Jupiters" —are expected to rotate synchronously with their orbits but to have strong longitudinal winds that redistribute heat from the day to the night side. The efficiency of redistribution is important for direct observables including infrared phase curves, the depth of secondary eclipses in transiting systems (i.e., eclipses of the planet by its star), and planetary spectra. Turbulence associated with the winds may contribute to chemical mixing of the atmosphere (Spiegel et al. 2009), and might even inject heat into the convective interior of the planet (Guillot & Showman 2002), thereby perhaps explaining why some of these planets have lower densities than expected by standard evolutionary models at their estimated ages.

In a previous paper (Goodman 2009, hereafter Paper I), one of us has argued that these thermally driven winds can be understood as natural heat engines, which convert a fraction of the thermal power into mechanical work: namely, the work expended to accelerate the wind. As with any other heat engine, the continual production of work must be balanced by mechanical dissipation, else the kinetic energy in the winds would grow without bound. Paper I offered a brief discussion of possible dissipation mechanisms. Here we explore those mechanisms in greater detail. The goal is to lay the foundation for subgrid models of the dissipative process suitable for use by global circulation codes.

The outline of this paper is as follows. Section 2 gives an overview of the dissipation mechanisms we consider. Section 3 presents a linear analysis of hydrodynamic Kelvin–Helmholtz (KH) instabilities of thermally stratified shear flows. We show that such flows can be destabilized by thermal diffusivity even at Richardson numbers greater than the well-known critical value Ri_crit = 1/4. Under conditions relevant to exoplanet winds, however, we estimate that the associated growth rates are too slow and the longitudinal wavelengths too long to be important for dissipation, unless Ri < 1/4, which requires transonic flow unless the vertical width of the shear layer is small compared to a pressure scale height. Hence shocks may be as important as shear-driven turbulence for dissipation. To investigate this, as described in Section 4, we have used the ZEUS code in two dimensions to simulate a part of the atmosphere with horizontal and vertical dimensions comparable to the pressure scale height. Thermal diffusion and viscosity have been added to the code, and thermal driving of the wind is simulated by a horizontal body force with nonzero curl. We study the velocity and dissipation rate of the wind in a statistical steady state as a function of the strength of the driving compared to the acceleration of gravity. Section 5 discusses our results in the context of previous work on winds in Jovian planets both within and beyond the solar system. Section 6 summarizes our main conclusions.

We begin by reviewing the physical conditions in those parts of the atmosphere where strong winds may occur. Because of the strong stellar irradiation, the temperature should be approximately constant with depth,

$$\begin{equation} T\approx 1300 \left(\frac{f L_*}{L_\odot }\right)^{1/4}\left(\frac{a}{10 \,R_\odot }\right)^{-1/2}{\rm \,K}, \end{equation} \tag{ 1 }$$

where L_* is the stellar luminosity, a the orbital semimajor axis, and f a dimensionless factor of order unity that depends upon the local zenith angle of the star, the albedo, and the efficiency of lateral heat redistribution (f = 1 for complete redistribution and negligible albedo). For a solar-mass primary, the orbital period is P_orb = 3.66(a/10 R_☉)^3/2 days. The sound speed in an atmosphere dominated by molecular hydrogen is c_s ≈ 2.75(T/1300 K)^1/2 km s^-1. A characteristic circulation timescale is

$$\begin{equation} t_{\rm circ}= \frac{2\pi R_P}{c_{\rm s}}\approx 1.85\left(\frac{T}{1300{\rm \,K}}\right)^{-1/2}\left(\frac{R_p}{R_J}\right){\rm \,d}\mbox{ays}\,, \end{equation} \tag{ 2 }$$

where R_p is the planetary radius. Assuming that the planet rotates synchronously, by this definition t_circ is comparable to the rotation period, as it is on Earth. Coriolis effects are likely to be important but not as dominant as on a more rapidly rotating and colder planet such as Jupiter itself.

Several considerations point to transonic wind speeds, as discussed in Paper I. First, absent dissipation, flow along isobars from the day to the night side, starting from rest, would accelerate to

$$\begin{equation*} U\sim \sqrt{g\Delta z}\ln \frac{T_{\rm day}}{T_{\rm night}}\,, \end{equation*}$$

where Δz ∼ H_p = c²_s/γg is the range of depths over which the thermal forcing changes sign (see Equation (16)). Indeed, repeated cycling between day and night could accelerate the wind indefinitely, were there no mechanical dissipation. Second, the strong entropy stratification tends to stabilize the flow against turbulent dissipation at small Mach numbers. And third, many numerical models of exoplanetary circulation do predict supersonic winds (see Paper I for references).

Another direct consequence of the radiatively induced temperature (1) is that the vertical density profile is approximately exponential, with a scale height $H_p=k_{\rm B}T/gm_{H_2}$ that is small compared to the planetary radius:

$$\begin{equation} \frac{H_p}{R_p}\approx 2.9\times 10^{-3}\left(\frac{T}{1300{\rm \,K}}\right) \left(\frac{R_p}{R_J}\right)^{-1}\left(\frac{g}{g_J}\right)^{-1}\,, \end{equation} \tag{ 3 }$$

where g = GM_p/R²_p is the surface gravity and $g_J\equiv {\it GM}_J/R^2_{J,\rm eq}\approx 2500\,{\rm cm\,s^{-2}}$.

The pressure at the infrared photosphere is

$$\begin{equation} p_{\rm ph}=\frac{2g_p}{3\kappa _{\textsc {r}}} \approx 1.7 \left(\frac{\kappa _{\textsc {r}}}{{\rm \,cm^{2}\,g^{-1}}}\right)^{-1}\left(\frac{g}{g_J}\right){\rm \,mbar}\,, \end{equation} \tag{ 4 }$$

where $\kappa _{\textsc {r}}$ is the Rosseland mean opacity appropriate to the temperature (1). A larger characteristic pressure is that at which the thermal time is comparable to t_circ. Taking t_th = H²_p/χ, where χ = 16σT³/κρ²c_P is the thermal diffusivity and $c_P\approx 3.5k_{\rm B}/m_{H_2}$ the specific heat capacity, we estimate the latter pressure at

$$\begin{equation} p_{\rm th}\approx 70\left(\frac{\kappa _{\textsc {r}}}{{\rm \,cm^{2}\,g^{-1}}}\right)^{-1/2}\!\left(\frac{T}{1300{\rm \,K}}\right)^{5/4}\! \left(\frac{g}{g_J}\right)\!\left(\frac{R_p}{R_J}\right)^{1/2}{\rm \,mbar}\,. \end{equation} \tag{ 5 }$$

The Reynolds number of these winds is enormous. The dynamic viscosity of molecular hydrogen is (Stiel & Thodos 1963) ρν = 1.9 × 10⁻⁴(T/1000 K)^0.65 g cm^-1 s⁻¹, and therefore

$$\begin{equation} \mbox{Re}\equiv \frac{c_{\rm s}H_p}{\nu }\approx 5\times 10^8\left(\frac{p}{\,{\rm mbar}}\right) \left(\frac{T}{1300{\rm \,K}}\right)^{-0.15}\left(\frac{g}{g_J}\right)^{-1}\,. \end{equation} \tag{ 6 }$$

At pressures above a microbar at most, therefore, it is appropriate to regard the winds as inviscid.

By other dimensionless measures, however, these flows are less ideal. At these temperatures and pressures, thermal ionization of alkali metals yields magnetic Reynolds numbers Rm ≡ c_sH_p/η ∼ O(1), where η is the magnetic diffusivity. Batygin & Stevenson (2010) and Perna et al. (2010) have suggested that non-ideal MHD effects could exert a significant drag on the winds and perhaps even contribute to heating the convective interior in the presence of a background planetary magnetic field ≳10 G.

A dimensionless inverse measure of thermal (rather than magnetic) diffusion is the Peclet number

$$\begin{equation} \mbox{Pe}\equiv \frac{c_{\rm s}H_p}{\chi }\approx 0.2\left(\frac{p}{\rm mbar}\right)^2 \!\left(\!\frac{\kappa _{\textsc {r}}}{{\rm \,cm^{2}\,g^{-1}}}\!\right)\!\left(\!\frac{T}{1300{\rm \,K}}\!\right)^{-7/2} \!\!\left(\!\frac{g}{g_J}\!\right)^{-1}\!\!. \end{equation} \tag{ 7 }$$

Thus, at the pressure (5) where thermal and advection timescales may be comparable, Pe ∼ 10³, with no direct dependence on $\kappa _{\textsc {r}}$.

Another important dimensionless quantity is the Richardson number, which characterizes the relative strength of stratification and shear. The local Richardson number at altitude z is Ri(z) ≡ N²/(dU/dz)², where N(z) is the Brunt–Väisälä frequency (Equation (11)). In the limit of infinite Re and Pe, a necessary condition for incompressible shear-driven instability is that Ri <  1/4 somewhere in the flow (Drazin & Reid 1981). For a flow that changes smoothly by ΔU over a vertical distance Δz, a characteristic value for Ri in a constant-temperature background is, using Equation (11) and noting H_p = c²_s/γg

$$\begin{equation} \mbox{Ri} \approx \frac{(\gamma -1)}{M^2\gamma ^2}\left(\frac{\Delta z}{H_p}\right)^2\,, \end{equation} \tag{ 8 }$$

where Ma = 2ΔU/c_s is the Mach number associated with the amplitude of the shear profile. Thus since γ ≈ 7/5 for an atmosphere dominated by molecular hydrogen, Ma ≳ 0.9(Δz/H_p) is required for instability. This linear stability criterion can be related to an energy principle and therefore should govern nonlinear instabilities also. However, when the stratification is primarily thermal rather than compositional, as we expect for the wind zone of these atmospheres, its stabilizing influence can be undercut by diffusion of heat. Thus instability may be possible at Ri>1/4 if the Peclet number (7) is not too large. One of the goals of the present study is to quantify this statement for hot-Jupiter winds.

From this survey of characteristic physical conditions and dimensionless numbers, three possible sources of mechanical dissipation for the winds come to the fore: (1) MHD effects, probably requiring B ≳ 10 G; (2) shear-driven turbulence, requiring Ri < 1/4 and/or moderate Pe; and (3) shocks, requiring M>1 as we have defined it. Another possibility, raised in Paper I, is the double-diffusive Goldreich–Schubert–Fricke (GSF) instability (Goldreich & Schubert 1967; Fricke 1968), which we do not consider here for two reasons: first, although the linear GSF instability is axisymmetric, we believe that its saturation can be studied reliably only with high-resolution three-dimensional (3D) simulations, whereas we limit ourselves to two dimensions (2D); and second, because the instability shuts off when (in the usual cylindrical coordinates) ∂v_ϕ/∂z = 0 and ∂(ϖv_ϕ)²/∂ϖ ⩾ 0, which conditions are not inconsistent with strong vertical shears at equatorial latitudes. Nonlinear MHD effects present even greater numerical challenges. Therefore, we concentrate on Richardson turbulence and shocks.

We are interested in the hydrodynamic stability of a horizontal shear flow. We first consider the Boussinesq limit in which we neglect density variations except where coupled to gravity. The equations of hydrodynamics become

$$\begin{eqnarray} \frac{{d}\boldsymbol{v}}{{\it dt}}&= -{\bnabla }\tilde{p} +\theta {\hat{\bm e}}_z+\nu \nabla ^2\boldsymbol{v}, \end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray} \frac{{d}\theta }{{d}t}&= -N^2{\hat{\bm e}}_z\boldsymbol{\cdot }\boldsymbol{v}+\chi \nabla ^2\theta, \end{eqnarray} \tag{ 9b }$$

$$\begin{eqnarray} \boldsymbol{\nabla \cdot }\boldsymbol{v}&=0. \end{eqnarray} \tag{ 9c }$$

Here kinematic viscosity ν is taken to be a constant,

$$\begin{equation} \theta \equiv -g\left(\frac{\partial \ln \rho }{\partial \ln T}\right)_p\frac{\delta T}{\bar{T}} \rightarrow g\frac{\delta T}{\bar{T}}\qquad \mbox{(ideal gas)} \end{equation} \tag{ 10 }$$

is proportional to the fractional temperature variation,

$$\begin{eqnarray} N^2 & \equiv & \frac{\rho g^2}{p}\left[\left(\frac{\partial \ln T}{\partial \ln P}\right)_S -\frac{{d}\ln T/{d}z}{{d}\ln P/{d}z}\right]_0 \nonumber\\ &&\rightarrow \frac{\gamma -1}{\gamma }\frac{g}{H_p}\quad \mbox{(ideal gas and $T_0=$ constant)}\quad \end{eqnarray} \tag{ 11 }$$

is the square of the Brunt–Väisälä frequency of the background hydrostatic state, and the pressure scale height

$$\begin{equation} H_p\equiv \left|\frac{{d}\ln p}{{d}z}\right|_0^{-1}\rightarrow \frac{k_{\rm B}T_0}{\mu m_p g}=\frac{c_{\rm s}^2}{\gamma g}\,. \end{equation} \tag{ 12 }$$

We assume the optically thick diffusion limit for radiative transfer in which χ = 16σT³/3κc_pρ². We further pick the opacity κ ∝ T³ρ⁻² so that χ = constant. Note that in our energy equation, we assume no background temperature gradient. Other authors have studied the linear stability of the Boussinesq equations (see, e.g., Gage 1973; Dudis 1974; Maslowe & Thompson 1971; Koppel 1964; Miller & Gage 1972; Gage 1972), but they do not consider in detail the destabilizing effects of thermal diffusivity in the viscous, optically thick limit. We are interested in the growth rates in the weak viscosity regime and their relevance to our numerical simulations.

Assuming a plane-parallel background state with horizontal flow

$$\begin{equation} \theta _0=0=\tilde{p}_0,\qquad \boldsymbol{v}_0=U(z){\hat{\bm e}}_x, \end{equation} \tag{ 13 }$$

we can expand all perturbations with the dependence exp(−iωt + ik_xx + ik_yy) to obtain the linearized sixth-order system of equations

$$\begin{eqnarray} v_{1z}^{\prime \prime \prime \prime }&\,{=}\,\left(\frac{ik^2\sigma }{\nu }-k^4 \,{-}\,\frac{ik_xU^{\prime \prime }}{\nu }\right)v_{1z} \,{+}\,\left(2k^2\,{-}\,\frac{i\sigma }{\nu }\right)v_{1z}^{\prime \prime }\,{+}\,\frac{(k^2)}{\nu }\theta _1,\nonumber\\ \end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray} \theta _1^{\prime \prime }&=\left(-\frac{i \sigma }{\chi }+k^2\right)\theta _1+\frac{N^2}{\chi }v_{1z}, \end{eqnarray} \tag{ 14b }$$

where σ = ω − k_xU, k ≡ (k²_x + k²_y)^1/2, and primes denote d/dz.

We use a relaxation method to solve this system of equations, and we test our results against limiting cases of background velocity profile U and parameters ν, χ, and N² for which there are simple analytic solutions. See Appendix A for the limits and analytical solutions. We can nondimensionalize our parameters by defining the Richardson number as Ri = N²H²/U²_max, Reynolds number Re' = HU_max/ν, Peclet number Pe' = HU_max/χ, and Prandtl number Pr = ν/χ = Pe'/Re', where H denotes a characteristic shear width. (In Section 2, we scaled Re and Pe by the sound speed rather than the flow velocity in the expectation of a Mach number of order unity, but that is not appropriate in the present Boussinesq context.) Our assumptions for energy transport and the shear profile allow us to give all our linear results in terms of these dimensionless numbers.

For the problem of interest, we specify boundary conditions for the system (14) as follows. We impose v_1z = 0 at z = ±z_max so that the fluid perturbations do not penetrate the boundary. We further impose v''_1z = −ik_xv'_1x − ik_yv'_1y = 0, corresponding to vanishing viscous stress, and θ₁ = 0, corresponding to fixed temperature, at both boundaries.

Consider now the limit ν → 0, χ → 0 of Equations (14)

$$\begin{equation} v_{1z}^{\prime \prime }=\left(-\frac{k_xU^{\prime \prime }}{\sigma }+k^2-\frac{(k^2)N^2}{\sigma ^2}\right)v_{1z}. \end{equation} \tag{ 15 }$$

From WKB arguments, the time taken for a wavelike disturbance to reach z = 0 from z = z₁ ≠ 0 is infinite if Ri>1/4. Further, from energy arguments the potential energy barrier separating two fluid elements at different heights exceeds the available kinetic energy that can be liberated if Ri>1/4 (Chandrasekhar 1961). From these considerations, we expect the stratified shear flow to be stable for Ri ≳ 1/4 in the inviscid case with no thermal diffusion. As we turn up the thermal diffusivity the potential energy barrier between fluid elements at different heights can be more easily overcome, and we expect instabilities to survive to higher Richardson number. We expect there to be a critical Richardson number, above which the stratified flow is stable and below which there exist unstable KH-like modes. We refer to this critical solution as the marginally stable mode.

For our stability analysis, we pick background velocity profile U = U₀tanh(z/H). Squire's theorem tells us that in the limit of vanishing thermal diffusivity the largest growth rates are obtained for k_yH = 0 (Drazin & Reid 1981), so we might expect this result to generalize to the case with thermal diffusion. In Figure 1, we show the marginally stable Richardson number as a function of k_xH for k_yH = 0, 0.1, 0.2 and Pe' = 10, Re' = 10³. As |k_yH| increases, the critical Richardson number decreases for all k_xH. Growth rates at fixed Richardson number correspondingly decrease. We set k_yH = 0 hereafter because we are interested in the fastest growing modes. Note that we find two distinct instabilities at long and short wavelengths. The long-wavelength instabilities have not converged with z_max. The solid lines indicate solutions computed with z_max = 3H, and the dashed curves indicate solutions computed at k_yH = 0 and z_max = 5H,  10H. The long-wavelength instability is driven to longer wavelengths as z_max increases. We run into similar issues encountered by Dudis (1973) in the long-wavelength limit. We must integrate to larger and larger z_max or use an appropriate asymptotic boundary condition as k_xH → 0 to obtain better convergence.

Zoom In Zoom Out Reset image size

Although our code has not converged for long wavelengths, our results suggest that there may be a real inviscid instability in the limit k_xH → 0. This is demonstrated analytically in Appendix B, The dimensionless form of the equations of motion given in Appendix B suggest that the long-wavelength modes should be sensitive to the combination of parameters Pe'Ri k z²_max/H, and the analysis guarantees the existence of the long-wavelength modes only where this combination is ≪1; some version of the long-wavelength modes may exist more generally, but we have not proven it. So it is not surprising that our numerical results for these modes have not converged with respect to z_max.

In Figure 2, we show the critical Richardson number as a function of k_xH at Re' = 10³ and at a sequence of Peclet numbers. As the Peclet number decreases below ∼40, the long-wavelength instability begins to strengthen and survive above Ri = 1/4. The dashed line is a run at higher z_max = 5H and shows the convergence of short-wavelength results. As the Peclet number decreases below ∼2.5, the flow at short wavelengths begins to destabilize above Ri = 1/4. For high Reynolds number and high Peclet number, our short-wavelength results are in rough agreement with the inviscid results of Dudis (1974). There are quantitative differences due to our assumption of an isothermal background temperature profile.

Zoom In Zoom Out Reset image size

In Figure 3, we show the growth rates of the modes with $\mbox{Ri}=\frac{1}{2}\mbox{Ri}_{\rm crit}$ as a function of k_xH, for fixed Re' = 10³ and at a sequence of Peclet numbers. The dashed curve shows the growth rate at fixed Richardson number Ri = 0.125 and Pe' = 2.5, Re' = 10³. The long-wavelength growth rates at k_xH = 3.2/Re' are well over two orders of magnitude smaller than shorter wavelength growth rates. This result agrees with our estimation that the inviscid instability in the limit k_xH → 0 has weak growth rate. For our numerical simulations, we will be primarily interested in the more unstable short-wavelength instabilities that have maximum growth at k_xH ∼ 0.5. Note that Ri_crit increases with decreasing Pe', so we are computing growth rates at higher Richardson number for the lower Peclet number curves. At fixed Ri, growth rates increase with decreasing Pe', as we expect from simple heat diffusion arguments.

Zoom In Zoom Out Reset image size

We conclude from Figures 2 and 3 that the linear growth rates of the fastest growing mode are unchanged in the parameter regime we are interested in, Pe' ∼ 10²–10⁴.

We use ZEUS-2D v2.0 (Stone & Norman 1992) to run hydrodynamical simulations of a horizontal wind profile in a local equatorial section of the atmosphere. Cartesian coordinates x and z correspond to longitude and altitude, respectively, each extending over a few pressure scale heights: typically 20H_p × 5H_p. Longitude x is periodic, though with a period much less than the true circumference of the planet. The simulations include vertical gravity, thermal diffusion of heat, shear viscosity, and a horizontal body force per unit volume $\rho f(z){\hat{\bm e}}_x$ (Equations (17) and (19)) rather than thermal forcing. A few words of justification are in order.

The pressure acceleration is nonpotential due to longitudinal entropy gradients; the component of its curl normal to the equatorial plane is, in spherical coordinates rθϕ,

$$\begin{eqnarray} {\hat{\bm e}}_\theta \cdot (\boldsymbol{\nabla }\rho ^{-1}\boldsymbol{\times }\boldsymbol{\nabla }p) &=& \frac{1}{r\rho ^2}\left(\frac{\partial \rho }{\partial T}\right)_p \left[\frac{\partial T}{\partial \phi }\frac{\partial p}{\partial r}- \frac{\partial p}{\partial \phi }\frac{\partial T}{\partial r}\right] \nonumber\\ &\approx & \frac{g}{r}\left(\frac{\partial \ln T}{\partial \phi }\right)_p\,, \end{eqnarray} \tag{ 16 }$$

in which the final form uses the ideal-gas law and assumes approximate vertical hydrostatic equilibrium. This curl must vanish in the (nearly isentropic) convection zone. Furthermore, since thermal forcing by itself does not add any net momentum to the flow,¹ the nonpotential longitudinal force must reverse sign with depth within the radiative atmosphere. In fact one expects this nonpotential "battery" to be concentrated at pressures ≲p_th (5), which is quite shallow compared to the radiative–convective interface, typically ∼kbar. The nonpotential force must also reverse sign at least twice in longitude because Equation (16) is periodic in ϕ, but the horizontal scale of variation of this force is probably small compared to its vertical scale since the radiative atmosphere is thin. Hence the profile (19), which is horizontally constant and extends over ∼2H_p vertically, reversing sign within this depth, is a plausible representation of the horizontal forcing within our local computational domain.

4.1. Setup

The source terms for momentum and energy in ZEUS are modified by the addition of the terms within square brackets below:

$$\begin{equation} \rho \left(\frac{\partial \boldsymbol{v}}{\partial t}\right)_{\rm sources}= -\boldsymbol{\nabla }p -\nabla \cdot \boldsymbol{Q}+ \left[\rho f(z){\hat{\bm e}}_x + \boldsymbol{\nabla \cdot \sigma ^{\prime }} -\rho g {\hat{\bm e}}_z\right] \end{equation} \tag{ 17 }$$

and

$$\begin{eqnarray} \left(\frac{\partial e}{\partial t}\right)_{\rm sources}&=& -p \boldsymbol{\nabla \cdot v} -\boldsymbol{Q\cdot \nabla v} \nonumber\\ &&+ [\boldsymbol{\nabla }\cdot (\chi \rho c_p \boldsymbol{\nabla }T) + \boldsymbol{\sigma ^{\prime }\cdot \nabla v}]\,, \end{eqnarray} \tag{ 18 }$$

in which $\boldsymbol{\sigma ^{\prime }}$ is the viscous shear tensor, with components

$$\begin{equation*} \sigma ^{\prime }_{ik}=\nu \rho \left(\frac{\partial v_i}{\partial x_k}+\frac{\partial v_k}{\partial x_i}-\frac{2}{3}\delta _{ik}\frac{\partial v_l}{\partial x_l}\right)\,, \end{equation*}$$

where νandg are taken to be constants. ZEUS' standard von Neumann & Richtmyer artificial viscosity tensor $\boldsymbol{Q}$ is applied with shocks spread over ≈2 zones. The thermal diffusivity χ is again given in the optically thick limit by χ = 16σT³/3κc_pρ² = const; constant χ would correspond therefore to κ ∝ T³ρ⁻². This again allows us to scale our results by the dimensionless Peclet number. The horizontal forcing profile is

$$\begin{equation} f(z)=\epsilon g[2\,\mathrm{sech}^2(z/H_p)\tanh (z/H_p)+\alpha \ \mathrm{sech}^2(z/H_p)], \end{equation} \tag{ 19 }$$

where H_p is the pressure scale height, α is a parameter that we adjust to ensure zero net momentum input on our grid, and ≪ 1 determines the overall strength of the forcing. In three dimensions a single latitude need not conserve momentum, but our simulations are limited to two dimensions. This form for the horizontal forcing is motivated by the viscous acceleration on a fluid with hydrostatic density profile and horizontal velocity profile given by v_x = U₀tanh(z/H_p). The viscous acceleration in this case is given by g = νU₀/H²_p and α = 1. We use an FTCS differencing scheme for our modified source terms, and the ZEUS time step constraint must satisfy δt < δt_ν = (Δx)²/4ν and δt < δt_χ = (Δx)²/4χ for stability. We are interested in Pr = ν/χ < 1, so the latter constraint is a tighter restriction on δt. But δt_χ/δt_CFL ≈ Pe/4N_p, where N_p is the number of grid points per pressure scale height and Pe is the Peclet number in the flow. Dimensionless parameters for our numerical results are scaled to the sound speed c_s and pressure scale height H_p as in Section 2. c_s here refers to the sound speed at the walls. Our simulations generally satisfy Pe ∼ 10–100 and N_p ∼ O(10), so δt_χ/δt_CFL ∼ O(1) and the diffusion and CFL time steps are usually comparable.

As stated, x is periodic. At z = z_min and z_max, we impose reflecting conditions for density and velocity and constant temperature T = T_wall. We initialize the density to ρ(x, z) = ρ₀ and the internal energy density to e(x, z) = e₀, allowing the flow to settle under the influence of gravity. This allows us to initialize the flow with velocity $\vec{v}=U_0\tanh (z/H_p){\hat{\bm e}}_x$ and vanishing net momentum. The wall temperature is set to the initial grid temperature, T_wall = T₀ = (e₀/ρ₀)μm_p(γ − 1)/k_B. We pick ρ₀ = 1, e₀ = 0.01, γ = 1.4, and g = 0.004. Our simulations are run over 5 vertical pressure scale heights, from z = −2.5H_ptoz = 2.5H_p. The density at the base of our simulation if it were in hydrostatic equilibrium at temperature T₀ would therefore be ρ_b = 5(1 − e⁻⁵)ρ₀. The code units for our problem are given in Table 1. We nondimensionalize parameters according to these units.

Table 1. Code Units

Unit	Length	Time	Density	Acceleration
Value	Hp = 1	t = Hp/cs = 13.36	ρb = 5	g = 0.004

Download table as:  ASCIITypeset image

We run our simulations over 20 scale heights in the horizontal direction. We are physically interested in Pe ≈ 10²–10⁴, Re⁻¹ → 0, but the viscosity must be large enough to maintain energy conservation, and ⩽ 10⁻². In the absence of a shear viscosity, the grid loses energy to numerical dissipation at the grid scale. We must pick ν large enough such that the viscous dissipation length l_ν ∼ (|∂_jv_i|/ν)^−1/2 is larger than the grid scale, and the dissipation of solenoidal turbulence is converted to internal energy. Our standard runs have N_p = 10 grid points per pressure scale height, but for certain runs which conserve energy less well we use up to N_p = 30 grid points per pressure scale height.

4.2. Tests

4.3. Diagnostics of Dissipation

We must also maintain energy conservation in our simulations. We are continually injecting energy onto the grid via the horizontal forcing, and to reach a statistical steady state this energy must first be converted into internal energy and then allowed to diffuse through the boundary walls. Because ZEUS is not based on a conservative form of the energy equation, numerical diffusion dissipates kinetic energy on the grid scale without converting it to internal energy; to avoid this, it is necessary to add an explicit shear viscosity (in addition to the artificial bulk viscosity used to mediate shocks).

We compute a fractional energy error

$$\begin{equation} \frac{\Delta \dot{E}}{\dot{E}_{\rm in}}=\frac{\dot{E}_{\rm in}-\dot{E}_{\rm kin}-\dot{E}_{\rm int}-\dot{E}_{\rm grav} -\dot{E}_{\rm out}}{\dot{E}_{\rm in}}, \end{equation} \tag{ 20 }$$

where

$$\begin{eqnarray*} \dot{E}_{\rm in}&=& \int\!\!\!\int \rho v_x f(z) dxdz \dot{E}_{\rm kin}=\frac{d}{dt} \int\!\!\!\int { {1 \over 2}}\rho v^2 dxdz \dot{E}_{\rm int}\nonumber\\ & =& \frac{d}{dt} \int\!\!\!\int e \,dxdz \\ [1ex] \dot{E}_{\rm out} &=& \int \left[-\rho c_p\chi \frac{\partial T}{\partial z}\right]_{z_{\min }}^{z_{\max }}\,dx \dot{E}_{\rm grav}= \frac{d}{dt}\int\!\!\!\int \rho g z\, dxdz. \end{eqnarray*}$$

We briefly mention here a number of strategies that we use to improve energy conservation. The fractional energy error will tend to decrease if we keep the amplitude of the forcing function lower so that there is less energy injection onto the grid. In a similar vein, if we localize the forcing function to the shear layer as opposed to distributing it over the entire vertical extent of the grid, then there will be less energy injection and better energy conservation. Other strategies include increasing the kinematic viscosity or making the grid resolution finer in order to capture kinetic energy dissipation above the grid scale via the explicit viscosity. Finally, we can increase the thermal diffusivity so as to allow excess internal energy to diffuse off the grid. These considerations must be balanced against keeping the computation time reasonably low and keeping the solutions physically relevant to the problems of interest.

We consider that a statistical steady state is reached when $\langle \dot{E}_{\rm int}\rangle _t$, $\langle \dot{E}_{\rm kin}\rangle _t$, and $\langle \dot{E}_{\rm grav}\rangle _t$ are all $\ll \dot{E}_{\rm in}$, and when the Mach number (Equation (24)) and the rms vertical velocity 〈v²_z〉^1/2_x,y are constant for at least a horizontal sound crossing time, 20H_p/c_s. In steady state, $\dot{E}/\dot{E}_{\rm in} \approx (\dot{E}_{\rm in}-\dot{E}_{\rm out})/\dot{E}_{\rm in}$. We maintain the fractional energy error to within 10 %.

Dissipation of kinetic energy involves production of entropy. If s is entropy per unit mass, then

$$\begin{eqnarray*} \rho \frac{Ds}{Dt} &=& \frac{-\boldsymbol{Q:\nabla v} +\boldsymbol{\sigma ^{\prime }:\nabla v}}{T} + \chi c_p\rho \left|\boldsymbol{\nabla }\ln T\right|^2\nonumber\\ && +\boldsymbol{\nabla \cdot }(\chi c_p\rho \boldsymbol{\nabla }\ln T). \end{eqnarray*}$$

When integrated over our computational domain, the final term is $-\dot{E}_{\rm out}/T_{\rm wall}$, $\dot{E}_{\rm out}$ being the rate of escape of heat. So

$$\begin{eqnarray} \dot{E}_{\rm out}& =& T_{\rm wall}\! \int [T^{-1}(-\boldsymbol{Q:\nabla v} + \boldsymbol{\sigma ^{\prime }:\nabla v}) + \chi c_p \rho |\boldsymbol{\nabla }\mbox{ln} T|^2]d^2\boldsymbol{x}\nonumber\\ &&+\frac{d}{dt}\int\!\!\!\int \rho s dxdy. \end{eqnarray} \tag{ 21 }$$

In statistical steady state, the time average of the last term vanishes.

While conversion of mechanical energy to heat ultimately depends upon the "microscopic" processes in Equation (21), we gain insight into the macroscopic processes that remove kinetic energy from the large-scale flow—shocks and turbulence—by monitoring the energy change in each operator split source step. The price paid for this approach is having to deal with a mixture of reversible and irreversible effects. Under adiabatic conditions, for example, the compressive heating term $\dot{E}_{pdV}=\int -p\boldsymbol{\nabla \cdot }\boldsymbol{v}\, d^2\boldsymbol{x}$ would be completely reversible, but shocks and thermal diffusion make $\langle \dot{E}_{pdV}\rangle \ne 0$ even in steady state. There is also a change in gravitational potential energy $\dot{E}_{\rm grav}=\int v_z \rho g d^2\boldsymbol{x}$ associated with the vertical gravity source step. In many of our runs, we find that $\dot{E}_{pdV}$ and $\dot{E}_{\rm grav}$ have large oscillations that nearly cancel. Using integration by parts, one can show that

$$\begin{equation} \dot{E}_{pdV}+\dot{E}_{\rm grav}= \int\!\!\!\int \boldsymbol{v}\cdot \left[\rho \boldsymbol{\nabla }(gz)+\boldsymbol{\nabla }p\right]\,dxdz\,. \end{equation} \tag{ 22 }$$

Thus the near cancellation suggests that our flows are close to vertical hydrostatic equilibrium despite large vertical motions. We find it easiest to monitor the combination $\dot{E}_{pdV}+\dot{E}_{\rm grav}$. The time average of this term is largely the work done to mix the fluid vertically; absent thermal diffusion, this could lead to an almost adiabatic temperature–pressure profile, ∂ln T/∂z ≈ −(γ − 1)g/c²_s, but in steady state diffusion tends to restore the stratification. Thus vertical mixing acts as a refrigerator—a heat engine in reverse—absorbing mechanical work to drive heat up the temperature gradient.

Shock dissipation is dominated by the artificial bulk viscosity, $\dot{E}_{\rm abv}=\int -\boldsymbol{Q\cdot \nabla v} d^2\boldsymbol{x}$ but receives small contributions also from the shear stress and from the $\dot{E}_{pdV}$ term.

Since we find that shocks are weak, we can estimate the turbulent energy dissipation rate as

$$\begin{equation} \langle \dot{E}_{\rm turb}\rangle _t\approx \langle \dot{E}_{pdV} +\dot{E}_{\rm grav}\rangle _t+\langle \dot{E}_\nu \rangle _t-\langle \dot{E}_{\rm lam}\rangle _t, \end{equation} \tag{ 23 }$$

where $\dot{E}_{\nu }=\int \boldsymbol{\sigma ^{\prime } \cdot \nabla v}$ and $\dot{E}_{\rm lam}=\int \nu \rho (\partial \bar{v}_x/\partial z)^2d^2\boldsymbol{x}$. Assuming an effective turbulent viscosity acting on the time-averaged mean shear profile, we estimate the turbulent dissipation rate as $\langle \dot{E}_{\rm turb}\rangle _t=\langle \int \nu _t \rho (\partial \bar{v}_x/\partial z)^2d^2\boldsymbol{x}\rangle _t$. The effective shear viscosity due to turbulence in a flow that does not resolve the turbulent energy cascade is then $\nu _{\rm t}=\nu \langle \dot{E}_{\rm turb}\rangle _t/\langle \dot{E}_{\rm lam}\rangle _t$.

4.4. Results

Figures 4, 5, and 6, respectively, give the run of Mach number, defined in Equation (24), rms vertical velocity, and fractional energy error versus Re for Pe = 11, = 0.1. The Prandtl number is less than one for all of these runs. Note the presence of four distinct solutions at Re < 10³. There is a fast solution, denoted by filled triangles, a medium velocity solution, denoted by filled squares, and two slow solutions, denoted by filled circles and open stars. Where solution branches overlap, we trace each branch by starting from a representative member and then gradually varying parameters. We use a similar parameter deformation to verify that each branch is limited to a certain range of Re as indicated in Figures 4–6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The solutions differ in the importance of laminar, turbulent, shock, and irreversible compressional dissipation. The fast solutions (triangles) have high laminar energy dissipation. Figure 5 indicates that there is a significant vertical component to the velocity field, but the bulk of the viscous dissipation is laminar rather than turbulent. This is reflected by wavelike oscillations in the flow as opposed to turbulent KH rolls. The medium and slow B solutions have stronger irreversible compressional heating and shocks, but they still have relatively high laminar dissipation. For the planetary wind problem of interest to us this laminar dissipation is unphysical, and as we increase the Reynolds number these solutions disappear. We conclude that these unphysical solutions are supported by the shear viscosity. Numerical diffusion can also support these unphysical solutions. At a grid resolution of N_p = 3 grid points per pressure scale height and parameters from Figure 4, we find unphysical solutions even for Re → ∞. We estimate that it is necessary to have N_p ≳ O(10) in order to enter the high-Re regime where there exists a single, physical solution. This unique solution is the slow A solution in Figures 4–6. It is characterized by weak laminar dissipation and strong irreversible compressional and shock dissipation. We check the convergence of this solution with runs at 1.5 and 2 times the standard resolution. Extra filled circles at Re = 3.03 in Figures 4–6 indicate higher resolution runs. The wind speeds and rms vertical velocities overlap, but the energy errors progress downward at higher resolution.

We summarize these results in Table 2 by providing quantitative diagnostics for each solution type. We give the shear viscosity, turbulent viscosity, relative importance of artificial bulk viscosity, relative importance of turbulent dissipation, Mach number, and Richardson number for a representative solution of each type. The laminar and turbulent energy dissipation rates are directly proportional to the quoted shear and turbulent viscosities, as described in Section 4.3. At this stage we are only providing rough diagnostics to characterize each solution type, but we infer that shocks contribute weakly to the irreversible compression because the artificial-bulk-viscosity dissipations are relatively small. The Mach number is computed as

$$\begin{equation} \mbox{Ma}\equiv \max _{z}\left[\frac{\int \rho (x,z)v_x(z)dx}{c_{\rm s,wall}\int \rho (x,z)\,dx}\right] -\min _{z}\left[\frac{\int \rho (x,z)v_x(z)dx}{c_{\rm s,wall}\int \rho (x,z)\,dx}\right]\,. \end{equation} \tag{ 24 }$$

We compute the minimum Richardson number in our flows via

$$\begin{eqnarray*} \mbox{Ri}&=& \min _z \frac{N^2}{(\partial \bar{U}/\partial z)^2} \nonumber\\ &=& \min _z\left[ \frac{\bar{p}}{\bar{\rho }}\left(\frac{\partial \ln \bar{p}}{\partial z}\right)^2\left(\frac{\partial \bar{U}}{\partial z}\right)^{-2}\left(\nabla _{\rm ad}-\frac{d\ln \bar{T}}{d\ln \bar{p}}\right)\right], \end{eqnarray*}$$

where overbar denotes horizontal average. The fastest and medium solutions are convective near the top wall. This is due to the higher wind speed and energy injection, which cause a steeper temperature gradient there.

Table 2. Viscosity, Turbulent Viscosity, Relative Importance of Artificial Bulk Viscosity, Relative Importance of Turbulent Dissipation, Mach Number, and Richardson Number for a Characteristic Solution of Each Type in Figure 4

Solution Type	log10Re	νt/csHp	$\langle \dot{E}_{\rm abv}\rangle _t/\langle \dot{E}_{\rm in}\rangle _t$	$\langle \dot{E}_{\rm turb}\rangle _t/\langle \dot{E}_{\rm in}\rangle _t$	Ma	Ri
Fastest (triangle)	2.17	0.006	0.07	0.43	7.6	<0
Medium (square)	2.57	0.010	0.22	0.60	5.5	<0
Slow A (circle)	3.03	0.015	0.13	0.79	2.3	0.033
Slow B (star)	2.40	0.017	0.12	0.67	2.9	0.024

Download table as:  ASCIITypeset image

Figure 4 shows that the fastest triangle solution in the laminar regime has wind speed scaling roughly as 1/ν, as we expect from our tests in Section 4.2. Further, the perfectly laminar solution at Re = 1.03 has Ri = 0.36. A flow with these dimensionless parameters is linearly stable in the Boussinesq limit. From our runs at lower forcing we determine that for the fastest (triangle) solution, the viscosity that maximizes the wind speed scales in proportion to the forcing: ν_peak/c_sH_p ∝ . The peak wind itself decreases only marginally. We can understand this behavior by assuming a steady-state balance between laminar dissipation of the mean flow and the forcing (19). These solutions are able to reach such a high wind speed because they are nearly laminar, until turbulent dissipation kicks in at low ν.

We map the solution space at lower forcing and thermal diffusivity and verify that there are four different solutions in the regime Pr ≲ 1. Specifically, we looked at Pe = 11 and = .025, .05, and0.1, and also at Pe = 21, = 0.1. These parameters correspond to the four values in the first three rows of Table 3, discussed in greater detail below. There are small quantitative differences in the results at Pe = 21, = 0.1. The fastest solution in particular has peak horizontally averaged Mach number ∼15% higher. From our linear analysis in the Boussinesq limit, we expect little difference in the results for any Peclet number in the range Pe'>10. In our simulations we are continually injecting energy into the fluid, however, and the diffusivity affects how this energy is redistributed. We attribute the differences in the nonlinear behavior for Pe ∼ O(10) to the important dynamical effect of this internal energy redistribution.

Table 3. Mach Number for a Sequence of Slow A Solutions at Different Peclet Number and Horizontal Forcing

	Pe
	11	21	75
0.1	2.3	2.3
0.05	1.9
0.025	1.7
0.01	1.7	1.6	1.7
0.005	1.5	1.7	1.7

Download table as:  ASCIITypeset image

The main takeaway point from these runs at different forcing and thermal diffusivity is that in the high-Reynolds-number limit there is a unique solution, for which the direct viscous dissipation of the mean flow is negligible compared to turbulent and shock dissipation. Further, if the grid resolution is insufficiently low, with N_p ≲ O(3), we may lock onto an artificially fast solution, even as Re → ∞. We check the latter result for the range of parameters we studied and not just for those relevant to Figure 4.

We now attempt to extract an effective turbulent viscosity from the physically relevant slow A solution. First, note that the energy errors in Figure 6 are largest for the slow A solutions. In order to decrease the numerical dissipation, we can either increase the shear viscosity ν or we can make the grid resolution finer. We cannot increase the shear viscosity indefinitely because we will enter the multiple-solution regime (Figure 4) or the laminar viscosity will become large enough to modify the properties of the solution. We must balance this against the increased computational cost of running the simulations at higher resolution. To calculate the turbulent viscosity, we use sufficiently high grid resolution and shear viscosity such that the numerical energy dissipation rate is ≲5% of the energy injection rate, but shear viscosity is low enough such that the gross properties of the solution are not modified.

In Tables 3 and 4, we give the time-averaged Mach number and turbulent viscosity for a series of slow A solutions at different horizontal forcing amplitudes and Peclet numbers. The exact value of the explicit viscosity we use is unimportant, so long as the viscous dissipation length is larger than the grid scale and the viscosity is small enough to leave the gross properties of the slow A solution unmodified. In general the runs satisfy Re>10³. The runs at Pe = 11 and = 0.025,  0.01 are done at the standard resolution of N_p = 10 grid points per pressure scale height. The rest are done at twice the standard resolution (i.e., twice as many points in each dimension), with the exception of the run at Pe = 75, = 0.01, which is done at three times the standard resolution, and the run at Pe = 11, = 0.005, which is done at 1.5 times the standard resolution. We emphasize that these values of N_p indicate the grid resolution required to conserve energy within 5% (Equation (20)), not the resolution required to find the slow A solution in the limit Re → ∞. At the standard resolution of N_p = 10, we still find the slow A solution in the limit Re → ∞ for all parameters surveyed, but the fractional energy error may be as large as 0.3 for the lowest diffusivity runs. The energy errors tend to increase with stronger forcing and with higher Peclet numbers, presumably because larger gradients of kinetic or thermal energy lead to larger spatial differencing errors in ZEUS.

Table 4. Turbulent Viscosity ν_t/c_sH_p for a Sequence of Slow A Solutions at Different Peclet Number and Horizontal Forcing

	Pe
	11	21	75
0.1	0.015	0.017
0.05	0.010
0.025	0.0061
0.01	0.0024	0.0023	0.0024
0.005	0.0012	0.0011	0.0010

Download table as:  ASCIITypeset image

The Mach number decreases gradually with forcing amplitude for = 0.1 down to = 0.01 and is roughly independent of for = 0.005–0.01. It is also independent of the diffusivity strength for Pe ∼ 10–100, at all forcing amplitudes surveyed. The latter is consistent with our linear analysis, and we attribute this result to the lower energy injection rate as compared to the faster solutions from Figure 4. There is less thermal redistribution of heat in the slow A solutions.

Figure 7 shows the Mach number as a function of time for the run at = 0.01 and Pe = 75. The run has reached a statistical steady state after time ∼10H_p/c_s. The steady state is characterized by a drive-up phase in which the Mach number gradually rises, and a short slowdown phase in which the Mach number rapidly drops. The slowdown occurs over a time period τ ∼ 0.5H_p/c_s, which is comparable to the growth time of the fastest growing linear instabilities in the Boussinesq limit (see Figure 3). In Figure 8, we show the relative contributions of dissipation due to artificial bulk viscosity, turbulence, and laminar drag as defined in Section 4.3. The first two dissipation terms are confined to the slowdown phases. Shocks appear briefly during slowdown, but the majority of the kinetic energy dissipation is due not to shocks but to the two terms in Equation (22): that is, the energy expended to mix the fluid vertically and overcome the stratification, which is restored by thermal diffusion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now characterize the nature of the sharp slowdown phases that limit the wind speed in our flow. The drive-up phase is relatively smooth with small changes in the state variables across the flow. Figure 9 shows the velocity and temperature fields for the starred point from Figure 7. At a short time before reaching the peak wind speed, t ≈ 0.2H_p/c_s, a vigorous instability develops with three wavelengths across the box. Figure 10 shows the velocity and temperature fields for the boxed point in Figure 7, which represents the slowdown phase. The unstable mode is clearly visible in the temperature profile. Note the shocks at z/H_p < −1 and x/H_p = −4.5, 2, 8.5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dimensions of our box enforce a periodicity under the transformation x → x + 20H_p. If the unstable mode is an ordinary KH mode, then from our linear analysis the only modes that can grow have either 1, 2, or3 wavelengths in the box. Apparently the mode with 3 wavelengths and wavenumber k_xH_p ≈ 0.94 is the fastest growing of the three. The Richardson numbers satisfy Ri < 1/4. We verify that the growth of this instability and its effect on the wind speed are independent of the dimensions of our domain by increasing the horizontal periodicity length (L_x), with a corresponding increase in the number of grid points (N_x) so as to maintain the same resolution. At L_x = 30H_p the fastest growing modes have 3 or 4 wavelengths in the box, and at L_x = 40H_p they have 4 or 5 wavelengths. These wave numbers are roughly consistent with the linear instabilities in the Boussinesq limit. More importantly though, the time-averaged Mach number of the flow remains unchanged.

Ultimately we would like to provide a dissipation prescription for global simulations that do not resolve the turbulent nature of the wind profile. We first note that from Figure 7 we find a horizontally and time-averaged Mach number Ma ∼ 1.7 and peak horizontally averaged Mach number Ma ∼ 2.7. The peak occurs just before the rapid growth of instabilities, and so the latter corresponds to the absolute peak Mach number in the flow. Global simulations are often more strongly supersonic. For example, the model of HD 209458b by Cooper & Showman (2005) finds v_max ≈ 9 km s^-1 at the 2.5 mbar level (see Paper I for a survey of the calculations prior to 2009).

Most of these studies of global circulation employ reduced equations of motion based on the approximation of vertical hydrostatic equilibrium: shallow-water, equivalent barotropic, or primitive equations. This is done in order to cope efficiently with the large ratio of horizontal to vertical scales (Equation (3)), but the approximations used also have the effect of filtering out sound waves and shocks (even horizontally propagating ones). Such codes cannot directly represent the shocks and turbulent vertical motions studied here, although they can represent horizontal shears (which we do not) and hydraulic jumps (e.g., Showman et al. 2009; Rauscher & Menou 2009), which bear some resemblance to shocks but do not conserve energy and may not be triggered in the same circumstances. Therefore no particular reason exists as to why our results for the wind speeds should agree quantitatively with those obtained from the primitive equations, for example.

Comparison with the calculations carried out by Burkert et al. (2005), Dobbs-Dixon & Lin (2008), and Dobbs-Dixon et al. (2010, hereafter DCL) is potentially more revealing, because these calculations, like ours, use the full compressible hydrodynamic equations. The most recent of these studies included an explicit shear viscosity. Simulations were run for ν = 10⁴, 10⁸, 10¹⁰, and 10¹² cm² s^-1. The peak wind speed was found to decrease monotonically with increasing viscosity, though the differences between the two least viscous cases were slight. For reference, ν = 10¹⁰ cm² s^-1 corresponds to Re ≈ 800 referred to the sound speed and scale height, which were approximately 10³ K and 360 km, respectively, near the western terminator for these simulations. Thus DCL's simulations straddle the regime in which we find multiple solution branches, with the wind speed actually increasing with decreasing Re on the fastest branch (Figure 4). While there is no clear sign of such behavior in DCL's results, this is perhaps not surprising in view of their coarse sampling in ν.

As with our own simulations, it is not clear whether DCL's wind speeds are limited fundamentally by Mach number (i.e., shocks) or Richardson number (i.e., shear instabilities). Although DCL do not report Richardson numbers, we have attempted to estimate these from the profiles of azimuthal velocity (v_ϕ) versus pressure in their Figure 4. The profile for $\nu =10^{8}{\,\rm cm^2\,s^{-1}}$ at latitude 0° and longitude 70° (counterclockwise from "noon") jumps from a small value to $3.2{\,\rm km\,s^{-1}}$ (Mach ∼1.4) near the 0.05 mbar level over a vertical distance ∼0.3H_p. From the information given in the paper, we estimate that this corresponds to Ri ≈ 10⁻²: quite small compared to the critical value Ri ≈ 1/4. The radial resolution appears also to be ∼0.3H_p, however, so the shear layer is underresolved; also, the jump occurs near the outer boundary, somewhat above the infrared photosphere, while the profile at greater depths is much smoother. Hence it is not clear that the minimum Richardson number is meaningful. This illustrates the difficulty in resolving the vertical structure of these atmospheres in a global, 3D simulation based on the full gasdynamic equations. It is nevertheless reassuring that DCL's peak Mach numbers are quite similar to what we find here, M_max ∼ 2, in their least viscous cases.

A limitation of our local, 2D models is the somewhat arbitrary assumption that the vertical width of the nonpotential horizontal force, and hence that of the shear layer, is Δz ≈ 2H_p. Other things being equal, one would expect the flow speed to scale ∝Δz at constant forcing amplitude if Richardson number rather than Mach number is the controlling parameter. Table 5 gives Mach and Richardson numbers for different shear widths calculated at Pe = 11,   = 0.01. The time-averaged Mach number scales roughly linearly with the shear width for the range of shear widths explored here. Both the time-averaged and minimum Richardson numbers are roughly constant for different shear widths. For shear widths outside the range shown in the table, it is difficult to compare results because the shape of the shear profile is strongly modified. As noted above, the shear profile in DCL's simulations is quite different: a jump by more than the sound speed across a single cell (corresponding to ∼0.3H_p), followed by a roughly linear decline with depth from 10⁻⁴ to 1 bar. A change in sign of the horizontal flow is not always obvious in their velocity plots but is inferred from their temperature profiles.

We can use our turbulent viscosities to estimate a turbulent diffusivity for composition, K_t, as may be required to keep the upper atmosphere turbulently mixed (Spiegel et al. 2009). The Schmidt number, defined by Sc_t = ν_t/K_t, in general has a complicated dependence on flow properties (see, e.g., Huq & Stewart 2008). For Ri ≲ O(10⁻¹), however, Sc_t is independent of Ri and of order unity. The data in Table 4 then gives us mass diffusivities pertaining to the shear layers in hot-Jupiter atmospheres: K_t ≳ 10⁻³c_sH_p ∼ 10¹⁰ cm² s^-1. Spiegel et al. (2009) estimated that K_t ∼ 10⁷–10¹¹ cm² s^-1 is required, depending on grain size, to maintain TiO as an optical absorber at p ≲ 1 mbar. Our values are toward the upper end of this range, but only in the shear layer (∼2H_p). To estimate K_t elsewhere, we would need to estimate how much of the mechanical power dissipated in the shear layer can propagate upward and downward in the form of internal waves, acoustic waves, or shocks. We hope to address this in a future paper.

We have investigated the balance between horizontal forcing and mechanical dissipation in model shear flows designed to imitate thermally driven winds of Jovian exoplanets. Our main conclusions are as follows.

• 1.  
  Following Paper I, some form of mechanical dissipation is required to offset the production of mechanical energy by the longitudinal entropy gradient. Plausible dissipation mechanisms include shocks, turbulence, and MHD drag.
• 2.  
  Entropy stratification has a stabilizing influence on KH modes, but thermal diffusivity (χ) tends to undercut this influence to a degree that can be quantified by the Peclet number Pe ≡ c_sH_p/χ based on the sound speed (c_s) and pressure scale height (H_p). We estimate Pe ∼ O(10³) at the altitude where the thermal timescale is comparable to the circumferential flow time at Mach 1, but Pe ≲ O(1) near the infrared photosphere.
• 3.  
  Linear analysis indicates that for Pe' ≡ HU_max/χ ≳ 10, where U is the horizontal wind speed, thermal diffusion has little effect on the stability of KH modes whose horizontal wavelength is comparable to the vertical thickness of the shear layer; the usual adiabatic Richardson criterion Ri ≡ N²/(∂U/∂z)²>1/4 is sufficient for the stability of these short-wavelength modes. For typical conditions on the day sides of hot Jupiters, Ri < 1/4 implies transonic or supersonic winds, presuming that the shear layer is no narrower than a pressure scale height.
• 4.  
  Nevertheless, at any finite Pe' and Ri, inviscid flows with an inflection point in the shear profile are unstable at sufficiently long horizontal wavelengths, at least when the flow is confined within a channel of finite vertical extent. We doubt that these long-wavelength KH modes are important for hot-Jupiter winds because their growth rates are very small (inversely proportional to wavelength) and because the planetary circumference is not large enough compared to H_p to accommodate them at Pe' ≳ 10³.
• 5.  
  Nonlinear, 2D (azimuth and depth), horizontally forced local simulations saturate at time-averaged Mach numbers ∼2 when the Reynolds number Re ≡ c_sH_p/ν ≳ 10³. In this regime, which is relevant to hot Jupiters, we find that the saturated wind speed depends only weakly on Pe and on the strength of the forcing. There is, however, a systematic increase in the time-averaged Mach number with the vertical width of the imposed forcing, although the time-averaged Richardson number is more nearly constant (see below).
• 6.  
  Turbulent mixing associated with the nonlinear development of KH modes is the primary agent for dissipating kinetic energy. The timescale for the onset of the instability is roughly 0.5H_p/c_s. The Richardson number plays an important role in controlling the onset of the instability, and we find that the time-averaged 〈Ri〉_t ≈ 0.15 ± 0.02 and the minimum Ri_min ∼ 0.03–0.05 in the physically relevant regime, independent of shear width in the range 1.6–3 H_p.
• 7.  
  Dissipation in the high-Re regime is due mainly to the work expended to overturn the stratification rather than shocks, though weak shocks are present. It is highly intermittent, at least in our 2D calculations, and appears to be triggered by a recurrent linear instability of KH type.
• 8.  
  At Re ≲ 10³ or at low numerical resolution (significantly fewer than 10 cells per pressure scale height), multiple sequences of solutions appear in overlapping ranges of Reynolds number, differing in their wind speed and predominant dissipation mechanisms.
• 9.  
  We estimate a turbulent diffusivity for composition in the wind zones of hot-Jupiter atmospheres. For Pe ≈ 100, we find K_t/c_sH_p ≈ 10⁻³, depending on the strength of the horizontal forcing. However, these values pertain to regions of maximum vertical shear; turbulent diffusivities may be much smaller elsewhere.

Points 5 through 8 suggest that some global simulations of hot Jupiters may have overestimated wind speeds due to incomplete physics (e.g., neglect of vertical accelerations and sound waves) and/or insufficient spatial resolution. In fairness to these efforts, it must be noted that the characteristic length and timescales of the shear-driven instabilities that appear to dominate the dissipation are quite small, of order H_p ∼ 200 km and H_p/c_s ∼ 10³ s, compared to the global length and timescales of the atmospheric motions (R_p ∼ 7 × 10⁴ km, R/c_s or Ω⁻¹ ∼ 10⁵ s). This makes it extremely challenging—or perhaps impossible with current computational resources—to resolve these instabilities directly in global simulations. On the other hand, local simulations such as those here are incapable of correctly representing the global boundary conditions and thermal forcing of the flow; in particular, they cannot properly represent the moderation of the forcing as the advection reduces the day–night temperature contrast. Until sufficient resources become available to treat the local and global scales simultaneously, we suggest that artificial dissipative terms be added to the equations of motion in the global simulations so as to limit the Richardson number (from below, perhaps Ri>0.1) and perhaps also Mach number (from above, perhaps Ma < 3).

We thank Jim Stone for help with ZEUS. This work was supported in part by the National Science foundation under grant AST-0707373.

We provide here four limits of Equations (14a) and (14b) for which there analytic solutions. We test our code against these solutions.

Taking N² = 0, ν → 0, Equations (14a) and (14b) admit solutions satisfying the inviscid perturbation equation with vanishing temperature perturbation,

$$\begin{eqnarray} v_{1z}^{\prime \prime }&=\left(-\frac{k_xU^{\prime \prime }}{\sigma }+k^2\right)v_{1z}, \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} \theta _1&=0, \end{eqnarray} \tag{ A2 }$$

in which k ≡ (k²_x + k²_y)^1/2. For a hyperbolic tangent velocity profile U = tanh(z), the inviscid perturbation equation for the marginal ω = 0 mode becomes

$$\begin{equation} v_{1z}^{\prime \prime }=-(2 \mathrm{sech}^2z\,-k^2)v_{1z}. \end{equation} \tag{ A3 }$$

A solution to this equation that vanishes as z → ±∞ is

$$\begin{equation} v_{1z}=\mathrm{sech}\ z\,,\qquad k^2=1. \end{equation} \tag{ A4 }$$

In the limit

$$\begin{equation} U=0,\nu \rightarrow 0,N^2=0, \end{equation} \tag{ A5 }$$

Equations (14a) and (14b) become

$$\begin{equation} v_{1z}^{\prime \prime }=(k^2)v_{1z}-\frac{i}{\omega }{k^2}\theta _1, \end{equation} \tag{ A6 }$$

$$\begin{equation} \theta _1^{\prime \prime }=\left(-\frac{i\omega }{\chi }+k^2\right)\theta _1, \end{equation} \tag{ A7 }$$

which for Re(ω) = 0 and Im(ω) < −(k²)χ admit the decaying mode solutions

$$\begin{equation} \theta _1=A\cos \left[\left(-k^2+\frac{i\omega }{\chi }\right)^{1/2}z +\phi \right], \end{equation} \tag{ A8 }$$

$$\begin{equation} v_{1z}=A\frac{\chi k^2}{\omega ^2}\cos \left[\left(-k^2+\frac{i\omega }{\chi }\right)^{1/2}z +\phi \right]+Be^{{k}z}+Ce^{-{k}z}. \end{equation} \tag{ A9 }$$

In the limit χ → ∞ with conditions (A5), we obtain solutions to Equations (A6) and (A7) of the form

$$\begin{equation} \theta _1=Ae^{{k}z}+Be^{-{k}z}, \end{equation} \tag{ A10 }$$

$$\begin{eqnarray} v_{1z}&=&\frac{i}{4\omega }(Ae^{{k}z}+Be^{-{k}z}+2{k}z(Ae^{{k}z}-Be^{-{k}z}))\nonumber\\ &&+Ce^{{k}z}+De^{-{k}z}. \end{eqnarray} \tag{ A11 }$$

Finally, in the limit

$$\begin{equation} U=0,\chi \rightarrow 0, \nu \rightarrow 0, \end{equation} \tag{ A12 }$$

Equations (14a) and (14b) combine to give

$$\begin{equation} v_{1z}^{\prime \prime }={k^2}\left(1-\frac{N^2}{\omega ^2}\right)v_{1z}. \end{equation} \tag{ A13 }$$

For ω real and N²/ω²>1, Equation (A13) has solutions of the form

$$\begin{equation} v_{1z}=A\cos \left[k\left(\frac{N^2}{\omega ^2}-1\right)^{1/2}z+\phi \right], \end{equation} \tag{ A14 }$$

$$\begin{equation} \theta _1=A\frac{iN^2}{\omega }\cos \left[k\left(\frac{N^2}{\omega ^2}-1\right)^{1/2}z+\phi \right]. \end{equation} \tag{ A15 }$$

Here we demonstrate analytically for the Boussinesq equations that there exist long-wavelength instabilities at arbitrarily large Richardson number if the Peclet number is finite. It is helpful to adopt dimensionless variables. Supposing U(z) ≈ ±U₀ at |z| ≫ H, where H is the width of the shear later, we define these variables as follows:

$$\begin{eqnarray} \begin{array}{rllrllrll} z_* &=& z/H & U_*(z_*) &=& U(z)/U_0& \\ v_*(z_*) &=& v_{1z}(z)/U_0 & \theta _* &=& \frac{\theta _1}{N^2H} \\ k_* &=& H\sqrt{k_x^2+k_y^2} & \cos \alpha &=& k_x/\sqrt{k_x^2+k_y^2} & \\ c_* &=& \omega /(k_xU_0)&&&&\\ \mbox{Re} &=& U_0H/\nu & \mbox{Pe} &=& U_0H/\chi &\\ \mbox{Ri} &=& ({\it NH}/U_0)^2\,.&&&& \end{array} \end{eqnarray} \tag{ B1 }$$

Equations (14a) & (14b) become, with primes² denoting d/dz_*

$$\begin{eqnarray} &&\left[v^{\prime \prime }_*+\frac{U^{\prime \prime }_*}{c_*-U_*}v_*-k_*^2v_*\right] (c_*-U_*)\cos \alpha \nonumber\\ &&\quad -\frac{i}{k_*\mbox{Re}}\left(\frac{d^2}{dz_*^2}-k_*^2\right)^2v_* = -\mbox{Ri}\, k_*\theta _*, \end{eqnarray} \tag{ B2a }$$

$$\begin{eqnarray} \theta ^{\prime \prime }_*+\left[i \mbox{Pe}\,k_*\cos \alpha \,(c_*-U_*)-k_*^2\right]\theta _* &= \mbox{Pe}\, v_*. \end{eqnarray} \tag{ B2b }$$

We take the inviscid, long-wavelength limit in such a way that k_*Re → ∞ while k_* → 0⁺, with Pe and Ri regarded as finite and nonzero. The wavelength $2\pi |\boldsymbol{k}|^{-1}$ is not infinite, but rather ≫2πH, so that we may neglect the explicit appearances of k_* in the equations above yet regard c_* as finite. If c_* should turn out to have an imaginary part, then the growth rate Im(ω) → k_xU₀Im(c_*) as k_* → 0. The potential temperature θ_* decouples:

$$\begin{eqnarray} v^{\prime \prime }_*+\frac{U^{\prime \prime }_*}{c_*-U_*}v_* &=0\,, \end{eqnarray} \tag{ B3a }$$

$$\begin{eqnarray} \theta ^{\prime \prime }_* &= \mbox{Pe}\, v_*. \end{eqnarray} \tag{ B3b }$$

The first of Equations (B3) is the long-wavelength limit of the inviscid, unstratified KH problem with a continuous shear profile U_*(z_*). Given boundary conditions v_*(±z_*,max) = 0 for some z_*,max that is finite and sufficiently large, and presuming that U_* is an odd function of z_* so that U''_*(0) = 0, then there is a solution to (B3) for which c_* has a positive imaginary part. In fact, although it does not satisfy the boundary conditions, v_* ∝ U_* − c_* is a solution to Equation (B3a) for any c_*, and a second solution v_* = f(z_*)(U_* − c_*) can be found by solving a first-order equation for f'(z_*), followed by a quadrature. Requiring that this solution vanish at both ±z_*,max determines c_* via

$$\begin{equation} \int _{-z_{*,\max }}^{z_{*,\max }}\frac{dz_*}{[U_*(z_*)-c_*]^2}=0. \end{equation} \tag{ B4 }$$

Since U_* is odd and ≈±1 at large |z_*|, it follows easily that c_* ≈ ±i when z_*,max ≫ 1.

This demonstrates the existence of inviscid instability when $|\boldsymbol{k}|^{-1}\gg z_{*,\max }^2\gg 1$. If z_*,max ≳ k^−1/2_*, however, then a separate analysis would be needed, because the solution of Equation (B3b) might yield θ_* ∼ Pe z²_*,maxv_*, and then the right-hand side of Equation (B2a) could be O(1) or larger. We have not done this analysis because the numerical evidence indicates that the long-wavelength growth rates are too slow to be important for our applications.