Sheardriven instabilities and shocks in the atmospheres of hot Jupiters
^{1} Laboratoire AIM,
CEA/DSMCNRSUniversité Paris 7, Irfu/Service d’Astrophysique,
CEASaclay, 91191 GifsurYvette, France
email: sebastien.fromang@cea.fr
^{2} Institut d’Astrophysique de Paris and
UPMC, CNRS (UMR 7095), 98bis
boulevard Arago, 75014
Paris,
France
^{3} Univ. Bordeaux, LAB, UMR 5804,
33270
Floirac,
France
^{4} CNRS, LAB, UMR 5804,
33270
Floirac,
France
^{5} University of Bern,
Center for Space and Habitability,
Sidlerstrasse 5, 3012
Bern,
Switzerland
Received:
19
October
2015
Accepted:
7
March
2016
Context. General circulation models of the atmosphere of hot Jupiters have shown the existence of a supersonic eastward equatorial jet. These results have been obtained using numerical schemes that filter out vertically propagating sound waves and assume vertical hydrostatic equilibrium, or were acquired with fully compressive codes that use large dissipative coefficients.
Aims. We remove these two limitations and investigate the effects of compressibility on the atmospheric dynamics by solving the standard Euler equations.
Methods. This was done by means of a series of simulations performed in the framework of the equatorial βplane approximation using the finitevolume shockcapturing code RAMSES.
Results. At low resolution, we recover the classical results described in the literature: we find a strong and steady supersonic equatorial jet of a few km s^{1} that displays no signature of shocks. We next show that the jet zonal velocity depends significantly on the grid meridional resolution. When this resolution is fine enough to properly resolve the jet, the latter is subject to a KelvinHelmholtz instability. The jet zonal mean velocity displays regular oscillations with a typical timescale of a few days and a significant amplitude of about 15% of the jet velocity. We also find compelling evidence for the development of a vertical shear instability at pressure levels of a few bars. It seems to be responsible for an increased downward kinetic energy flux that significantly affects the temperature of the deep atmosphere and appears to act as a form of drag on the equatorial jet. This instability also creates velocity fluctuations that propagate upward and steepen into weak shocks at pressure levels of a few mbars.
Conclusions. We conclude that hotJupiter equatorial jets are potentially unstable to both a barotropic KelvinHelmholtz instability and a vertical shear instability. Upon confirmation using more realistic models, these two instabilities could result in significant time variability of the atmospheric winds, may provide a smallscale dissipation mechanism in the flow, and might have consequences for the internal evolution of hot Jupiters.
Key words: hydrodynamics / instabilities / shock waves / methods: numerical / planets and satellites: atmospheres
© ESO, 2016
1. Introduction
In the past decade, hotJupiter observations have evolved from characterizing only the orbital (semimajor 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 temperaturepressure 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 wellunderstood 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; DobbsDixon & 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 hotJupiter 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 semimajor axis of their orbit. This has been given a theoretical basis by the semianalytical 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 threedimensional 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 hotJupiter atmospheres, such as the eastward shift of the hot spot at the photosphere (Knutson et al. 2007) or the daynight temperature contrast measurements (PerezBecker & 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 hotJupiter 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 socalled 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 (DobbsDixon & Lin 2008; DobbsDixon et al. 2010; Mayne et al. 2014a,b). These groups solved the NavierStokes 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, DobbsDixon et al. (2010) reported shocklike 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 hotJupiter 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)^{1}. 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 hotJupiter 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 finitevolume shockcapturing 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, finitevolume 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 hotJupiter 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 (DobbsDixon & 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 wellknown approximation in atmospheric dynamics that has been successfully used recently in the context of hotJupiter 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 hotJupiter 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 lowresolution 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 highresolution simulation in Sect. 5. We show that the flow features variability at different spatial and temporal scales that we relate to wellknown 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 xdirection 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 pressuredependent 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 xdependence is reminiscent of their longitudinal profile, a similar dependence would have little meaning in the ydirection. 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 finitevolume scheme based on the MUSCLHancock Godunov method (Toro 1997). By systematically upwinding all the waves that enter the problem, finitevolume codes are intrinsically stable (provided they satisfy the CourantFriedrichLax, 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 wellknown 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 CrankNicholson 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 MUSCLHancock 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 xdirection. 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 ydirection 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 hotJupiter 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 lowresolution 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.
Fig. 1 Color contours of the timeaveraged zonal mean zonal velocity (left panel) and temperature (right panel) in the (y,P) plane for the lowresolution model. The time average is performed from t = 200 days to t = 500 days using 300 snapshots. 

Open with DEXTER 
Planet parameters:
we followed the shallow hotJupiter model of Heng et al. (2011b) and considered a planet with a radius a_{p} = 10^{8} 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 Ω = Ω_{p}sin(φ), 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 ydependence. 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^{6} 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 hotJupiter atmosphere simulations. This is because the typical pressure scaleheight H is large in this case (typically 10^{3} 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 Earthlike simulation.
3. Lowresolution 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 onetoone 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., volumeintegrated) 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).
Fig. 2 Spacetime diagram (i.e., time evolution of the pressure profile) of the zonally averaged zonal wind at the equator in the lowresolution model. 

Open with DEXTER 
Fig. 3 Timeaveraged temperature (color contour) and horizontal wind (arrows) for the lowresolution model at pressure level P = 1.66 mbar (top panel), P = 97 mbar (middle panel) and P = 4.4 bar (bottom panel). The time average is performed from t = 200 days to t = 500 days using 300 snapshots. 

Open with DEXTER 
The daytonight 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
Model properties.
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 timeaveraged (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 largescale longitudinal structure. In model LowRes, grid cells have a meridional extent of dy = 10^{7} 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
Fig. 4 Left axis: variations of the 50 mbar zonally and timeaveraged zonal wind at the equator with N_{y}. The blue squares correspond to the models with a zonal resolution N_{x} = 128, the red circle shows model LowRes, and the black plusses show the results of the models with a meridional resolution N_{y} = 65 and varying number of cells in the zonal direction. Right axis (green stars): variations of the amplitude of the highfrequency fluctuations of the zonally averaged zonal wind fluctuations (see text for details on its calculation) with N_{y} for the series of models with N_{x} = 128. 

Open with DEXTER 
Fig. 5 Spatial distribution of η at the pressure level P = 1.66 mbar in model 128 × 65 (color contours). The arrows show the amplitude and direction of the hozizontal wind. 

Open with DEXTER 
As discussed in the introduction, it has been a longstanding problem to determine whether shocks exist in hotJupiter 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 shockcapturing 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 finitevolume codes such as RAMSES, such a threshold indeed corresponds to a velocity jump larger than the sound speed across a discontinuity (or, using the RankineHugoniot 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 zonalwind 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 highfrequency 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 lowfrequency 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. Highresolution model
Fig. 6 Top: evolution of the vertical profile of the zonally averaged zonal wind at the equator in a timepressure plane for model HighRes. The horizontal dashed line marks the 50 mbar level. The solid lines are contours of the zonally averaged Richardson number (see text). The contours are for 0.1 (interior contour) and 0.25 (exterior contour). Bottom: time variations of the zonally averaged zonal wind at the equator at the pressure level P = 50 mbar for model HighRes (see text for a discussion of the three phases displayed in both panels). 

Open with DEXTER 
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 quasiperiodic 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 spinup 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 quasisteady 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 quasiperiodic oscillations around a welldefined 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 quasiperiodic 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^{2} stands for the BruntVä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 hotJupiter 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 KelvinHelmholtz instability (phase II)
Fig. 7 Color contours showing the horizontal distribution of the zonal wind at time t = 208 (top panel) and t = 213 (bottom panel) in model HighRes at P = 10 mbar. The arrows display the horizontal wind vectors. Note the clear meandering of the equatorial jet in the bottom panel as opposed to the more zonal structure of the equatorial jet in the top panel. 

Open with DEXTER 
Fig. 8 Top: color contours show the rms of the eddies specific kinetic energy (u^{′2} + v^{′2}) spatial distribution in the horizontal plane at a pressure level P = 50 mbars for model HighRes (during phase II of the evolution), while arrows display the steady component of the wind for the same time interval. The thick solid line plots the locus of the points where the meridional gradient of total vorticity (planet+flow) vanishes. Bottom: yprofile of the timeaveraged zonal velocity (during phase II) at the vicinity of the equator and at pressure level P = 50 mbars for model HighRes (blue crosses) compared with a Bickney jet profile (red line – see text for details). The two vertical arrows mark the location of vanishing meridional gradient of total vorticity. 

Open with DEXTER 
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 timeaveraged 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 timeaveraged 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 wellknown 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 twodimensional 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 twodimensionality 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 KelvinHelmholtz instability that originates from the velocity meridional gradient and drives the oscillations seen during phase II.
Fig. 9 Time evolution of the root mean square of the meridional velocity fluctuations v′ (blue curve – see text for details) and mean jet velocity (black curve) at the equator. The red circles indicate the different time at which the jet structure is plotted in Fig. 10. 

Open with DEXTER 
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 welldefined 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)^{2}. 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 largescale structure of the equatorial jet.
Fig. 10 Color contours showing the horizontal spatial distribution of the zonal wind at P = 50 mbars at times indicated with red dots in Fig. 9, namely t = 223.3 (left panel) and t = 224.7 (right panel). Black contours plots the meridional velocity fluctuations v′. Contours are shown every 50 m s^{1} from −100 to 100 m s^{1} on the left hand side panel and every 100 m s^{1} from −400 to 400 m s^{1} on the right hand side panel. In both panels, positive (resp. negative) contours are shown with solid (resp. dashed) lines and the zero contour is omitted. 

Open with DEXTER 
Fig. 11 Left: amplitude of the meridional velocity fluctuations Fourier transform in the (k_{x}time) plane for model HighRes, spatially averaged over the region y < 5 × 10^{7} m. Contours are for , , and . The region surrounding the maximum is filled in gray. Right: same as the left panel, but showing the time evolution of particular modes: k_{x} = 5 (thick solid line), k_{x} = 7 (dashed line,) and k_{x} = 9 (dotdashed line). The straight dashed line represents an exponential growth with a timescale of one day. 

Open with DEXTER 
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 wellstudied example of such flows is the socalled 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 KelvinHelmoltz instability, and the most unstable wavelength satisfies the relation k_{x}L_{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^{7} 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_{x}U_{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 threedimensional 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^{8} and y ~ −0.1 × 10^{8}, 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 nonsteady 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.
Fig. 12 Zonal velocity (top panel) and temperature (bottom panel) at time t = 224.7 in model HighRes (i.e., when velocity fluctuations reach their highest value) at pressure level P = 50 mbar. Jet meanderings upstream of the substellar point, located at (x,y) = (0,0), are clearly associated with temperature fluctuations. 

Open with DEXTER 
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)^{3}. 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 highfrequency variations in the xdirection that seem to extend over the entire atmosphere.
5.2.1. Identification of the instability
Byeye 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_{x}H ~ 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 spacetime diagram of the highfrequency component of the zonal wind (hereafter noted δu). The latter is calculated by applying a highpass 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 lowpass filter and subtracted the lowpassfiltered 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.
Fig. 13 Hovmöller spacetime diagram of the meridional velocity fluctuations v′ at the equator in model HighRes at pressure level P = 50 mbar. The thin black lines have a slope of 6 km s^{1} . 

Open with DEXTER 
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 highfrequency 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^{4} 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^{8} 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} ~ 1.5. The associated temperature discontinuity shown in Fig. 17 (right panel) can then be compared with the prediction of the RankineHugoniot condition: (19)where T_{1} and T_{2} are the upstream and downstream gas temperatures, respectively. For the case considered here (M_{1} = 1.5 and γ = 1.4), Eq. (19) gives T_{2}/T_{1} = 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.
Fig. 14 Zonal velocity at the equator in a (x,P) plane at time t = 319 in model HighRes (left panel). Note the highfrequency variations of the velocity (particularly easy to see in the atmosphere upper layers), superposed on the largescale zonal and vertical variations of the equatorial jet velocity. The right panel shows an enlargement of the dashed box shown in the left panel. The horizontal bar has a length of 10 000 km. 

Open with DEXTER 
Fig. 15 Spacetime (timepressure) diagram of the zonally averaged highfrequency component of the zonal wind (see text for details) in color contours (note that the color table is saturated at 50 m s^{1}). Contours of Ri (solid lines) and T (dashed lines) are overplotted. For Ri, contours are for Ri = 0.1 and 0.25. For the temperature, contours are for T = 1900, 2000, 2100, and 2200 K. 

Open with DEXTER 
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 smallscale 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 BruntVä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^{2}. 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 ydirection. 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 smallscale 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 hotJupiter atmospheres, their longterm 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.
Fig. 16 Zonal velocity in model HighRes at t = 335 at the pressure level P = 1.66 mbar. Arrows indicate velocity vectors and the square marks the region studied in Fig. 17. 

Open with DEXTER 
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 timeaveraged 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 smallscale 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
Fig. 17 Left: color contours of the zonal wind (top panel) and temperature (bottom panel) at time t = 335 in the rectangle box depicted in Fig. 16. The contours show the distribution of η, as given by Eq. (12). Levels are for η = 0.1 and 0.2. Right: profiles of the zonal velocity (black curve) and temperature (red curve) along the dashed line shown in the left panels. For both curves, the empty squares mark the locations of the cell centers. 

Open with DEXTER 
Using a series of highresolution idealized simulations of hotJupiter atmospheres that solve the Euler equations with a finitevolume shockcapturing 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 KelvinHelmholtz instability that results in quasiperiodic 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^{4} 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.
Fig. 18 Vertical profile of the kinetic energy flux during phase III of model HighRes (green curve) and in model LowRes (blue curve), calculated according to Eq. (20). The dashed line marks the location of the top of the inert layer at 10 bar. Both curves are time averaged between t = 350 and t = 450. 

Open with DEXTER 
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 hotJupiter atmospheres (P ≤ 10–100 bar). For example, Showman et al. (2008) use τ_{rad} = 10^{8} 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 longterm use of such codes as the primary tool to model hotJupiter 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 “Waveflow interaction in Geophysics, climate, astrophysics and plasmas”, as well as enlightening discussions with Gwendal Rivière.
References
 Agol, E., Cowan, N. B., Knutson, H. A., et al. 2010, ApJ, 721, 1861 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., Chabrier, G., & Barman, T. 2010, Rep. Prog. Phys., 73, 016901 [NASA ADS] [CrossRef] [Google Scholar]
 Bending, V. L., Lewis, S. R., & Kolb, U. 2013, MNRAS, 428, 2874 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J. Y.K., Menou, K., Hansen, B. M. S., & Seager, S. 2008, ApJ, 675, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J. Y.K., Polichtchouk, I., & Thrastarson, H. T. 2015, MNRAS, submitted [arXiv:1508.05719] [Google Scholar]
 de Wit, J., Gillon, M., Demory, B.O., & Seager, S. 2012, A&A, 548, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DobbsDixon, I., & Lin, D. N. C. 2008, ApJ, 673, 513 [NASA ADS] [CrossRef] [Google Scholar]
 DobbsDixon, I., Cumming, A., & Lin, D. N. C. 2010, ApJ, 710, 1395 [NASA ADS] [CrossRef] [Google Scholar]
 Drazin, P. G., & Reid, W. H. 1981, NASA STI/Recon Technical Report A, 82 [Google Scholar]
 Ginzburg, S., & Sari, R. 2015, ApJ, 803, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J. 2009, ApJ, 693, 1645 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., & Showman, A. P. 2002, A&A, 385, 156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heng, K. 2012, ApJ, 761, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., & Workman, J. 2014, ApJS, 213, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., & Showman, A. P. 2015, Ann. Rev. Earth Planet. Sci., 43, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Frierson, D. M. W., & Phillipps, P. J. 2011a, MNRAS, 418, 2669 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Menou, K., & Phillipps, P. J. 2011b, MNRAS, 413, 2380 [NASA ADS] [CrossRef] [Google Scholar]
 Käppeli, R., & Mishra, S. 2014, J. Comput. Phys., 259, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Knutson, H. A., Lewis, N., Fortney, J. J., et al. 2012, ApJ, 754, 22 [NASA ADS] [CrossRef] [Google Scholar]
 LeVeque, R. J. 1998, J. Comput. Phys., 146, 346 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Li, J., & Goodman, J. 2010, ApJ, 725, 1146 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, B., & Showman, A. P. 2013, ApJ, 770, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Madhusudhan, N., Knutson, H., Fortney, J. J., & Barman, T. 2014, Protostars and Planets VI, 739 [Google Scholar]
 Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014a, A&A, 561, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014b, Geoscientific Model Development, 7, 3059 [NASA ADS] [CrossRef] [Google Scholar]
 Menou, K., & Rauscher, E. 2009, ApJ, 700, 887 [NASA ADS] [CrossRef] [Google Scholar]
 Menou, K., Cho, J. Y.K., Seager, S., & Hansen, B. M. S. 2003, ApJ, 587, L113 [NASA ADS] [CrossRef] [Google Scholar]
 PerezBecker, D., & Showman, A. P. 2013, ApJ, 776, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Polichtchouk, I., & Cho, J. Y.K. 2012, MNRAS, 424, 1307 [NASA ADS] [CrossRef] [Google Scholar]
 Polichtchouk, I., Cho, J. Y.K., Watkins, C., et al. 2014, Icarus, 229, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Rauscher, E., & Menou, K. 2010, ApJ, 714, 1334 [NASA ADS] [CrossRef] [Google Scholar]
 Rauscher, E., & Menou, K. 2012, ApJ, 750, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2010, Geophys. Res. Lett., 37, 18811 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thrastarson, H. T., & Cho, J. Y. 2010, ApJ, 716, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. 1997, Riemann solvers and numerical methods for fluid dynamics (Springer) [Google Scholar]
 Tsai, S.M., DobbsDixon, I., & Gu, P.G. 2014, ApJ, 793, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Vallis, G. K. 2006, in Atmospheric and Oceanic Fluid Dynamics (Cambridge University Press), 745 [Google Scholar]
 Watkins, C., & Cho, J. Y.K. 2010, ApJ, 714, 904 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Mitchell, J. L. 2010, ApJ, 721, 1113 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Stone, J. M., & Rafikov, R. R. 2013, ApJ, 768, 143 [NASA ADS] [CrossRef] [Google Scholar]
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_{0} = 1 bar, P_{1} = 10 bar, τ_{0} = 10^{5} and τ_{1} = 10^{7.5} s. The dimensionless exponent α_{0} and α_{1} 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_{1} and T_{2} are set to 1100, 1800, and 3000 K, respectively. The slopes of the linear relations are given by γ_{2} = 100 K, γ_{1} = 233 K, and γ_{2} = 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.
Fig. A.1 Pressure variations of the equilibrium temperature T_{eq} (top panel) and the Newtonian relaxation time τ_{rad} (bottom panel) used to calculate the cooling function ℒ in the deep hotJupiter model. 

Open with DEXTER 
Appendix B: Numerical tests
B.1: Baroclinic instability in an adiabatic atmosphere
Fig. B.1 Initial (i.e., at t = 0) zonally averaged zonal wind used for the baroclinic instability simulations (colors and white contours). The black contour show the zonally averaged zonal wind at the end of the simulation for the model with resolution (N_{x},N_{y},N_{z}) = (256,128,32). Contours are drawn every 100 m s^{1} from 100 to 900 m s^{1}. 

Open with DEXTER 
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 finitevolume 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 onetoone 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_{0} = 1000 m s^{1}, and L_{z} = 2 × 10^{6} 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_{0} 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.
Fig. B.2 Time evolution of the volumeintegrated eddy kinetic energy (per unit area) in the baroclinic instability test. The solid and dotted lines both share the resolution (N_{x},N_{y},N_{z}) = (256,128,32) and correspond to the model with and without an initial temperature perturbation, respectively. The dashed and dotdashed lines shows the results of the models with resolution (N_{x},N_{y},N_{z}) = (128,64,32) and (N_{x},N_{y},N_{z}) = (512,256,32), respectively. 

Open with DEXTER 
Fig. B.3 Relative vorticity (top panel) and temperature (bottom panel) distribution in the horizontal plane at 930 mbar at time t = 32 showing the development of a baroclinic wave that grows on top of a zonal equatorial jet. The parameters are chosen after Polichtchouk & Cho (2012) and the resolution is (N_{x},N_{y},N_{z}) = (256,128,32). 

Open with DEXTER 
Fig. B.4 Relative vorticity at 930 mbar at time t = 32 for the baroclinic instability growth with resolution (N_{x},N_{y}) = (128,64) (top panel) and (N_{x},N_{y}) = (512,128) (bottom panel). 

Open with DEXTER 
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} = 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 welldefined 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 hotJupiter model
Parameters used for the shallow hotJupiter 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 hotJupiter 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^{6} 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 timeaveraged fields. The grid resolution was chosen to be (N_{x},N_{y},Nz) = (64,32,24).
As for the deep hotJupiter 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 timeaveraged 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 chevronlike shape that has been identified to be the consequence of planetary scale waves by Showman & Polvani (2010, 2011).
Fig. B.5 Zonally averaged zonal wind (top panel) and temperature (bottom panel) in the shallow hotJupiter model. The raw simulation data are averaged over 200 planet days starting at t = 100. See text and Table B.1 for the model parameters. 

Open with DEXTER 
Fig. B.6 Temperature (color contours) and wind velocities at the 750 mbar pressure level in the shallow hotJupiter model, averaged in time over 200 planet days starting at t = 100. 

Open with DEXTER 
The good agreement of our results with previously published simulations of shallow hotJupiter atmospheres demonstrates that the equatorial β plane model is appropriate for studying slowly rotating tidally locked gas giant planets.
All Tables
All Figures
Fig. 1 Color contours of the timeaveraged zonal mean zonal velocity (left panel) and temperature (right panel) in the (y,P) plane for the lowresolution model. The time average is performed from t = 200 days to t = 500 days using 300 snapshots. 

Open with DEXTER  
In the text 
Fig. 2 Spacetime diagram (i.e., time evolution of the pressure profile) of the zonally averaged zonal wind at the equator in the lowresolution model. 

Open with DEXTER  
In the text 
Fig. 3 Timeaveraged temperature (color contour) and horizontal wind (arrows) for the lowresolution model at pressure level P = 1.66 mbar (top panel), P = 97 mbar (middle panel) and P = 4.4 bar (bottom panel). The time average is performed from t = 200 days to t = 500 days using 300 snapshots. 

Open with DEXTER  
In the text 
Fig. 4 Left axis: variations of the 50 mbar zonally and timeaveraged zonal wind at the equator with N_{y}. The blue squares correspond to the models with a zonal resolution N_{x} = 128, the red circle shows model LowRes, and the black plusses show the results of the models with a meridional resolution N_{y} = 65 and varying number of cells in the zonal direction. Right axis (green stars): variations of the amplitude of the highfrequency fluctuations of the zonally averaged zonal wind fluctuations (see text for details on its calculation) with N_{y} for the series of models with N_{x} = 128. 

Open with DEXTER  
In the text 
Fig. 5 Spatial distribution of η at the pressure level P = 1.66 mbar in model 128 × 65 (color contours). The arrows show the amplitude and direction of the hozizontal wind. 

Open with DEXTER  
In the text 
Fig. 6 Top: evolution of the vertical profile of the zonally averaged zonal wind at the equator in a timepressure plane for model HighRes. The horizontal dashed line marks the 50 mbar level. The solid lines are contours of the zonally averaged Richardson number (see text). The contours are for 0.1 (interior contour) and 0.25 (exterior contour). Bottom: time variations of the zonally averaged zonal wind at the equator at the pressure level P = 50 mbar for model HighRes (see text for a discussion of the three phases displayed in both panels). 

Open with DEXTER  
In the text 
Fig. 7 Color contours showing the horizontal distribution of the zonal wind at time t = 208 (top panel) and t = 213 (bottom panel) in model HighRes at P = 10 mbar. The arrows display the horizontal wind vectors. Note the clear meandering of the equatorial jet in the bottom panel as opposed to the more zonal structure of the equatorial jet in the top panel. 

Open with DEXTER  
In the text 
Fig. 8 Top: color contours show the rms of the eddies specific kinetic energy (u^{′2} + v^{′2}) spatial distribution in the horizontal plane at a pressure level P = 50 mbars for model HighRes (during phase II of the evolution), while arrows display the steady component of the wind for the same time interval. The thick solid line plots the locus of the points where the meridional gradient of total vorticity (planet+flow) vanishes. Bottom: yprofile of the timeaveraged zonal velocity (during phase II) at the vicinity of the equator and at pressure level P = 50 mbars for model HighRes (blue crosses) compared with a Bickney jet profile (red line – see text for details). The two vertical arrows mark the location of vanishing meridional gradient of total vorticity. 

Open with DEXTER  
In the text 
Fig. 9 Time evolution of the root mean square of the meridional velocity fluctuations v′ (blue curve – see text for details) and mean jet velocity (black curve) at the equator. The red circles indicate the different time at which the jet structure is plotted in Fig. 10. 

Open with DEXTER  
In the text 
Fig. 10 Color contours showing the horizontal spatial distribution of the zonal wind at P = 50 mbars at times indicated with red dots in Fig. 9, namely t = 223.3 (left panel) and t = 224.7 (right panel). Black contours plots the meridional velocity fluctuations v′. Contours are shown every 50 m s^{1} from −100 to 100 m s^{1} on the left hand side panel and every 100 m s^{1} from −400 to 400 m s^{1} on the right hand side panel. In both panels, positive (resp. negative) contours are shown with solid (resp. dashed) lines and the zero contour is omitted. 

Open with DEXTER  
In the text 
Fig. 11 Left: amplitude of the meridional velocity fluctuations Fourier transform in the (k_{x}time) plane for model HighRes, spatially averaged over the region y < 5 × 10^{7} m. Contours are for , , and . The region surrounding the maximum is filled in gray. Right: same as the left panel, but showing the time evolution of particular modes: k_{x} = 5 (thick solid line), k_{x} = 7 (dashed line,) and k_{x} = 9 (dotdashed line). The straight dashed line represents an exponential growth with a timescale of one day. 

Open with DEXTER  
In the text 
Fig. 12 Zonal velocity (top panel) and temperature (bottom panel) at time t = 224.7 in model HighRes (i.e., when velocity fluctuations reach their highest value) at pressure level P = 50 mbar. Jet meanderings upstream of the substellar point, located at (x,y) = (0,0), are clearly associated with temperature fluctuations. 

Open with DEXTER  
In the text 
Fig. 13 Hovmöller spacetime diagram of the meridional velocity fluctuations v′ at the equator in model HighRes at pressure level P = 50 mbar. The thin black lines have a slope of 6 km s^{1} . 

Open with DEXTER  
In the text 
Fig. 14 Zonal velocity at the equator in a (x,P) plane at time t = 319 in model HighRes (left panel). Note the highfrequency variations of the velocity (particularly easy to see in the atmosphere upper layers), superposed on the largescale zonal and vertical variations of the equatorial jet velocity. The right panel shows an enlargement of the dashed box shown in the left panel. The horizontal bar has a length of 10 000 km. 

Open with DEXTER  
In the text 
Fig. 15 Spacetime (timepressure) diagram of the zonally averaged highfrequency component of the zonal wind (see text for details) in color contours (note that the color table is saturated at 50 m s^{1}). Contours of Ri (solid lines) and T (dashed lines) are overplotted. For Ri, contours are for Ri = 0.1 and 0.25. For the temperature, contours are for T = 1900, 2000, 2100, and 2200 K. 

Open with DEXTER  
In the text 
Fig. 16 Zonal velocity in model HighRes at t = 335 at the pressure level P = 1.66 mbar. Arrows indicate velocity vectors and the square marks the region studied in Fig. 17. 

Open with DEXTER  
In the text 
Fig. 17 Left: color contours of the zonal wind (top panel) and temperature (bottom panel) at time t = 335 in the rectangle box depicted in Fig. 16. The contours show the distribution of η, as given by Eq. (12). Levels are for η = 0.1 and 0.2. Right: profiles of the zonal velocity (black curve) and temperature (red curve) along the dashed line shown in the left panels. For both curves, the empty squares mark the locations of the cell centers. 

Open with DEXTER  
In the text 
Fig. 18 Vertical profile of the kinetic energy flux during phase III of model HighRes (green curve) and in model LowRes (blue curve), calculated according to Eq. (20). The dashed line marks the location of the top of the inert layer at 10 bar. Both curves are time averaged between t = 350 and t = 450. 

Open with DEXTER  
In the text 
Fig. A.1 Pressure variations of the equilibrium temperature T_{eq} (top panel) and the Newtonian relaxation time τ_{rad} (bottom panel) used to calculate the cooling function ℒ in the deep hotJupiter model. 

Open with DEXTER  
In the text 
Fig. B.1 Initial (i.e., at t = 0) zonally averaged zonal wind used for the baroclinic instability simulations (colors and white contours). The black contour show the zonally averaged zonal wind at the end of the simulation for the model with resolution (N_{x},N_{y},N_{z}) = (256,128,32). Contours are drawn every 100 m s^{1} from 100 to 900 m s^{1}. 

Open with DEXTER  
In the text 
Fig. B.2 Time evolution of the volumeintegrated eddy kinetic energy (per unit area) in the baroclinic instability test. The solid and dotted lines both share the resolution (N_{x},N_{y},N_{z}) = (256,128,32) and correspond to the model with and without an initial temperature perturbation, respectively. The dashed and dotdashed lines shows the results of the models with resolution (N_{x},N_{y},N_{z}) = (128,64,32) and (N_{x},N_{y},N_{z}) = (512,256,32), respectively. 

Open with DEXTER  
In the text 
Fig. B.3 Relative vorticity (top panel) and temperature (bottom panel) distribution in the horizontal plane at 930 mbar at time t = 32 showing the development of a baroclinic wave that grows on top of a zonal equatorial jet. The parameters are chosen after Polichtchouk & Cho (2012) and the resolution is (N_{x},N_{y},N_{z}) = (256,128,32). 

Open with DEXTER  
In the text 
Fig. B.4 Relative vorticity at 930 mbar at time t = 32 for the baroclinic instability growth with resolution (N_{x},N_{y}) = (128,64) (top panel) and (N_{x},N_{y}) = (512,128) (bottom panel). 

Open with DEXTER  
In the text 
Fig. B.5 Zonally averaged zonal wind (top panel) and temperature (bottom panel) in the shallow hotJupiter model. The raw simulation data are averaged over 200 planet days starting at t = 100. See text and Table B.1 for the model parameters. 

Open with DEXTER  
In the text 
Fig. B.6 Temperature (color contours) and wind velocities at the 750 mbar pressure level in the shallow hotJupiter model, averaged in time over 200 planet days starting at t = 100. 

Open with DEXTER  
In the text 