© ESO, 2016

1. Introduction

In the past decade, hot-Jupiter observations have evolved from characterizing only the orbital (semi-major axis, period, eccentricity) and structural (mass and radius) properties of the newly discovered planets to also constraining the physical properties of the atmosphere such as the composition and the temperature-pressure profile (Madhusudhan et al. 2014; Heng & Showman 2015), as well as observables that are directly connected to the dynamics such as 2D brightness maps (Knutson et al. 2007; de Wit et al. 2012) and even (although maybe still tentatively) atmospheric winds (Snellen et al. 2010). With new instruments such as the James Webb Space Telescope (JWST) soon available to the community, future observations will become more and more constraining. Their interpretation and understanding will require elaborate and well-understood models of the dynamical structure of their upper atmosphere. This need explains the recent development by several groups of numerical models describing the atmospheric dynamics of hot Jupiters (among which Showman & Guillot 2002; Menou et al. 2003; Cho et al. 2008; Showman et al. 2008, 2009; Dobbs-Dixon & Lin 2008; Rauscher & Menou 2010; Heng et al. 2011b), a comprehensive review of which has recently been published by Heng & Showman (2015). Taken as a whole, these global circulation models (GCM) have produced a few robust results that characterize the dynamics of hot-Jupiter atmospheres. Probably the most notable is the existence of an equatorial jet with typical eastward gas velocity on the order of a few kilometers per second, as first shown by Showman & Guillot (2002). This is faster than the speed of sound by a factor of a few (i.e., the Mach number of the jet, defined as the ratio between the wind and the sound speed, is larger than unity). This remarkable result comes from the fact that hot Jupiters are both tidally locked and strongly irradiated by their central star as a result of the small semi-major axis of their orbit. This has been given a theoretical basis by the semi-analytical calculations of Showman & Polvani (2010, 2011), who were able to explain the origin of the equatorial jet by the excitation of a stationary planetary scale Rossby wave by stellar irradiation, a result recently confirmed and extended to include three-dimensional effects by Tsai et al. (2014). One reason of the success of these models is that they appear to be consistent with some of the observations of hot-Jupiter atmospheres, such as the eastward shift of the hot spot at the photosphere (Knutson et al. 2007) or the day-night temperature contrast measurements (Perez-Becker & Showman 2013).

Despite these significant findings, several questions remain. First, the supersonic nature of the equatorial jet has raised the question whether shocks exist in hot-Jupiter atmospheres. The presence of such structures could have important consequences for the atmosphere time variability and act as an efficient drag mechanism on the flow (Rauscher & Menou 2012). This question is made even more acute because many of the published results have been obtained using approaches borrowed from the community that studies the Earth atmosphere climate and uses approximations that are well adapted to study the Earth atmospheric dynamics. The most widely used approach is based on solving a reduced set of equations, the so-called primitive formulation of the hydrodynamics equations, which filters out vertically propagating sound waves and assumes that the gas is in hydrostatic equilibrium in the vertical direction. Such codes are only able to capture hydraulic jumps as opposed to real shocks (see Rauscher & Menou 2012, for a discussion of this issue). Using such approaches, several groups have reported the existence of regions where significant compression occurs and takes the form of standing discontinuities (Showman et al. 2009; Rauscher & Menou 2010), suggesting that shocks might indeed develop, in agreement with theoretical arguments such as given by Heng (2012). In addition to these models, a few papers have been published and presented simulations that do not rely on these standard GCM approximations and thus do not suffer from these limitations (Dobbs-Dixon & Lin 2008; Dobbs-Dixon et al. 2010; Mayne et al. 2014a,b). These groups solved the Navier-Stokes equations on the sphere and explicitly included some form of dissipation to stabilize the scheme and/or to account for the effect of unresolved features of the flow. Even though the issue of shocks is not the main focus of either of these papers, Dobbs-Dixon et al. (2010) reported shock-like features in some of their models. However, the use of high dissipation coefficients in these particular simulations significantly changes the mean flow so that the question of shock formation in hot-Jupiter atmospheres remains unanswered.

In addition, it is possible that the equatorial jet is destabilized by various hydrodynamic instabilities. For example, Li & Goodman (2010) have recently shown using analytical arguments coupled with idealized numerical simulations that fast supersonic equatorial jets are vulnerable to a vertical shear instability, a possibility first mentioned by Showman & Guillot (2002)¹. As shown by Li & Goodman (2010), shocks might develop during the nonlinear evolution of this instability and have significant consequences for the jet mean velocity. As explained by Li & Goodman (2010), properly resolving this instability in a GCM is a tremendous task because it requires the numerical grid to be sufficiently fine to resolve the vertical pressure scale height of the atmosphere in the horizontal direction. Given these difficulties and because of the very idealized nature of the numerical setup used by Li & Goodman (2010), it is perhaps not surprising that the question of the possible development and potential consequences of a vertical shear instability in hot-Jupiter atmospheres has never been addressed in any GCM so far.

Besides their role as possible drag mechanisms, the questions of the existence of shocks and/or of the growth of hydrodynamic instabilities such as described above may have some potential observational consequences in terms of the flow variability, but also to explain the origin of the inflated hot Jupiters (Baraffe et al. 2010). Indeed, it has been suggested that downward transport of kinetic energy possibly associated with these instabilities (Showman & Guillot 2002; Guillot & Showman 2002) along with turbulent mixing of heat (Youdin & Mitchell 2010) might transport and deposit a sufficient amount of thermal energy into the deep layers of the atmosphere to account for such highly inflated radii (for a detailed discussion, see Ginzburg & Sari 2015).

The purpose of this paper is to develop an idealized model for investigating the problems described above. For this purpose, we solve the compressible Euler equations. Since we are interested in shocks, we use the finite-volume shock-capturing scheme RAMSES (Teyssier 2002) that is based on the Godunov method (Toro 1997, see also Sect. 2.3). In addition to being well adapted to resolving shocks in supersonic flows such as encountered here, finite-volume codes like RAMSES conserve the gas total energy, the sum of its kinetic and thermal energy, even in the absence of explicit dissipation. This means in particular that all the kinetic energy that might be numerically dissipated (as would be the case, for example, in a turbulent flow) returns into heat. As noted by Goodman (2009), this is of particular importance in the context of hot-Jupiter atmospheres. Our strategy in this first paper is to keep the numerical setup as simple as possible while retaining the basic ingredients that drive the flow dynamics. We therefore simplify the thermodynamics and geometry in two important aspects. While most published results now use a complex but realistic treatment of radiative transfer effects, we model these effects with a simple parametrized cooling function. For practical purpose, we use a form of the cooling function that is linear in the temperature, also known as Newtonian cooling. This has proved very useful in modeling atmospheric flows in general and the general circulation of hot Jupiters in particular (Showman et al. 2008; Rauscher & Menou 2010; Heng et al. 2011b). Published simulations using this approximation compare well with more realistic models that include a detailed treatment of radiative transfer effects (Dobbs-Dixon & Lin 2008; Showman et al. 2009; Heng et al. 2011a; Rauscher & Menou 2012). Our second approximation is to solve the equations in a Cartesian geometry using the equatorial beta plane model. This is a well-known approximation in atmospheric dynamics that has been successfully used recently in the context of hot-Jupiter atmospheric flows (Showman & Polvani 2010, 2011; Tsai et al. 2014; Heng & Workman 2014). The idea is to focus on the equatorial region of the planet and to expand the vertical projection of its angular velocity linearly to the first order in the distance y to the equator: (1)where β is a constant and a free parameter of the model, while e_z stands for the unit vector in the vertical direction. For hot-Jupiter atmospheres, this is a particularly useful approximation because of their slow rotation and since the dynamics mainly develops in the vicinity of the equator. Using these simplifications but including the effect of compressibility, the goal of this paper is to investigate questions such as the occurrence of shocks, the stability of the equatorial jet, and the origin of the flow variability, if any.

The plan of the paper is as follows. In Sect. 2 we detail our numerical setup and the parameters we use. The idea is to use a set of numerical parameters as close as possible to the benchmark models presented by Heng et al. (2011b) and to investigate the potential effects induced by compressibility in this model. In Sect. 3 we present a low-resolution model and show that we recover many of the features described in the literature despite our simplistic geometry. These results validate our approach and our approximations. In Sect. 4 we make a brief resolution study before presenting a high-resolution simulation in Sect. 5. We show that the flow features variability at different spatial and temporal scales that we relate to well-known hydrodynamics instabilities. We then conclude and discuss the limitations of our work and the perspectives it opens in Sect. 6.

2. Physical model and numerical implementation

2.1. Equations and notations

As explained in the introduction, we solved the hydrodynamic equations in a Cartesian coordinate system denoted by (e_x,e_y,e_z), in a frame rotating with angular velocity Ω. We considered a grid that extends over the ranges [−L_x/ 2,L_x/ 2 ], [−L_y/ 2,L_y/ 2 ] , and [ 0,L_z ] in the x, y, and z directions, respectively. The x-direction should be thought of as representing longitudes in traditional GCM models, y as a proxy for latitudes and z would stand for the radial direction (aligned with the gravitational acceleration). Since hydrostatic equilibrium is not built into the equations, we did not use pressure coordinates in the vertical direction as is common in atmospheric sciences when solving the primitive formulation of hydrodynamics. The time evolution for the gas density ρ, the velocity v and the total energy E writes

where g is the constant vertical acceleration due to the planet gravitational field, ℒ is the cooling function (see Sect. 2.2), and P denotes thermal pressure. This is related to the total energy with the relation (5)In the above equation, we have introduced the adiabatic exponent of the gas γ and assumed an ideal equation of state in writing the internal energy e as a function of thermal pressure. In this case, P also relates directly to the gas temperature T through the following equation: (6)Here k is the Boltzmann constant, μ is the mean molecular weight, m_H the mass of the hydrogen atom, and ℛ is the specific gas constant. Throughout this paper, we assume that the atmosphere is composed of a mixture of hydrogen and helium and take for ℛ the same numerical value as Heng et al. (2011b): ℛ = 3779 J kg^-1 m^-1. We also set γ = 1.4. Finally, as explained in the introduction, we work in the framework of the equatorial β-plane model, so that Ω is given by Eq. (1). In the remainder of this paper, we therefore refer to the y = 0 plane as “the equator”.

2.2. Cooling function ℒ

ℒ varies linearly with the departure from an equilibrium temperature, a prescription that is also known as Newtonian cooling: (7)where T_eq and τ_rad are the radiative equilibrium temperature and radiative cooling timescale, respectively. They both depend a priori on the location in the atmosphere. ℒ is meant to provide a rough but computationally straightforward description of the balance between heating due to stellar irradiation and radiative cooling. Our choice for the spatial dependence of τ_rad is discussed in Sect. 2.4 and Appendix A. To calculate T_eq, we followed a procedure that is largely inspired by (but not identical to, essentially because we use a different coordinate system) the benchmark calculation presented by Heng et al. (2011b). We describe it here briefly for completeness. We started from a reference pressure-dependent temperature profile that is usually computed from 1D radiative transfer calculations. We next defined (8)where ΔT is a constant and a free parameter of the simulation. For simplicity, we did not prescribe any pressure variation for ΔT, which is different from papers published in the literature that used Newtonian cooling to model radiative effects. Because the planet is tidally locked, the substellar point is fixed in time. We chose its horizontal coordinates to be (x,y) = (0,0). Using these definitions for T_day and T_night, we calculated T_eq as follows: on the night side of the planet (|x| >L_x/ 4), we followed Heng et al. (2011b) and set T_eq = T_night. On the planet day side (|x| <L_x/ 4), we used the relation (9)We note that this profile is different from the functional form of T_eq introduced by Heng et al. (2011b). While the sinusoidal x-dependence is reminiscent of their longitudinal profile, a similar dependence would have little meaning in the y-direction. It is more natural to use an exponential variation in that direction, as done recently by Showman & Polvani (2011). The drawback of this formulation is that it introduces an additional length scale L_th in the problem that is another free parameter of the model. The specific value we chose in this paper is discussed below in Sect. 2.4.

2.3. Numerical implementation

We solved the above equations using a uniform grid version of the community code RAMSES (Teyssier 2002). RAMSES uses a finite-volume scheme based on the MUSCL-Hancock Godunov method (Toro 1997). By systematically upwinding all the waves that enter the problem, finite-volume codes are intrinsically stable (provided they satisfy the Courant-Friedrich-Lax, or CFL, condition) and do not require using any explicit dissipation. In fact, the algorithm is constructed to add the minimum amount of dissipation that is necessary to stabilize the numerical scheme. This is important for one of the problems we are interested in, namely the issue of shock formation. However, a well-known difficulty of Godunov codes is their inability to properly handle vertical hydrostatic equilibrium (see, e.g., LeVeque 1998). This is because the latter results from a balance between two terms that are treated differently by the numerical algorithm, namely the pressure gradient (treated as a flux of momentum computed as part of the solution of a Riemann problem) and gravity (treated as a source term by a simple split Crank-Nicholson algorithm). Starting from a solution initially in exact equilibrium balance, this situation immediately creates a mass flux at cell interfaces that rapidly leads to spurious vertical velocities (even in 1D) and compromises the simulation. Recently, Käppeli & Mishra (2014) developed an algorithm for the MUSCL-Hancock scheme that is especially designed to fix this problem. The idea is to assume that the density and pressure profiles within each cell are isentropic and to use the extrapolated values at cell faces as inputs for the Riemann problem. We have implemented their solution in RAMSES and found that it gives very satisfactory results, in the sense that the code is now able to keep a stratified atmosphere in hydrostatic equilibrium with almost vanishing deviations.

The boundary conditions (BC) we used are periodic in the x-direction. Because the Coriolis force increases in amplitude with y, the dynamics is confined to the vicinity of the equator. For this reason, the boundary conditions in the y-direction are not critical provided the domain is chosen wide enough. In the simulations presented in this paper, we set the horizontal velocity to zero there and assumed a zero gradient for all other variables. The vertical boundary conditions are more subtle to implement and somewhat arbitrary. We empirically found that the following produces good results, in the sense that artifacts introduced by the boundaries could not be detected: we extrapolated both density and pressure assuming an isentropic vertical profile, imposed a zero horizontal velocity gradient BC and reflective BC on the vertical momentum. In addition, we forced the mass, momentum, and energy fluxes to vanishes at the vertical boundaries. We note that the two conditions (on the variables themselves and on the fluxes) are not necessarily consistent with each other, but did not lead to any problem in the simulations.

The numerical implementation described above was tested by reproducing standard results of the literature. We present two of them in the appendix of the paper, namely the growth of a baroclinic wave in an adiabatic atmosphere as described by Polichtchouk et al. (2014) in Appendix B.1 and a shallow hot-Jupiter model such as presented by Menou & Rauscher (2009), Heng et al. (2011b) and Mayne et al. (2014a) in Appendix B.2. Although these tests do not have an analytical solution, the similarity between our results and the published calculations, along with the analysis of a low-resolution deep model that we present in Sect. 3, give confidence in our setup.

2.4. Model parameters

The physical model presented above contains several free parameters: g and β characterize the planet, while , τ_rad, ΔT, and L_th describe the heating and cooling processes. Instead of presenting a complete survey of the associated parameter space, the idea of the present paper is to choose a unique set of parameters that matches those of Heng et al. (2011b) as closely as possible and to show that the flow properties are similar to the results presented in that paper. We defer a detailed study of the influence of these free parameters on the results to a future publication.

Planet parameters:

we followed the shallow hot-Jupiter model of Heng et al. (2011b) and considered a planet with a radius a_p = 10⁸ m for which g = 8 m s^-2. We assumed that the planet rotates with a frequency Ω_p = 2.1 × 10^-5 s^-1. Expanding Ω = Ω_psin(φ), where φ is the latitude, close to the equator, and using Eq. (1), we thus have (10)

Newtonian cooling parameters:

we next need to prescribe the parameters entering in the definition of T_eq and τ_rad. To do so, we used the analytic relation between and P and between τ_rad and pressure. These two functions are given in Appendix A and have numerical parameters chosen to give an approximate match to the profiles of Heng et al. (2011b). We also set ΔT = 300 K in all our simulations. The curves so obtained for T_night and T_day are plotted in Fig. A.1 and are close to the profiles used by Heng et al. (2011b, see their Fig. 7). It is worth noting that τ_rad becomes infinite below 10 bar. This choice creates an “inert layer” at the location where P> 10 bar, in the sense that gas cannot cool radiatively in this region. As we shall see, this point turned out to have rather important consequences in the simulations. Finally, the cooling function ℒ depends on the parameter L_th that governs its y-dependence. Apart from being on the order of a_p, there is no obvious way of choosing L_th, and its influence on the flow properties deserves further investigations. After a few tests at low resolution, we found a good agreement with published results for L_th = 0.7a_p. We use this value in the remainder of the paper.

Computational domain and initial conditions:

we need to specify the extent of the computational domain and the initial conditions. We used L_x = 2πa_p in the x direction and found that L_y = 2.5a_p was wide enough for the lateral boundaries to have no effect on the flow topology. As noted by Mayne et al. (2014a), the vertical extent of the box has to be chosen carefully because hot gas of the planet day side is inflated and expands upward. We thus followed these authors and used L_z = 9 × 10⁶ m. Except for one model (see Sect. 5), our calculations were all initiated from an atmosphere at rest (v = 0 everywhere), with a uniform temperature T = 1800 K and in vertical hydrostatic equilibrium. We set the density in the bottom layer (i.e., at z = 0) so that P = 220 bar. With these parameters, we typically found that the pressure at the top of the domain ranges between 6 × 10^-3 and 0.1 mbar, i.e., well below the typical pressures we are interested in in the remainder of the paper. We note that there is a current debate in the literature regarding the sensitivity of the flow to the initial conditions (Thrastarson & Cho 2010; Liu & Showman 2013), particularly when explicit dissipation is small or vanishing, as we consider here (Cho et al. 2015). Here, we leave these questions aside and consider the same initial conditions as used by Heng et al. (2011b).

Simulation post treatment:

unless otherwise stated, we typically produce snapshots of all the physical variables every planet day during the simulation. To facilitate comparison with published papers, we present all of our results interpolated in the vertical direction on a pressure grid that extends from 200 bar at the bottom to 1 mbar at the top using 48 levels uniformly spaced in the logarithm of the pressure.

Finally, we measured time in units of one planet day, which corresponds to 299 200 s with our choice of parameters, or about 3.5 Earth days. In the simulations we present below, the time step is always limited by sound waves propagating in the vertical direction. It is typically on the order of 25 s in model LowRes (Sect. 3) and 15 s in model HighRes (Sect. 5). This should be compared with the typical time steps of 120 s used by, for example, Heng et al. (2011b). Although markedly larger when using the primitive set of hydrodynamics equations, the difference with a fully compressible set of equations is not large enough to be prohibitive when performing hot-Jupiter atmosphere simulations. This is because the typical pressure scale-height H is large in this case (typically 10³ km) so that the ratio between H and the horizontal scales involved in the problem is not as small as it would be for an Earth-like simulation.

3. Low-resolution model

To validate our numerical implementation and our choice for the numerical parameters entering the problem, we start by presenting a model LowRes. The grid spatial resolution was set to (N_x,N_y,N_z) = (64,33,48) in this section. This is smaller by about a factor of two than the resolution used by the finite difference core of Heng et al. (2011b). However, we did not use any explicit dissipation here, so that a one-to-one correspondence is difficult to establish a priori. The model was integrated for 500 planet days, which corresponds to about 1700 Earth days for our choice of Ω_p.

We now describe the climatology of the atmosphere. As previously documented in the literature, the flow develops a strong eastward jet around the planet equator (Fig. 1, left panel). In about 200 days, it reaches steady state at pressure levels above 1 bar, while higher pressure levels gradually accelerate until the end of the simulation (Fig. 2). Based on this result, we computed the mean properties of the flow by calculating time averages from t = 200 until the end of the simulation. The zonally averaged jet velocity at the equator reaches about 6.6 km s^-1 at the top of the atmosphere (i.e., at a pressure level of about 1 mbar), which corresponds to a mean Mach number of about 2. As a result of the short radiative timescale in the atmosphere (compared to the dynamical timescale), the temperature structure displays a structure that is close to the zonally averaged radiative equilibrium temperature, except in the deep layers of the atmosphere (at P = 10 bar and below) where a hot spot is seen with T up to 2000 K. We find very little velocity fluctuations in the flow above 1 bar after 200 days, which is consistent with, for example, Showman et al. (2008). Even if we tend to find slightly faster equatorial jets, the zonally averaged structure described above compares favorably with results published in the literature that use a similar cooling function (Showman et al. 2008; Rauscher & Menou 2010; Heng et al. 2011b; Mayne et al. 2014a). One difference is, however, that we find weaker westward jets at high latitudes (with maximum westward winds of 100 m s^-1 compared to typical published values that are about one order of magnitude higher). This difference most likely arises because we used a Cartesian coordinates system, as opposed to the spherical geometry that is commonly used. It needs to be kept in mind when comparing our results with previously published models. We checked that the total (i.e., volume-integrated) angular momentum is conserved in our calculations. In the equatorial β-plane approximation, the latter takes the form (11)where the integral should be taken over the computational domain. We found that M is conserved to within 2% in model 128 × 33 over the 500 days of integrations (with the small change being associated with a small leakage through the meridional boundaries of the domain).

The day-to-night heating contrast results in strong zonal asymmetries in both the dynamical fields and the temperature. These asymmetries are stronger in the atmosphere upper layers and gradually decrease downward. This is illustrated in Fig. 3, which shows the temperature and horizontal winds in the (x,y) plane at the pressure levels 1.66 mbar, 100 mbar, and 4.4 bar. The zonal asymmetries in the first two panels are large because of the short radiative timescale at those locations. For example, at 1.66 mbar, the highest (resp. lowest) velocity amounts to 7.1 km s^-1 (resp. 5.5 km s^-1) and temperatures range from about 700 K to 1400 K. By contrast, the equatorial jet is essentially zonally symmetric at P = 4.4 bar, in agreement with previous results. At 100 mbar, we recover the chevron shape structure in the temperature reminiscent of previously published calculations. The structure of the flow at 1.66 mbar is interesting because it illustrates one of the consequence of using the equatorial β plane model: while simulations performed on the sphere typically display a planetary scale Hadley cell, with upward motions in the substellar point and downward motions at the antistellar point, this flow is impossible to obtain by definition in the framework of the equatorial β-plane model. Instead, meridional velocity are always deflected back to the equator since the Coriolis force becomes gradually stronger as y increases.

To conclude on Sect. 3, the good agreement between the flow structure in our simulation, as shown in Figs. 1 and 3, with equivalent figures reported by Heng et al. (2011b), along with the fact that these results broadly agree with several studies of the same kind published in the literature, validate both our numerical implementation and the choice of our model parameters.

4. Resolution study

We next investigated the sensitivity of the properties described above to variations in the horizontal spatial resolution. All of the other parameters of the simulations were kept fixed. We varied N_x from 64 to 1024 and N_y from 33 to 195. The simulation parameters and their outcomes are summarized in Table 1. All runs but one were integrated for 400 days. Model HighRes, with a resolution of (N_x,N_y) = (1024,195), is computationally demanding. To reduce part of the burden, we restarted model 1024 × 65 at t = 100 and multiplied the number of cells in the y direction by a factor of three. The equations were then integrated for another 350 days.

4.1. Effect on the climatology

All models qualitatively display the same climatology as described above. However, quantitative measures vary. As an example, we report in Table 1 (Col. 4, see also Fig. 4) the zonally and time-averaged (between t = 200 and t = 400) zonal wind at the equator at the pressure level P = 50 mbar. We note that the equatorial jet is still being slightly accelerated over that period and has not yet reached a perfect equilibrium. For example, in model 128 × 65, the zonal wind at 50 mbar increases from 7.1 to 7.4 km s^-1 during that period (which represents less than 5%, a rather small acceleration that justifies the claim made above that the flow is close to steady state). Taken as a whole, this resolution study shows a strong sensitivity to the meridional resolution N_y (see the blue squares), but a very weak sensitivity to the zonal resolution N_x (see the symbols that correspond to N_y = 65 and varying N_x, all clustered at a mean zonal wind of about 7.2–7.3 km s^-1). This difference comes from the different spatial scales of the equatorial jet in the x and y directions, namely a sharp meridional velocity gradient, but a large-scale longitudinal structure. In model LowRes, grid cells have a meridional extent of dy = 10⁷ m, which means that the jet is only resolved with a handful of grid points (see the jet meridional size in Fig. 1) and is severely affected by numerical dissipation. Quite differently, the zonal variations of the jet have a typical scale similar to the zonal computational domain itself (simply because thermal forcing is modulated on that scale) and is easily resolved with a few tens of cells in all of our models.

4.2. Shocks

As discussed in the introduction, it has been a long-standing problem to determine whether shocks exist in hot-Jupiter atmospheres. We did not find any signatures of shocks in our simulations that would be traced by a sudden decrease of the local Mach number or a rapid increase of the density. As a more quantitative diagnostic, we computed a dimensionless measure of the flow divergence, as done recently by Zhu et al. (2013): (12)where dx denotes the cell size in the radial direction. In a shock-capturing scheme such as RAMSES, the numerical algorithm is designed to spread shocks over only a handful of cells, so we would expect η to reach at least a few tens of percent in shocks and be independent of spatial resolution. Zhu et al. (2013) argued that η> 0.2 is a good criterion for detecting strong shocks in their simulation. Since shocks are spread over a few cells in finite-volume codes such as RAMSES, such a threshold indeed corresponds to a velocity jump larger than the sound speed across a discontinuity (or, using the Rankine-Hugoniot conditions, to an upstream Mach number larger than ~1.5).

In general, we found that η only amounts to a few percent at a pressure higher than a few tens of mbar, and decreases downward. For example, in model 128 × 65, its highest value reaches 0.045, 0.034 and 0.022 at pressure levels P = 50 mbar, 100 mbar, and 1 bar, respectively. η only takes significant values at pressure of a few mbar, as illustrated in Fig. 5. At such a low pressure, three zones feature values of η in excess of 10%. The first two are symmetric around the equator and located on the jet wings for positive values of x. We note that they also manifest themselves by an increased temperature (see Fig. 3). It is likely that these two compressive regions are an artifact of the equatorial β-plane model that prevents a hemispheric Hadley cell from developing around the planet. The last zone with a large η is located at x ~ −0.8 around the equator and corresponds to the region where the zonal wind is strongly decelerated when entering the day side of the planet. η reaches 23% at this location in model 128 × 65. As emphasized by Heng (2012), this is the location where shocks would most easily form. However, here, it appears that this zone is only a region of strong compression, but not a shock. Indeed, we found that η gradually decreases when the zonal resolution is increased, reaching 12% in model 256 × 65 and only 6% in model 512 × 65. This gradual decrease shows that the zonal-wind decline when entering the planet day side is smooth and converged upon increasing the resolution, as opposed to a shock that would steepen and be caracterized by a constant high η value. We note that the two symmetric high η regions in Fig. 5 show the same trend with resolution, which means that these two zones are also regions of adiabatic compression. We conclude from this analysis that the mean flow in the atmosphere of hot Jupiters displays no shocks.

4.3. Time variability

We also computed an estimate of the high-frequency rms fluctuations of the zonally averaged equatorial zonal wind at 50 mbar (see last column in Table 1 and green symbols in Fig. 4 for the models with N_x = 128). To do so, we first subtracted the low-frequency component of the wind by smoothing the raw data with a Hamming function with a width of 20 days. is then simply the root mean square of that signal. While it is very small when N_y = 33, we find that increases steadily and amounts to more than 100 m s^-1 for model 128 × 195 and about twice that value for our highest resolution model HighRes (not shown in Fig. 4). We now focus on this last model to investigate the physical origin of this variability.

5. High-resolution model

The nature of the zonally averaged zonal wind variability in model HighRes is best illustrated by showing the time variation of the vertical profile of the zonal wind at the equator (Fig. 6, top panel). The difference with Fig. 2 is striking: shows significant variations during the entire duration of the simulation. These variations take the form of quasi-periodic oscillations of the zonal wind, with a period of about ten days, during the first half of the simulation (t< 300) and are less regular, but still significant, at later times. The zonal wind variations are in phase from the top of the planet atmosphere down to the inert layer below 10 bar, which is most likely a result of the fast communication timescale across the atmosphere provided by sound waves. The time history of at P = 50 mbar is shown in Fig. 6 (bottom panel). We can distinguish three phases of the zonal wind evolution (indicated with labels in the upper part of both panels): phase I (t = 100 to t = 170) is a spin-up phase during which the flow adjusts to the sudden change in meridional resolution. In agreement with the results of Sect. 4, the zonal wind rapidly increases by about 1 km s^-1 over a period of just a few days. We call the period that extends between t = 170 and t = 270 phase II. The flow appears to have reach a quasi-steady state during this period (see also the top panel of the same figure, which suggests that this is the case at all levels above a few bar), and the zonal wind displays quasi-periodic oscillations around a well-defined mean value. At 50 mbar, the zonally averaged zonal wind fluctuations are significant and amount to almost 1 km s^-1 between its highest and lowest values. We investigate the flow properties during phase II in detail in Sect. 5.1. At time t ~ 270, this quasi-periodic evolution of the wind seems to come to an end and the mean zonal wind decreases. We call this phase III of the flow evolution. As shown by the solid contours in Fig. 6 (top panel), the beginning of this phase coincides with a significant decrease of the Richardson number Ri. The latter is defined by the relation (13)where N² stands for the Brunt-Väisälä frequency and is calculated according to (14)is which the first term inside the square bracket denotes the isentropic gradient. As discussed by Showman & Guillot (2002) in the context of hot-Jupiter atmospheres and later investigated in more details by Li & Goodman (2010), sheared flows in stratified atmospheres are prone to a vertical shear instability when Ri is smaller than a threshold value that depends on the gas thermodynamics but is on the order of 1 / 4. In our simulation, the coincidence between the end of phase II and the appearance of zonally averaged Richardson numbers smaller than 0.1 suggests that such an instability develops at that time and starts perturbing the flow. In Sect. 5.2 we investigate phase III in more detail and give compelling evidence that this is indeed the case.

5.1. Barotropic Kelvin-Helmholtz instability (phase II)

5.1.1. Equatorial jet variability

During phase II of model HighRes, the equatorial jet periodically alternates between two states, during which it is either almost perfectly zonal in the vicinity of the equator or display meanders (Fig. 7). The jet velocity increases when the jet is zonal and slows down when the meander amplitude saturates and decreases. To proceed, we next decomposed the horizontal velocity as the sum of a time-averaged component and a fluctuating part: (15)where the time averaged is performed using 100 dumps evenly spaced between t = 170 and t = 270. In the following, we use the notations (resp. ) and u′ (resp. v′) to denote the x (resp. y) component of these velocities. Using this decomposition, we define the specific kinetic energy of the fluctuations as (16)and the meridional gradient of the total vorticity of the flow (including the planetary scale vorticity associated with Ω_p), (17)%

Taking advantage of the barotropic nature of the flow (Fig. 6, top panel), we now focus on its structure at the 50 mbar pressure level until the end of this section.

Figure 8 (top panel) shows that the time-averaged value of the rms of E_K during phase II reaches values on the order of a few hundred m s^-1 at 50 mbar in the jet core, interior to the region where ξ_y vanishes, and drops rapidly outward of this region. It is a well-known result of hydrodynamics that the existence of two extrema in the flow vorticity profile (such as shown in Fig. 8) is a necessary condition to destabilize a two-dimensional incompressible flow (see, for example, Vallis 2006). Here, the flow is compressible, stratified, and three dimensional, and such a simple criterion does not apply, strictly speaking. Nevertheless, we can still expect that it remains a good guide for linear stability because the deviations from both compressibility and two-dimensionality are small. The rapid drop of E_K outside of the two locations where ξ_y vanishes also argues in favor of this interpretation and strongly suggests that the jet is subject to a Kelvin-Helmholtz instability that originates from the velocity meridional gradient and drives the oscillations seen during phase II.

5.1.2. Linear instability properties

To investigate this possibility in more detail, we next focus on a single oscillation of the jet. It is highlighted as the region shaded in light gray in Fig. 6 (bottom panel). In practice, the analysis that follows was conducted by restarting model HighRes at t = 220 for about seven days, saving the simulation data every 4000 s to sample the fluid evolution at high frequency. Figure 9 shows that the rms fluctuations of v′ at the equator (solid blue curve), initially on the order of 20 m s^-1, grows during about three days until they reach an amplitude of ~350 m s^-1, after which they decay in roughly one day. The spatial structure of these fluctuations is further illustrated in Fig. 10 using two snapshots that illustrate the growing (left panel) and maximum (right panel) phases of the evolution. During the growing part (here plotted at t = 223.3), v′ clearly displays a regular and oscillating pattern at the equator, with a well-defined spatial period of one - fifth of the domain size, suggestive of a linear modal growth. The background jet structure is only weakly modified during this phase. We note that the jet velocity reaches its maximum at this point (black curve in Fig. 9). By time t = 224.7, the velocity fluctuations have grown significantly and show clear signs of saturation, as evidenced by the colliding positive and negative contours. Before examining this saturation phase in more detail in the next paragraph, we first quantify the properties of the fastest growing mode of the instability by performing a spectral analysis of the flow properties. This was done by computing the amplitude of the Fourier coefficients of v′(x,y,t) in the zonal direction. We then spatially averaged these coefficients in the vicinity of the planet equator and show their variations as a function of the normalized zonal wavenumber k_x and time in Fig. 11 (left panel)². From t = 222 to t = 224, the flow evolution is dominated by the k_x = 5 mode. Consistent with the idea of a linear instability, this mode grows exponentially with a typical timescale of about one day (Fig. 11, right panel). It reaches its maximum amplitude at t ~ 224, shortly after higher k_x modes amplitudes also start to grow, presumably as a result of nonlinear interactions with the k_x = 5 mode. This later growth of higher k_x modes is illustrated in the right panel of Fig. 11 for the particular cases of the modes k_x = 7 and k_x = 9. At t = 225, a wide range of spatial scales, with wave number k_x up to 15 (see left panel of Fig. 11), have reached a significant amplitude. For the entire duration of this additional simulation, we also find that larger scale modes (with k_x ranging from 1 to 3) have significant amplitudes. Their evolution, however, is quite different from that of the modes with k_x ≥ 5. They do not display any signature of an exponential growth, and their amplitude is modulated by no more than a factor of two during the simulation. It is likely that these properties reflect the nonlinear feedback of the instability on the large-scale structure of the equatorial jet.

To summarize, we have determined the wave number of the most unstable mode of the instability (k_x = 5) and its growth rate (σ = 1 day^-1 = 3.5 × 10^-6 s^-1). In principle, these properties (as well as the mode phase velocity, see below) could be compared to the result of a numerical linear instability analysis. However, for a 3D compressible flow, such an analysis is tedious and beyond the scope of this paper. Here, we only provide a very crude estimate of these properties based on 2D incompressible flows. A useful and well-studied example of such flows is the so-called Bickney jet, for which the jet profile is given by (18)where the parameters U_B and L_B caraterize the jet. The Bickney jet is known to be unstable to the Kelvin-Helmoltz instability, and the most unstable wavelength satisfies the relation k_xL_B ~ 1 (Drazin & Reid 1981). The core of the jet in our simulations (i.e., within the region where the vorticity is maximum) can be well fit by the Bickney jet meridional profile (see Fig. 8, bottom panel) for which L_B = 2 × 10⁷ m. Given the size of our computational box and the above scaling, this is consistent with the instability having a most unstable wavenumber k_x = 5, as seen in the simulations. However, in this case, the predicted growth rate σ_th is on the order of ϵk_xU_B, with ϵ ~ 0.15. This translates into a growth rate that is more than an order of magnitude faster than measured in the simulation. This large disagreement illustrates the limit of a naive 2D reasoning and suggests that three-dimensional effects, compressibility, or the finite resolution of our simulations likely affect the flow. More work is needed to clarify the properties of the instability, and, more importantly, the conditions under which it develops.

5.1.3. Nonlinear saturation (due to shocks?)

As mentioned above, the background jet becomes significantly distorted when the amplitude of the perturbations reaches large amplitudes (see right panel of Fig. 10). The nonlinear interaction that results quickly dampens the flow velocity fluctuations, the jet slows down (black curve in Fig. 9) and returns to a more zonal structure such as shown in Fig. 7 (top panel). We have found temperature fluctuations of a few hundred Kelvin during this phase (Fig. 12, bottom panel). They tend to be associated with the region of the jet that displays the largest meanders (Fig. 12, top panel). In addition, sharp features in the zonal velocity (such as seen at the location x ~ −0.2 × 10⁸ and y ~ −0.1 × 10⁸, for example) are ubiquitous in snapshots of the flow during the nonlinear stage of the instability.

For these reasons, and because of the supersonic nature of the jet velocity, it is natural to ask whether saturation occurs through shocks. Detecting shocks in 3D non-steady flows is known to be extremely difficult to achieve in a systematic manner. Here, we only give some simple arguments that suggest that shocks play no role in the instability saturation. First, as in Sect. 4.2, we computed the distribution of η given by Eq. (12). We again found that the highest value of η never exceeds a few percent at 50 mbar (and, in fact, at all pressure levels). The sharp gradients mentioned above in the zonal velocity are compensated for by regions with significant vertical downward flow. This suggests that the nonlinear saturation of the instability is not associated with large compressive events.

A second argument is that in the frame where the shock would be stationary, the flow is just about sonic. This can be seen by plotting the equatorial meridional velocity fluctuations in a spacetime diagram in the (x,t) plane, also known as a Hovmöller plot (see Fig. 13). It shows that the eastward velocity of the growing wave pattern amounts to about 6 km s^-1 (see solid black line)³. Because the flow velocity in the frame rotating with the planet reaches at most 8 km s^-1 (see upper panel of Fig. 12), this means that the zonal wind velocity upstream of a putative shock would be lower than 2 km s^-1 in the frame in which the flow structures are stationary. This is similar to (and in fact, slightly lower than) the sound speed, so that the flow Mach number is at best on the order of unity. If shocks exist, they are only weak shocks, with an upstream Mach number on the order of 1 at best. More work is needed to characterize the dynamical mechanism responsible for the saturation and to investigate whether this result holds across the entire parameter space, but the present simulation seems to rule out the presence of shocks, at least for the set of parameters we considered.

5.2. Vertical shear instability (phase III)

The nature of the flow changes during phase III of the simulation. A typical snapshot of the equatorial zonal wind in the (x,P) plane during this phase is plotted in Fig. 14 (left panel, here shown at t = 319). It demonstrates that the zonal wind at the equator displays high-frequency variations in the x-direction that seem to extend over the entire atmosphere.

5.2.1. Identification of the instability

By-eye measurement of the oscillations spatial period (Fig. 14, right panel, see the black horizontal line) gives a typical value of about 8000–10 000 km. By comparison, the most unstable mode of the vertical shear instability has a typical wave number k_x that satisfies the relation k_xH ~ 0.5 (Li & Goodman 2010). For a temperature of about 1800 K such as we used here, the associated wavelength amounts to about 11 000 km, remarkably close to our measurement. As an additional diagnostic, we plot in Fig. 15 a space-time diagram of the high-frequency component of the zonal wind (hereafter noted δu). The latter is calculated by applying a high-pass filter to individual snapshots such as shown in Fig. 14. To do so, we used a Hamming function with a half width of 30 zonal cells as a low-pass filter and subtracted the low-pass-filtered data from the raw data. The comparison with the Richardson number (solid contours in Fig. 15) shows a clear positive correlation: at the bottom of the atmosphere (P ~ 1–10 bar), δu displays a local maximum at pressure depth for which Ri drops below one - fourth. In addition, this maximum is higher when Ri is smaller (see the difference between phases II and III of the simulation, or the increase of δu at t ~ 150 that is coincident with a period of lower Ri values). Taken together, these results are consistent with the analysis of Li & Goodman (2010) and strongly suggest that the vertical shear instability operates at pressure levels where the Richardson number reaches zonally averaged values lower than one - fourth.

5.2.2. Formation of shocks in the upper atmosphere

Figure 15 also shows that δu increases strongly in the atmosphere upper layers, reaching values of up to a few hundred m s^-1 at pressure levels of a few mbar. It is also clear that these velocity fluctuations are connected to the atmospheric activity at 10 bar. Future work is needed to investigate their exact origin, but it is plausible, as recently suggested by Cho et al. (2015) in another context, that they are gravity waves triggered by the vertical shear instability that propagate upward and amplify as the density decreases.

Because these fluctuations reach large amplitudes in the upper atmosphere, we closely investigated whether they eventually steepen into shocks. In agreement with the findings of the preceding sections, the flow remains devoid of such shocks at pressures higher than ~10 mbar, with η being on the order of a few percent at most. However, we found signatures of weak shocks in the top layers of the atmosphere. As an example, we show in Fig. 16 that the zonal wind at P = 1.66 mbar displays strong and high-frequency variations near the equator, accross which the zonal velocity appears to decrease by a few km s^-1. These fluctuations are associated with values of η up to 35% and also display significant increase of the temperature by a few hundred Kelvin (see Fig. 17, left panel). Given the short radiative timescale at this pressure (~10⁴ s), significant temperature fluctuations like this are indications of fast dynamics. We thus focused on one of these structure (indicated as a dashed line in the left panels of Fig. 17) and plotted the profiles of the zonal wind and of the temperature along this line (Fig. 17, right panel). There is a clear discontinuity in both profiles at x ~ 1.72 × 10⁸ m that we identify as a shock. At this location, the zonal velocity rapidly decreases from 7 km s^-1 to 4 km s^-1 within a handful of cells. Simultaneously, the temperature increases from ~800 K to about 1050 K. We measured the shock front velocity V_sh by restarting the simulation for a short time, saving output data at a high frequency. We found for V_sh a value on the order of 3.9 km s^-1 (not shown), which, along with a typical sound speed of 2 km s^-1 at this location, indicates that the shock upstream Mach number is roughly M₁ ~ 1.5. The associated temperature discontinuity shown in Fig. 17 (right panel) can then be compared with the prediction of the Rankine-Hugoniot condition: (19)where T₁ and T₂ are the upstream and downstream gas temperatures, respectively. For the case considered here (M₁ = 1.5 and γ = 1.4), Eq. (19) gives T₂/T₁ = 1.32, in very good agreement with the simulation (whose ratio is about 1.35). This confirms that the numerous structures identified in Fig. 17 (left panel) with the contours of η are indeed shocks. We did not attempt to quantify the frequency, Mach number distribution, and preferred locations of these shocks. In addition to being a difficult task on its own to be carried in a systematic manner, there are also limitations and artifacts of the present setup (see below) that have led us to postpone a detailed characterization of these shocks to future work.

5.2.3. Dissipation in the deep atmosphere and implications for the interior

Figure 15 also shows a significant increase of the temperature at P ~ 10 bar (see the white dashed contours) up to a temperature of about 2400 K (which should be compared to T_eq = 1800 K at that level). Again, there is a strong correlation between this rise and the flow activity, which suggests that the former is a consequence of the latter. This is expected and only a result of small-scale kinetic energy being dissipated and transformed into heat. In RAMSES, such a thermalization is naturally captured because we solve the total energy equation and do not need to use an explicit dissipation coefficient (in this sense, the dissipation is solely numerical in origin). At 10 bar and below, the radiative timescale τ_rad goes to infinity and heated gas cannot cool radiatively anymore: in other words, heat deposited in the inert layer as a result of kinetic energy dissipation accumulates and temperature rises. As shown in Fig. 15, this increase in temperature at 10 bar also increases the vertical temperature gradient (in absolute values) and helps decrease the Richardson number through its dependence on the Brunt-Väisälä frequency. This effect can be quantified by measuring the kinetic energy flux in the vertical direction. To do so, we computed the kinetic energy flux per unit area according to the relation (20)where the integral is calculated at the pressure level P. In both models LowRes and HighRes, Fig. 18 shows that we recover a negative vertical flux of kinetic energy, in quantitative agreement with Showman & Guillot (2002) and with a similar value of about −2000 to −3000 W/m². We note that this agreement is somewhat fortuitous and should not be taken too seriously: in the equatorial β-plane model such as used here, motions gradually go to zero away from the equator, so that F_KE depends on the arbitrary location of the boundaries in the y-direction. However, this qualitative agreement with Showman & Guillot (2002) supports the fact that the existence of this flux is reliable.

Figure 18 also shows that a higher negative flux is present in the deep atmosphere (P ~ 5–10 bar) in model HighRes than in model LowRes, precisely at the location where we see the temperature increase during phase III. This additional flux is likely associated with the vertical shear instability and illustrates the idea of Showman & Guillot (2002) that downward kinetic energy transport can be associated with extra heating in the deep layers of the atmosphere. In addition to the increased kinetic energy flux suggested by Fig. 18, turbulence in a stably stratified atmosphere also transports heat downward (Youdin & Mitchell 2010) and is likely at play here as well.

As first recognized by Guillot & Showman (2002), such downward energy fluxes are similar to, and even greater than, the internal cooling flux of a typical inflated hot Jupiter. This mechanism, that is, conversion of stellar energy into kinetic energy that is transported and dissipated deep in the atmosphere, is thus still a viable explanation of the problem of the large radius of highly irradiated giant planets. The main difficulty, which is faced by most proposed mechanisms, is that this energy must be deposited deep enough inside the planet to significantly affect its thermal evolution (Ginzburg & Sari 2015), possibly in the deep adiabat (Guillot & Showman 2002). Depending on the precise planet and its age, Guillot & Showman (2002) predicted this level to be in the 100–1000 bar range. Although Fig. 18 seems to imply that not much energy is deposited below the 10–20 bar level, it should be clear that in our model this barrier is most probably an artifact of the numerical setup. As discussed above, the presence of an inert layer below 10 bar, where τ_rad becomes infinite, provides a strong positive feedback on the wave activity and the heating at this level. If we had a more realistic, gentler increase of the radiative timescale, the level of maximum energy deposition would probably be pushed deeper. We remark that Fig. 15 displays clear evidence of significant small-scale activity at pressure levels of a few bar, that is, largely above the inert layers, early in the simulation (t< 300 days). This is suggestive evidence that an infinite radiative timescale is not needed to trigger the vertical shear instability. But we cannot yet assess the maximum depth at which energy deposition will occur.

Therefore, while our model illustrates the possibility that vertical shear instability can develop in hot-Jupiter atmospheres, their long-term effect on the flow remains to be determined with a more realistic treatment of the radiative properties of the deep atmosphere. Only then will we be able to properly quantify whether the mechanism first proposed by Guillot & Showman (2002) can realistically account for the inflated radius of strongly irradiated giant planets.

Before ending this section, and with the above caveats in mind, we stress that the effect of the vertical shear instability on the equatorial jet velocity is significant. While its time-averaged value over phase II of the flow evolution is about 7700 m s^-1 at 50 mbar, the equatorial wind decrases to 6700 m s^-1 when averaged between t = 400 and t = 450. We have checked that the spatial distribution of the temperature above 1 bar is almost identical during the two phases (not shown) and cannot be responsible for the jet velocity decrease. It is instead likely that the small-scale velocity fluctuations associated with the vertical shear instability act as a form of drag that slows down the jet, possibly mediated by gravity waves (see Watkins & Cho 2010) and, as shown above, by shocks, both excited by the vertical shear instability. Given the limitations of the numerical setup used here, a detailed and quantitative investigation of this possibility is beyond the scope of this paper, but opens up the possibility of developing a physically motivated subgrid scale model for the dissipation in the atmosphere of hot Jupiters that could be incorporated into traditional GCMs.

6. Discussion and conclusion

6.1. Main results

Using a series of high-resolution idealized simulations of hot-Jupiter atmospheres that solve the Euler equations with a finite-volume shock-capturing scheme, we have found the following results:

• Numerical simulations performed in the framework of theequatorial β-plane model agree well with results published in the literature using a wide range of models and elaborate treatments of the radiative effects.

• A supersonic, equatorial, eastward jet forms quickly in the upper layers of the atmosphere. At P ~ 50 mbar, its zonally averaged velocity reaches about ~7 km s^-1 and is found to be sensitive to the meridional resolution (or, equivalently, to the dissipation).

• At high enough spatial resolution (or low enough dissipation), the jet displays strong velocity fluctuations that can be attributed to meander formations upward of ~1 bar and smaller scale fluctuations with an amplitude of a few tens of m s^-1 at pressure levels P ~ 1–10 bar.

• The meander formations are clearly associated with a barotropic Kelvin-Helmholtz instability that results in quasi-periodic modulations of the jet velocity with a typical period of aboutten days. The properties of the flow are broadly consistent with the expectations of linear stability analysis. Temperature fluctuations of a few hundred Kelvin are found at the photosphere of the planet at the peak of the instability, most likely a result of adiabatic compression associated with the equatorial jet meanderings. Future work is needed to determine whether these variations are compatible with the observed upper limit of 2.7% of the dayside variability of HD 189733b (Agol et al. 2010; Knutson et al. 2012) and whether they could be observable with the JWST using brightness mapping such as discussed by de Wit et al. (2012).

• The smaller scale fluctuations are likely associated with a vertical shear instability. They correlate nicely with the locations where the Richardson number is smaller than one - fourth. They create zonal variations of the jet velocity with a spatial scale of ~10⁴ km that is consistent with the most unstable mode predicted by a linear analysis (Li & Goodman 2010).

• The dissipation of the kinetic energy associated with the vertical shear instability results in a substantial increase of the temperature at P ~ 10 bar. This thus confirms that the atmosphere converts stellar energy into kinetic energy that is transported downward to be deposited at deeper levels (Showman & Guillot 2002). A better treatment of the lower boundary is needed to know whether the deposition level can be deep enough to affect the interior and help explain the radius anomaly (Guillot & Showman 2002).

• We find weak shocks in the upper layers of the atmosphere (P< 10 mbar). They have typical upstream Mach numbers of between 1 and 2 and create temperature fluctuations of a few hundred Kelvin. At higher pressure (P> 10 mbar), we find no shocks despite the supersonic nature of the equatorial jet.

6.2. Limitations and future work

These results should not hide the limitations of the work presented here that are as many avenues for progress. As noted in Sect. 2, the model depends on a number of free parameters. The motivation of the present paper was to choose a unique set of these parameters so that the flow properties match those of the benchmark calculation described by Heng et al. (2011b). Even if the dynamical mechanisms responsible for the instabilities highlighted here are fairly general, future work is needed to systematically investigate the sensitivity of the results presented here to each of these parameters and to make quantitative predictions that would be valuable for more realistic atmospheres. For example, we note that the zonal wind velocity we obtained tends to be higher than commonly found by other authors. Is this because we work in the framework of the equatorial β-plane model? Or is it related to the form of the thermal forcing we use, and in particular to the value of L_th? Another question is that of convergence with spatial resolution: it is possible (and even likely!) that numerical dissipation still affects the instabilities described in Sects. 5.1 and 5.2. This may affect their saturation properties and their quantitave effect of the mean flow, and particularly the velocity of the jet. The present work shows that a proper convergence study (systematically varying both the meridional and the vertical resolution), even though it is very computationally demanding, is clearly needed.

Perhaps the most stringent limitation of the present work comes from the presence of an inert layer below 10 bar where the radiative timescale τ_rad goes to infinity. The strong increase in temperature we see at this location is probably overestimated because of the inability of the gas to cool radiatively. This problem may be somewhat mitigated upon noting that τ_rad is expected to increase rapidly in the deep layers of hot-Jupiter atmospheres (P ≤ 10–100 bar). For example, Showman et al. (2008) use τ_rad = 10⁸ s at 20 bar. This is more than 300 planet days and much faster than the timescale of a few days that is associated with the vertical shear instability (see Fig. 15), so that the effect of choosing an infinite value for τ_rad may not be as severe as naively expected. Nevertheless, it is clear that a strong increase in temperature such as found here acts as a positive feedback onto the vertical shear instability by reducing the Richardson number Ri. Whether it affects the findings presented in this paper, and if it does, by how much, remains to be clarified. One possibility to do so is to replace the Newtonian relaxation scheme with a simplified radiative transfer scheme (Heng et al. 2011a; Rauscher & Menou 2012). This may alleviate the problem in the deep atmosphere.

The need to strengthen these results is made even more important because standard GCM codes that use the primitive formulation of the hydrodynamic equations are unable to account for such vertical shear instabilities. By assuming hydrostatic equilibrium in the vertical direction, they assume de facto that such an instability does not exist. Instead, they rely on subgrid scale modeling to include their effect on the flow. Future work is therefore needed to validate these subgrid models and/or develop more appropriate approaches if necessary. This is an important step to validate the long-term use of such codes as the primary tool to model hot-Jupiter atmospheric flows.

---

Acknowledgments

S.F.’s research is funded by the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC Grant agreement No. 258729. S.F. acknowledges the hospitality of the Kavli Institute of for Theoretical Physics, Santa Barbara, where some of this work was completed during the program “Wave-flow interaction in Geophysics, climate, astrophysics and plasmas”, as well as enlightening discussions with Gwendal Rivière.

References

Appendix A: Equilibrium temperature and cooling timescale

We give here for reference the pressure profiles of and τ_rad that we used here. They are a simplified version of the profiles presented by Heng et al. (2011b). The thermal timescale τ_rad is calculated with the following relation: (A.1)where P₀ = 1 bar, P₁ = 10 bar, τ₀ = 10⁵ and τ₁ = 10^7.5 s. The dimensionless exponent α₀ and α₁ amounts to −0.41 and −2.5, respectively. Likewise, is given as a combination of linear functions of the logarithm of the pressure: (A.2)where P_-2 = 10^-2 bar and T_-2, T₁ and T₂ are set to 1100, 1800, and 3000 K, respectively. The slopes of the linear relations are given by γ_-2 = 100 K, γ₁ = 233 K, and γ₂ = 983 K. The relation between τ_rad and P and between and P are shown in Fig. A.1, where we also display the variation of T_day and T_night such as defined in Sect. 2.2. All curves are meant to be compared with, for example, Fig. 7 of Heng et al. (2011b), with which they agree well.

Appendix B: Numerical tests

B.1: Baroclinic instability in an adiabatic atmosphere

As argued by Polichtchouk et al. (2014), the growth of a baroclinic wave is a severe test for codes that pretend to accurately describe atmospheric flows because it grows slowly from infinitesimal perturbations. This is even more so for finite-volume codes such as RAMSES that have problems handling hydrostatic equilibria, and we found that problem to be very helpful in assessing the reliability of our setup. In this appendix, we thus qualitatively reproduce one of the models presented by Polichtchouk & Cho (2012), namely their equatorial jet case, since the base flow is close to the jet configuration we studied here. As noted by Polichtchouk et al. (2014), the detailed evolution of the atmosphere is very sensitive to the exact structure of the jet and the initial perturbation, so that our goal here is not to reproduce the results of Polichtchouk & Cho (2012) quantitatively, but instead to show that we obtain the same qualitative evolution of the flow. Indeed, since we neither use the same equations (Euler vs. primitive) nor the same geometry (Cartesian vs. spherical), a one-to-one quantitative comparison is not possible.

The structure of the atmosphere at t = 0 is computed assuming thermal wind balance after specifying the jet zonal velocity. The latter is computed using the following relation, adapted to the equatorial β-plane geometry from Polichtchouk & Cho (2012): where G is defined as The previous relations depend on a number of free parameters. In this appendix, we use z_J = 0.8L_z, Δz_J = 0.2L_z, U₀ = 1000 m s^-1, and L_z = 2 × 10⁶ m s^-1. The initial structure of the jet is shown in Fig. B.1 and resembles that shown in their Fig. 4 (right panel) by Polichtchouk & Cho (2012). Next, a localized temperature perturbation of amplitude δT₀ is added at all pressure levels in the atmosphere. It takes the form (B.1)so that it is localized initially at x = 0 and y = L_y/ 4, that is, on the jet northern flank.

We first present and compare the results of two simulations with resolution (N_x,N_y,N_z) = (256,128,32) for which we set δT₀ = 0 (control run) and 1 K (perturbed run). We computed the volume integral of the specific eddy kinetic energy (EKE) as a diagnostic of the baroclinic instability growth: (B.2)As described by Polichtchouk & Cho (2012), the specific EKE grows after a few days in the perturbed case (see Fig. B.2) as a result of a baroclinic wave, but stays very small in the control run during the entire simulation. The specific EKE remains at the level of 10^-10 J/m^-2 (i.e., 17 orders of magnitude smaller than the perturbed model after the growth of the baroclinic wave!), indicating that our numerical scheme accurately conserves the symmetry of the flow and the initial atmospheric hydrostatic equilibrium (the highest meridional velocity at the end of that run at 975 mbar is only on the order of 8 cm s^-1!). In the perturbed run, the specific EKE peaks between t = 30 and t = 35. The flow at 975 mbar shows five well-defined cyclones that roll on both sides of the jet (see Fig. B.3), clearly correlated with temperature fluctuations. This again agrees qualitatively very well with the results of Polichtchouk & Cho (2012). We note that the number of cyclones and the amplitude of the vorticity perturbations they display also agrees quantitatively well with these results. Finally, after a few tens of days, the specific EKE saturates and decays. At the end of the model evolution, the equatorial jet structure is modified and displays a vertical structure that is more barotropic (see Fig. B.1, black contours), in agreement with the finding of Polichtchouk & Cho (2012).

We also performed a resolution study, keeping the number of vertical levels fixed to N_z = 32, and gradually varying the horizontal resolution from (N_x,N_y) = (128,64) to (N_x,N_y) = (512,256), although the highest resolution simulation is only integrated to t = 45. The vorticity distribution at 975 mbar (see Fig. B.4) shows that the minimum resolution required to capture the instability is (N_x,N_y) = (256,128) and that the structure of the flow is qualitatively captured at this resolution, although quantitatively, Fig. B.2 shows that the growth rate has not yet converged. This resolution study gives a useful comparison with the analysis of Polichtchouk & Cho (2012), who used a spectral numerical scheme, and somewhat mitigates their claim that finite difference shemes require resolution higher by an order of magnitude to capture the growth of the baroclinic wave. Here, it appears that a difference of about a factor of two is sufficient.

Overall, the results presented here agree qualitatively well with the results of Polichtchouk & Cho (2012) and Polichtchouk et al. (2014) and give credit to our implementation.

Appendix B.2: A shallow hot-Jupiter model

We next reproduce the benchmark simulation of a shallow hot Jupiter, such as presented by Menou & Rauscher (2009), Heng et al. (2011b), Bending et al. (2013), and Mayne et al. (2014a). As for the case of the baroclinic instability presented above, we adapted the setup described by these authors to the equatorial β-plane. It differs from the deep hot-Jupiter model by the depth of the atmosphere, which extends downward only to 1 bar, and by the parameters entering in the cooling function ℒ. More specifically, τ_rad is a constant equal to half a planet day, and T_eq is given by the function (B.3)where T_perp and β_trop depend on z according to (B.4)and (B.5)The different parameters of the problem are taken to be identical to the original papers mentioned above and are recalled for completeness in Table B.1. In our case, we also need to specify the location of the upper boundary of the domain, which we take to be L_z = 3 × 10⁶ m, while we retain for L_x and L_y the same values as for the deep model described in Sect. 2.4. L_th and β also take values identical to those of the main body of the paper. We started our simulation with an atmosphere initially at rest and in hydrostatic equilibrium with T = T_perp. We integrated the hydrodynamics equations for 300 planet days, of which the last 200 days were used to produce time-averaged fields. The grid resolution was chosen to be (N_x,N_y,Nz) = (64,32,24).

As for the deep hot-Jupiter model described in Sect. 3, we find that a fast equatorial wind develops with typical amplitude of 1000 m s^-1. Here, however, the flow displays significant fluctuations with typical wind velocity fluctuations of a few tens of m s^-1 (not shown). These fluctuations have also been reported by previous authors and probably result from the interaction between the atmosphere and the bottom domain boundary. The time-averaged spatial distributions of the zonally averaged zonal wind and temperature both agree well with the papers mentioned above (see Fig. B.5). Likewise, the temperature distribution at 750 mbar displays the familiar chevron-like shape that has been identified to be the consequence of planetary scale waves by Showman & Polvani (2010, 2011).

The good agreement of our results with previously published simulations of shallow hot-Jupiter atmospheres demonstrates that the equatorial β plane model is appropriate for studying slowly rotating tidally locked gas giant planets.

All Tables

All Figures