Simulating Jupiter’s weather layer. Part I: Jet spin-up in a dry atmosphere

We investigate the dynamics of Jupiter’s upper troposphere and lower stratosphere using a General Circulation Model that includes two-stream radiation and optional heating from below. Based on the MITgcm dynamical core, this is a new generation of Oxford's Jupiter General Circulation Model [Zuchowski, L.C. et al., 2009. Plan. Space Sci., 57, 1525–1537, doi:10.1016/j.pss.2009.05.008]. We simulate Jupiter’s atmosphere at up to 0.7° horizontal resolution with 33 vertical levels down to a pressure of 18 bar, in configurations with and without a 5.7 W m 2 interior heat flux. Simulations ran for 130000–150000 d to allow the deep atmosphere to come into radiative equilibrium. Baroclinic instability generates alternating, eddy-driven, midlatitude jets in both cases. With interior heating the zonal jets migrate towards the equator and become barotropically unstable. This generates Rossby waves that radiate away from the equator, depositing westerly momentum there via eddy angular momentum flux convergence and spinning up a super-rotating 20 m s 1 equatorial jet throughout the troposphere. There are 30–35 zonal jets with latitudinal separation comparable with the real planet, and there is strong eddy activity throughout. Without interior heating the jets do not migrate and a divergent eddy angular momentum flux at the equator spins up a broad, 50 m s 1 sub-rotating equatorial jet with weak eddy activity at


Introduction
Despite more than 40 years of investigation by space-borne instrumentation, the nature of the atmospheric circulation of the Solar System's giant planets remains one of the most enigmatic problems in planetary atmospheric science. The ubiquitous cloud cover effectively prevents remote sensing more than a few tens of km below the tropopause except at radio wavelengths, yet the observable tropospheric meteorology -in the form of intense eastward and westward jets, travelling and near-stationary waves, and long-lived coherent vortices (Ingersoll et al., 2004;Vasavada and Showman, 2005) -appears to penetrate (and perhaps also to originate from) deep below the visible cloud tops.
Jupiter's composition implies that its H-He atmosphere must remain continuously fluid to great depths, but until very recently it was not even possible to determine whether the observed large-scale meteorological circulations penetrate deep into the planet (Busse, 1976) or are relatively shallow, as observed on Uranus and Neptune (Kaspi et al., 2013). This conundrum has apparently been resolved by Kaspi et al. (2018) using measurements from the Juno spacecraft, in which the jet depth is inferred to be 2000-3000 km, although the uniqueness of this interpretation of the Juno measurements has been challenged by Kong et al. (2018). Such considerations are further complicated by the demonstration (Lian and Showman, 2008) that relatively shallow energetic processes, confined to the 'weather layer' around and above the clouds of H O, its detailed application to the gas giant planets is uncertain. In particular, it is unclear whether such a cascade involves spectrally local energy transfer between adjacent wavenumbers, or direct transfer to the zonal jets.
Of particular interest is Jupiter's prograde equatorial jet. Equatorial super-rotation requires a net transfer of eastward angular momentum (AM) into the equatorial region, i.e. a convergence of angular momentum fluxes there. This must be due to eddy processes because in an axisymmetric flow on a sphere eastward AM flux into low latitudes is impossible (Hide, 1969). While the equatorial flow is retrograde at the visible cloud level on Uranus and Neptune, there is super-rotation (prograde flow) on Jupiter and Saturn. Furthermore, the observed zonal flow on Jupiter appears to be stable, but does not satisfy the barotropic stability criterion that u y / 2 2 remain the same sign everywhere (Read et al., 2006, Fig. 1).
The source(s) of energy responsible for the small-scale forcing required to drive the planetary-scale zonal jets on the giant planets has been a rapidly-advancing field in recent years. Over the last several years a number of groups have used numerical models of varying sophistication to explore the processes relevant to the generation of the planetary-scale zonal jets and large eddies in these planets' atmospheres.
In simplified domains, Showman et al. (2006) examined zonal jet formation using the linearised primitive equations on an f-plane, generating deep zonal jets from shallow forcing, and Sayanagi et al. (2008) were able to produce vertically coherent jets in a β-plane simulation with small-scale turbulent forcing using the EPIC model (Dowling et al., 1998). Other simplified models include Li et al. (2006) and Showman (2007), who incorporated simple parametrisations of moist convection into quasi-geostrophic and shallow-water models respectively, under conditions relevant to Jupiter and Saturn, and were able to generate zonal jets. Warneford and Dellar (2014) applied the thermal shallowwater equations to the Jupiter case and, unlike in the standard shallowwater system, were able to generate a strongly prograde equatorial jet. Most recently, Thomson and McIntyre (2016) used a 1 1/2-layer model and a stochastic representation of moist convection to reproduce a statistically-steady weather layer without resorting to large-scale friction.
Using more sophisticated General Circulation Models (GCMs), various groups have investigated the relative importance of effects such as interior heat fluxes, differential radiative forcing, and latent heat release on the formation of zonal jets (Williams, 2003b;Lian and Showman, 2008;Schneider and Liu, 2009;Liu and Schneider, 2010;Lian and Showman, 2010;Liu and Schneider, 2011;Medvedev et al., 2013). Various combinations of these effects have been able to produce flow fields qualitatively similar to observations in some cases. Others have attributed the observed zonal flow to a surface manifestation of deep convective motion (Heimpel et al., 2005(Heimpel et al., , 2015Kaspi et al., 2009).
Work on equatorial super-rotation has shown that it can occur due to the relative importance of generation and dissipation of Rossby waves at different latitudes (Held and Hoskins, 1985;Saravanan, 1993;Schneider and Liu, 2009;Mitchell and Vallis, 2010;Laraia and Schneider, 2015;Polichtchouk and Cho, 2016). Rossby waves deposit prograde AM in their source region (seen as convergence of eddy AM fluxes) and retrograde AM where they dissipate (divergence of eddy AM fluxes). Equatorial super-rotation can occur when the AM flux convergence from waves generated at the equator exceeds the AM flux divergence near the equator generated by dissipation of equatoriallypropagating Rossby waves from off-equatorial latitudes. On Jupiter, there is some evidence for Rossby waves in both the north and south branches of the equatorial jet, the most clear being wavenumber 11-13 waves embedded in Jupiter's equatorial plumes around 8°N, with phase speeds of 25-75 m s −1 and equivalent depth 2 km (Allison, 1990).
Large equator-to-pole temperature gradients (Polichtchouk and Cho, 2016) and fast-rotating planets  act to suppress equatorial super-rotation by increasing extra-tropical baroclinicity and hence the propagation of westward momentum into low latitudes by Rossby waves. The giant planets are typically fast-rotating but with weak equator-to-pole temperature gradients, and so there is a trade-off between these two effects. Observations of Jupiter show horizontal thermal contrasts (along constant pressure surfaces, so also of potential temperature) of order 5 K around 500 mbar; see, for example, latitudeheight maps of temperature derived from Voyager (Conrath et al., 1998) or Cassini (Flasar et al., 2004;Read et al., 2006). This can be compared with vertical potential temperature variations of order a few 10 K between 500-150 mbar , which means potential temperature surfaces in the upper troposphere are close to flat.
Here we investigate the spin-up of Jupiter's jets using a weatherlayer General Circulation Model to simulate Jupiter's troposphere and lower stratosphere. This GCM has been developed over the last 15-20 years based on the dynamical core of the UK Met-Office External Unified Model (ExtUM) (Yamazaki et al., 2004;Yamazaki and Read, 2006). Zuchowski et al. (2009b,c) used it to develop a simple cloud parametrisation over a limited area of the southern hemisphere, and a single-column moist convection parametrisation (Zuchowski et al., 2009a).
This work presents a major upgrade of this model, re-coding it for use in a global configuration, with updated radiation and other parametrisations. The new generation of this model uses the MITgcm (Marshall et al., 1997a,b), which has also been used for giant planet modelling by Lian and Showman (2008) (Jupiter dynamics), Kaspi et al. (2009) (deep convective model of Jupiter), Showman et al. (2009) and Polichtchouk et al. (2014) (extra-solar planets).
The aims of this paper are twofold. First, to present the new generation of our Jupiter model, which has undergone several major changes since it was last used by Zuchowski et al. (2009a). Second, we use the model to study the spin-up of flow in Jupiter's atmosphere, particularly the presence or absence of equatorial super-rotation. In Young et al. (2018b, called Part II hereafter), we study the global distribution of cloud condensates using the new generation of the model.
The new generation of the Jupiter model is described in Section 2, followed by a description of the simulations in Section 3. The atmospheric dynamics in these experiments are analysed in Section 4, followed by a particular focus on equatorial jet spin-up in Section 5, and we conclude in Section 6. Appendix A contains full details of our Jupiter model parametrisations.

Model description
The current generation of our Jupiter model is based on the MITgcm (Marshall et al., 1997a;1997b;Adcroft et al., 2004), a highly-customisable oceanographic and atmospheric GCM with a large user base. It is a major update of the Jupiter version of the Oxford Planetary Unified (Model) System (OPUS), as used by Zuchowski et al. (2009a), which was a limited area model based on the dynamical core of the Met-Office External Unified Model (ExtUM) v4.5. Our initial global version of the Jupiter model used the ExtUM v8.2, but we had problems with gridscale noise near the poles that we could not eliminate using the polar filter or hyper-diffusion. Therefore we eventually rewrote the model using a new dynamical core, and the MITgcm was a sensible choice as other groups had successfully used it to simulate Jupiter (e.g. Kaspi et al., 2009;Lian and Showman, 2010). All the existing physical parametrisations for Jupiter described in earlier work were ported across to the new model.

MITgcm dynamical core
The MITgcm is a finite volume GCM that solves the primitive equations for an incompressible fluid on a rotating sphere. We use the hydrostatic version, without vertical Coriolis terms. It is primarily an oceanographic model, but the vertical coordinate can be swapped from height z to pressure p to produce an atmospheric model. In atmosphere mode it solves the following equations (Adcroft et al., 2018, Sects 1.3.6.2 and 1.4.1.2) h is the horizontal velocity on pressure surfaces, where u and v are the zonal and meridional velocities respectively, ∇ p is the gradient operator in pressure coordinates, = Dp Dt / is the vertical velocity in pressure coordinates (positive downwards), = f 2 sin is the Coriolis parameter, = gz is the geopotential, F includes all additional mechanical forcing, = 1/ is the specific volume, p s is the pressure at the bottom of the domain, c p is the specific heat capacity of dry air at constant pressure, = p p ( / ) 0 is the Exner function, where p 0 is a reference pressure and = R c / , d p with R d the specific gas constant for dry air, θ is potential temperature, and Q is a heating rate per unit mass due to external forcing. Absolute temperature is defined by . Primed quantities are 2D perturbations from 1D reference profiles (Adcroft et al., 2018, Sect. 1.4.1.2).
Eqs. 1 and 2 contain the horizontal and vertical (hydrostatic) components of the Navier-Stokes equation, Eq. 3 is the continuity equation, Eq. 4 describes the evolution of the pressure at the bottom of the domain, Eq. 5 is the ideal gas law, and Eq. 6 is the thermodynamic energy equation.
In our setup, the model solves these equations globally as a function of longitude λ and latitude ϕ, discretised on an Arakawa-C grid. We used the vector-invariant form of the momentum equation, centred second order advection-diffusion for potential temperature, with no harmonic or bi-harmonic viscosity, Adams-Bashforth III time-stepping with time step δt, and a linear free surface. Grid-scale noise is damped in the velocity and temperature fields over a δt timescale using a 4thorder (∇ 8 ) Shapiro filter, and the zonal filter is applied poleward of 45°t o damp small-scale longitudinal waves near the poles.

Domain
We used a latitude-longitude grid at two horizontal resolutions: low (L) with 256 × 128 grid points ( ∼ 1.4°) and = t 600s, and medium (M) with 512 × 256 ( ∼ 0.7°) and = t 300s. Runs began at resolution L and after equilibrating were increased to resolution M to resolve the deformation radius where N is the buoyancy frequency and H is the pressure scale height. According to Showman et al. (2011,  there. This means the lower boundary is impermeable (no mass flux through it) but physically moves up and down. The choice of bottom boundary is always tricky for Jupiter GCMs, as the flow at depth is not yet well known. There is nothing special about 18 bar, and the model is not directly forced towards observations, so a free-slip (barotropic) boundary condition is best in the absence of more information.
It is possible to run the MITgcm in "cubed-sphere" mode, and this has been done for giant planet and exoplanet work in the past Lian and Showman, 2010). This avoids grid point convergence at the poles, but has problems with energy and angular momentum conservation for Jupiter-like planets (Polichtchouk et al., 2014). Energy and AM conservation are particularly important for the long run time required by our radiation scheme (below), so we adopted the longitude-latitude grid in this case.

Jupiter parametrisations
We took the MITgcm atmospheric dynamical core and added a number of physical processes and parametrisations specific to Jupiter. Full technical details are given in Appendix A. Each of these adds a term into the right hand side of the prognostic equations, i.e. to F or Q in Eqs. 1 or 6. Physical parameters used in the model are listed in Table A.1.
We upgraded the radiation scheme (which contributes to Q) from the Newtonian cooling used by Zuchowski et al. (2009a) to a more realistic and flexible two-stream scheme. This was because Newtonian cooling could not decouple the effects of differential and interior heating, and also to provide more flexibility for future model development, such as radiatively active clouds. The model must run for much longer than before, however, as the two-stream approach implies realistic radiative time scales, which are O(decades) near the model base. The scheme simulates the absorption of shortwave (SW) radiation by the atmosphere, the absorption and emission of longwave (LW) radiation, and an optional interior LW heat source, under the assumption that the peaks of the SW and LW black body curves for the Sun and Jupiter are sufficiently far apart to be treated separately. There is no scattering. When included, the interior heating of 5.7 W m 2 is uniform everywhere. Differential heating between equator and pole is represented by a realistic annual-mean incoming solar radiation (ISR) profile ( Fig. 1). Optical depth profiles are prescribed as a function of p, linear in SW and quadratic in LW. The two-stream approximation is a fairly standard approach to simple radiation modelling; our scheme is largely based on Heng et al. (2011, for shortwave) and Mendonça et al. (2015, for longwave), and was developed for use in an idealised GCM; Tabataba-Vakili et al. (2019) describes the scheme in full detail. The Jupiter implementation has a few changes from that version, and in many respects is similar to the method used by Schneider and Liu (2009).
Vertical mixing of horizontal velocity components (which contributes to F) is included as a vertical diffusion term, where the Fig. 1. Annual mean incoming stellar radiation at Jupiter, the boundary condition for downwards shortwave flux at the top of the atmosphere. This profile already takes the planetary (bond) albedo into account. R.M.B. Young, et al. Icarus 326 (2019) 225-252 diffusion coefficient is based on the local Richardson number, Ri. This scheme parametrises sub-grid scale vertical mixing processes such as Kelvin-Helmholtz (Ri < 0.25) and symmetric (Ri < 1) instabilities. Adjustment of the temperature column due to dry convection is represented by a simple scheme derived from Manabe et al. (1965) (which contributes to Q). Where convective instability is found (dθ/dp > 0), an unstable column is relaxed towards a neutrally stable profile over a prescribed timescale, conserving the total enthalpy/ potential energy of the column. The Galileo entry probe found temperatures close to the dry adiabat below the cloud deck (Seiff et al., 1998), at least at the entry point, which was absent of water clouds (Ragent et al., 1998). There are no active moist processes in this version of the model.
The timescale for both vertical mixing and dry convection is the inertial timescale 1/f. Guided by Cessi (1996), who showed that instantaneous dry convective adjustment leads to unstable grid modes, we follow Young (1994) and adopt a realistic mixing timescale; he argued that for the mixed layer in the ocean = µ f ( ) 1 relates the timescale of vertical momentum mixing to the inertial frequency, and μ ∼ 1, so the vertical momentum mixing timescale is f 1 . On Jupiter, 1/f is typically around 2 h in midlatitudes.
At p > 0.8p base there is a linear drag at all latitudes towards rest (which contributes to F). Kinetic energy lost through drag is recycled into the temperature field (which contributes to Q). Ideally we would prefer not to need this drag, as the true lower boundary condition is poorly constrained by observations. The Galileo entry probe (albeit not at a representative point on the planet) did not find a reduction in flow speeds at this depth (Atkinson et al., 1997). However, we found that without this drag the zonal velocity at the equator accelerated to supersonic speeds and crashed the model.
While this drag was motivated by the representation of possible magnetohydrodynamic (MHD) drag from the planet's interior due to Ohmic dissipation of induced current by Liu et al. (2008) and Schneider and Liu (2009), we use it purely to ensure numerical stability. We note that Williams (2003b) did the same, albeit with a longer timescale. Using MHD drag as a physical basis for damping the flow well above the ∼ 100 kbar pressure where Ohmic dissipation is likely to be important has been controversial (Liu et al., 2008).
There is also a sponge layer in the top few layers of the model domain, where the eddy parts of the horizontal velocities and potential temperature are damped towards zero (which contributes to F and Q).

Description of runs and model spin-up
This paper focuses on two simulations, one with interior heating and the other without. The parameters for these runs are listed in Table 1.
The runs were initialised from rest and a θ profile generated by a radiative-convective model (RCM) containing the radiation and dry convection parts of the code only. This single-column model was run at each latitudinal grid point, starting with an isothermal column and finishing when the difference between the global mean flux through the top and bottom of the domain fell below ± 0.01 W m 2 . Typically this took 10 10 5 7 d, and the final profiles were independent of the initial temperature. All "days" are Earth days, exactly 86400 s. To break the horizontal symmetry, the θ field was then perturbed by uniformly distributed random numbers at all grid points with amplitude 0.01 K.
The simulations were first run at resolution L to a state of radiativeconvective-dynamical equilibrium. Once the runs had equilibrated at resolution L, to resolve small-scale features we doubled the resolution to M, halved the time step, and continued the run until the flow had equilibrated again. The model was considered equilibrated once (1) the net radiative flux into and out of the top of the model domain and (2) the total atmospheric KE were approximately constant in time. Fig. 2 shows the top-of-atmosphere net radiative flux and the total atmospheric KE over time for runs A and B. The model runs required are long because the radiative time constant varies approximately linearly with pressure, and the model base is at 18 bar. Schneider and Liu (2009), who put their model base at 3 bar, needed about 30000 d for their runs to equilibrate, which is broadly consistent with this trend.
The simulations used 64 and 128 cores for the resolution L and M runs respectively. The MITgcm is documented to scale to thousands of cores (Hill et al., 2007), but unfortunately our current configuration scales much more poorly than this. Because we use the longitude-latitude grid, we had to include a zonal filter to satisfy the Courant-Friedrichs-Lewy condition near the poles, where the longitudinal grid size approaches zero; without filtering the model crashes within a handful of timesteps. The zonal filter requires variables along each latitude circle to be gathered to a single processor and this scales extremely poorly with the number of longitudinal cores. We intend to return to this as our model could be run much more efficiently if this is resolved. Unlike Dowling et al. (1998), who faced similar challenges with the EPIC giant planet model, we found that it was more efficient to use two or four cores along each latitude circle, rather than one.   . Note the colour scales are not necessarily the same in the two runs, as the ranges are quite different in most quantities. The vertical velocity is = w g /( ). R.M.B. Young, et al. Icarus 326 (2019) 225-252 The optimal use of resources (about 50% of linear speedup) at the two resolutions is shown in Table 2. The total cost of a run was typically about 75000 core-hours. The general appearance of both Runs A and B is a banded structure with eastward jets at some latitudes and westward jets at others, extending throughout the vertical domain. The meridional circulation extends throughout the domain in both cases, with significantly stronger circulation in the equatorial troposphere than elsewhere (note the logarithmic colour scale in Fig. 4b).

General appearance of the flow
Run A is characterised by a strong retrograde jet at the equator, which extends almost to the bottom of the domain, while Run B has a prograde equatorial jet in the troposphere and a retrograde stratospheric jet ( Fig. 3a-b, 4a). The jets in Run A are significantly stronger than in Run B. In Run B the magnitude of the meridional circulation is weak above around 300 mbar, but is approximately constant with height in Run A, except near the equator above 1 bar, where it peaks. The shortwave optical depth equals one at 1 bar, so the peak in solar heating will be around this pressure, triggering convective instability above and driving a mean meridional circulation. The flow has a strong barotropic component (Fig. 4a) and only prograde flow in the stratosphere away from the equator. Most of the overturning circulation (Fig. 4b) is confined to the neutrally stable part of the atmosphere below ∼ 300 mbar, which can be seen in Fig. 5 and 6. In this region there is overturning circulation whose vertical extent is bounded by the tropopause at the top. The strongest meridional circulation is confined to the equatorial region in Run B, which is about 100 times stronger than at other latitudes. The off-equatorial meridional circulation in Run B is comparable with the equivalent circulation in Run A. At low latitudes in Run A the meridional circulation is generally weak, confined to a region near the tropopause with weak  R.M.B. Young, et al. Icarus 326 (2019)   circulation above and below. Note that, while there is velocity damping at the bottom level of the model, the zonal velocities do not actually go to zero at the bottom of the domain. There are vortices away from the equator in Run A, and over the whole globe in Run B . Their size distribution is fairly narrow, and more restricted than one might expect for 2D turbulence and from observations. There are no large-scale isolated vortices like the Great Red Spot (GRS, which covers a 30°× 20°area) in either run. We have not seen such vortices generated in any of our runs to date, and are not aware of any GCM of Jupiter's atmosphere that can generate vortices the size of the GRS spontaneously outside the polar regions. While the GRS has been present for possibly up to 350 yr (Hooke, 1665), models generally over-damp the atmosphere such that large vortices like the GRS dissipate more quickly than in the real atmosphere. A vortex the size of the GRS requires a potential vorticity (PV) anomaly that bends the lines of constant PV along latitude circles considerably (Marcus and Shetty, 2011). In models with significant diffusion and/or strong meridional shear, the lines of constant PV may be too stiff to maintain such an anomaly against diffusion for the time required for a large vortex to develop spontaneously, at least at low latitudes. Modelled vortices typically decay over several turnaround periods or radiative timescales, which is much less than the observed lifetimes of some of these features, and this phenomenon is not yet well understood.
Clearly some of these vortices are part of longitudinal wave trains, particularly at midlatitudes in Run A and throughout Run B. The vertical velocity field in Run A (Fig. 3e) shows several particularly clear wave trains in the prograde jets at midlatitudes. These waves have the so-called "chevron" pattern characteristic of waves in Jupiter's atmosphere (Simon-Miller et al., 2012), which are indicative of Rossby wave activity at those latitudes (Vallis, 2006, Ch. 12). Such patterns can also be seen throughout the relative vorticity field in Run B (Fig. 3d). The 512 × 256 resolution used in these runs is not high enough to resolve the ∼ 300 km mesoscale waves seen by Flasar and Gierasch (1986), Reuter et al. (2007), and Simon et al. (2015).
The "chevron" pattern is characterised by vortices aligned with the latitudinal shear such that eastward eddy AM is carried into latitudes with eastward jets, and westward eddy AM is carried away (vice-versa for westward jets). Fig. 7 zooms into the midlatitudes of both runs. Eddy momentum flux (EMF) 〈u′v′〉 and zonal mean zonal velocity 〈u〉 show eastward momentum being advected into the eastward jet, aligned with the "chevron" pattern. 〈·〉 denotes a zonal mean, i.e., along a latitude circle, and ′ the deviation therefrom. The flow is cyclonicallybiased north of the eastward jet, and anticyclonically-biased south of the jet. On average, therefore, eastward eddy momentum converges at both eastward jet flanks, accelerating both flanks of the jet.

Thermal structure and radiative-convective model
The temperature field is also banded with a strongly zonal structure in both cases . However, the vertical structures are quite different. Fig. 4c shows zonal mean temperatures for the two runs, and The upper atmosphere (above ∼ 200 mbar) is similar in both runs. In this region, radiation dominates over convection. The global minimum in temperature occurs over the poles around 0.1 bar, where the temperature falls as low as 50-60 K. The tropopause is close to 0.1 bar in Run A at all latitudes and in Run B at low latitudes; it is slightly higher at high latitudes in Run B, and points close to the pole have a tropopause around 0.03 bar. This tropopause height is consistent   showing the link between the "chevron" patterns and the zonal jets. Each panel shows the zonal mean zonal velocity (solid), eddy momentum flux (dashed), and vertical velocity. Dark is upwards, with a thick white contour at zero. R.M.B. Young, et al. Icarus 326 (2019) 225-252 with most other thick atmospheres in the Solar System (Robinson and Catling, 2014).
In the lower atmosphere, however, the interior heating in Run B has a dramatic effect on the temperature structure. Without interior heating the initial radiative-convective equilibrium is statically stable at the bottom of the atmosphere, leading to a temperature inversion (∂T/∂p < 0) in Run A in the deep atmosphere (Fig. 5, left). The interior heating in Run B pushes the lower atmosphere towards convective instability, however, and so the deep atmosphere becomes neutrally stable. As a result the latitudinal dependence of the temperature structure is quite different. Fig. 6 shows profiles of buoyancy frequency at a selection of latitudes in both Runs A and B. When the atmosphere is heated from below, dry convection is triggered and there is a latitudinally uniform tendency of the deep atmosphere in Run B towards neutral stability. This results in a smaller equator-to-pole temperature difference (ΔT ep < 10 K) in the deep atmosphere in Run B, as the atmosphere largely ignores solar heating in maintaining the deep temperature structure -a consequence of the so-called Jovian "thermostat" (Ingersoll, 1976;Ingersoll and Porco, 1978). By contrast, there are large tropospheric equator-to-pole temperature differences in Run A.
Our temperature profile in run B is closer to neutral stability than the Galileo probe data shows (Seiff et al., 1998;Allison and Atkinson, 2001;Magalhães et al., 2002). The probe found an atmosphere about 0.1 K km 1 more stable than the dry adiabat between 3-8 bar. A slightly unstable region between 8-14 bar would probably become stabilized if contributions from molecular weight gradients due to water were fully accounted for, and by the probe's final depth of 22 bar the atmosphere was significantly stable. Evidence from wavefront expansions after the Shoemaker-Levy 9 impact also suggests that the deep troposphere is statically stable (Ingersoll and Kanamori, 1995). A lack of deep static stability in run B may explain some of the differences between the model and the real planet. We shall explore forcing the deep troposphere towards a slightly stable profile in future.
While our model's deep temperature structure has a weak equatorto-pole temperature difference, Jupiter's observed outgoing longwave radiation and associated brightness temperatures are an even flatter function of latitude (Ingersoll, 1976). In this model we have imposed uniform interior heating everywhere. We expect that, if we used the latitudinally-varying interior heat flux deduced from Voyager observations by Pirraglia (1984), where the extrapolated interior heat flux increases from 4 to 12 W m 2 from equator to pole, we might obtain an even flatter tropospheric temperature structure in Run B. We notice there is a zonal velocity maximum in the upper atmosphere near the poles in Run B. This is consistent with thermal wind balance, as ∂T/∂y is small near the poles in the lower stratosphere but increases above (where = y a is the absolute distance in the meridional direction). If the interior heat flux varied with latitude with a consequent temperature profile closer to the observed profile, however, then these jets might disappear. Fig. 8 presents a comparison between the initial radiative-convective model state and the final temperature structure. In both cases the dynamics have a noticeable effect. The contribution to the radiative-convective-dynamical equilibrium at different levels is highlighted in the latitudinal profiles in Fig. 8b,d. There is an important dynamical component at high latitudes in all cases, with the amount of heating due to the dynamics increasing deeper into the atmosphere. For example, in Run A at 0.1 bar the temperature profile follows the radiative equilibrium profile everywhere except near the poles where there is 30 K of heating due to the dynamics. This increases to about 80 K in the deep atmosphere. In Run B, however, the dynamics affect the temperature profile at all latitudes, with an increasing cooling effect at low latitudes and an increasing heating effect at high latitudes the deeper one goes.
Both radiative-convective equilibrium profiles have large negative values of ∂T/∂p between the top two levels of the model (Fig. 8a,c). This inversion persists in the GCM runs in both cases. This appears to be a consequence of heating in the top model layer by the radiation scheme when the top of the uppermost grid cell is at vacuum. We experimented with the positions of the uppermost few model levels in the RCM, but with the top of the domain at vacuum we were unable to get rid of this inversion. Unfortunately the MITgcm in its current form does not support an atmospheric model with its top at p > 0, as far as we are aware, but this could be looked at in future. As the top few levels are within a sponge layer the uppermost few model levels should be treated with caution anyway.

Characteristics of the zonal jets
At the cloud level (0.3-3.0 bar) there are about 34 jets in Run B (identified by eye), of which 18 are eastward, and about 24 in Run A, of which 14 are eastward. Fig. 9 shows zonal mean zonal velocity profiles for the two runs at the same pressure as Fig. 3.
In Run A, the latitudinal jet scale varies strongly with latitude, with wider jets at low latitudes. The jet width in Run B varies more slowly with latitude. The jets in A are much straighter, in the sense that the percentage of KE in the zonal mean zonal flow as a fraction of the total KE is much higher in Run A (98.5%, at the level in Fig. 3) than in Run B (76.0%). The measured figure at Jupiter's cloud tops is 92.9% (Galperin et al., 2014), somewhere between the two.
The jet scale is often compared with the Rhines scale, the scale at which the RMS velocity and the Rossby wave phase speed are equal (Rhines, 1975). It can be defined in various ways; two are rms 2 2 suitably averaged, is a typical wind speed (Galperin et al., 2014), and 2 is the eddy kinetic energy (Chemke and Kaspi, 2015b). On Jupiter the Rhines scale is some 20000 km using the first definition, and 12000 km using the second (using figures for total and eddy kinetic energy in Galperin et al., 2014), and these are comparable with the observed jet scales.
We estimate the jet scale here as the meridional distance between adjacent peaks in zonal mean zonal velocity (Fig. 9). We estimated this at each latitude for both the eastward and westward peaks separately, and took the average. We also computed the two Rhines scales as a function of latitude, vertically averaged (in ln p) up to 0.1 bar, and these three scales are shown in Fig. 10. The two definitions of the Rhines scale generally agree except in the equatorial jet in Run A, where eddy activity is very weak. In both cases the Rhines scale increases with latitude, although L Rh,1 is approximately flat with latitude (ignoring jetby-jet variation) except near the poles. In Run A the jet scale and Rhines scales generally agree at midlatitudes, but near the equator the jet scale is significantly larger than the Rhines scale. In Run B the jet scale peaks at midlatitudes, falling towards both the equator and the poles, and while generally larger than both L Rh,1 and L Rh,2 , they agree to within a factor two.
This model is not a complete representation of Jupiter's atmosphere, and is not forced by observations explicitly, so one does not expect quantitative correspondence between the model and observations, but nevertheless qualitative comparisons can be made based on our more realistic case (Run B). In general, the number of jets is a reasonable match to observations. However, jet speeds are weaker than observed, particularly at low latitudes, where observations show jets with speeds 100 m s 1 (Porco et al., 2003, and Fig. 9). Furthermore, the real planet has strong hemispherically asymmetric jets around 20°N/S. As the forcing of our model is symmetric about the equator, we are not surprised that it does not reproduce this strong asymmetry. R.M.B. Young, et al. Icarus 326 (2019) 225-252 Flasar et al. (2004) found a 140 m s 1 prograde equatorial stratospheric jet, but this is above the top of the domain in our model. Offequator the observed stratospheric jets are primarily barotropic, with weakening of the jets above 100 mbar equatorward of 30°N/S. In our model the jets are also primarily barotropic, but generally the zonal velocity increases into the stratosphere. In the real atmosphere stratospheric temperatures increase towards the pole (Flasar et al., 2004, Fig. 1a), so zonal velocity goes down with height. But this is likely due to physical processes that are not included in our model (e.g. gravity wave breaking or joule heating). Here stratospheric temperatures decrease towards the pole as only solar heating is included as a physical parametrisation in the stratosphere, and so we expect zonal velocity to increase with height according to the thermal wind equation (and not necessarily so near the equator, where that balance breaks down).
While Flasar et al. (2004) do not see a retrograde equatorial jet in the lower stratosphere, Read et al. (2006) do using the same dataset (although they used a simpler version of the thermal wind equation and applied this balance closer to the equator). They see a retrograde equatorial jet of unknown magnitude between 30-100 mbar, which has similar latitudinal width to the equatorial stratospheric jet in our   (Porco et al., 2003, dotted, red). Model profiles are averaged over the same times as Fig. 4 (100 d, at or near the end of the run). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)  R.M.B. Young, et al. Icarus 326 (2019) 225-252 model. The jet right at the top of the model domain is likely to be affected by both the top of the model domain and the sponge layer, which suppresses eddy activity in the uppermost few model levels.
Below the clouds, the Galileo probe showed increasing wind speed down to 4 bar, and approximately constant (barotropic) wind between there and 20 bar. Our model shows an approximately barotropic structure in the lower part of the domain, but not the rapid increase in the few bars below the cloud level.
We calculated the net power from the eddies into the zonal jets as a function of latitude and height. The power from the eddies into the jets is 〈u′v′〉∂〈u〉/∂y, where 〈u′v′〉 is the EMF and ∂〈u〉/∂y is the latitudinal gradient of zonal mean zonal velocity. When these two quantities are correlated the tendency is for energy to be transferred from the eddies to the zonal jets, and vice versa when they are anti-correlated (Salyk et al., 2006).
In Run A, the eddy to zonal power is weak at low latitudes, and almost all of the KE is in the zonal flow (Figs 3a, 4a). Averaging globally we obtained a positive net flux of eddy-to-zonal flow power over the final 5000 d of the two runs, with variability over ∼ 300 d in both cases. In Run A this power is 0.06 ± 0.01 W m −2 , mostly confined to the mid-and high-latitudes. In the interior heated case Run B the net power from the eddies into the zonal jets, averaged over the whole domain, is 0.115 ± 0.009 W m −2 , or about 0.9% of the emitted infrared heat flux. This power occurs primarily in the upper troposphere at latitudes within 60°of the equator, which also occurs on Earth (Ait-Chaalal and . However, the total power per unit area is somewhat less than the 0.7-1.2 W m −2 inferred by Salyk et al. (2006). The difference between the modelled and observed winds, which differ by a factor of five at the equator, is the likely cause. Maximum transfer from eddies into the jets occurs in the regions of strongest zonal shear (both cyclonic and anticyclonic), as in observations (Salyk et al., 2006, Figs 4-5), and in those regions the typical power is ± × (4 2) 10 W kg , 5 1 which is within a factor of two of the observed values of × (7 12) 10 W kg 5 1 (Salyk et al., 2006). The stability of a flow on a rotating sphere with respect to symmetric instability is diagnosed by the relative steepness of the potential temperature (θ) and absolute angular momentum (M) contours, where = + M r u r ( cos ) cos and r is the radius of the planet. Interest in this instability in the context of giant planets has been recently revived through studies of 'slant-wise' convection by O'Neill and Kaspi (2016) and O'Neill et al. (2017). If the potential vorticity and the Coriolis parameter have opposite sign, then the flow is unstable to symmetric instability. This condition is equivalent to θ contours being steeper than M contours.
We computed M and θ at the end of Runs A and B, and these are shown in Fig. 11. Run B is strongly stable with respect to symmetric instability, with the angle between the M and θ contours close to the latitude angle throughout the domain; only right at the equator do the M contours approach the θ contours. The peak in absolute AM around 0.5-1.0 bar at the equator is due to the peak in the zonal velocity there (Fig. 4a). Run A is similar over most of the lower troposphere and stratosphere, but between 0.3-1.0 bar the θ contours steepen considerably, such that on the flanks of the subtropical eastward jet the flow is close to becoming symmetrically unstable. The same occurs in the southern hemisphere. Such an instability would cause motion along is entropes in this region to become unstable, mixing momentum in a close-to-vertical direction, and therefore might contribute to the weak vertical gradients in the jets (Fig. 4a). We note that the timescale for symmetric instability and convective adjustment are similar (Emanuel, 1985), so the dry convective adjustment in our model may circumvent any symmetric instability that might occur in Run B.

Meridional structure of the midlatitude jets
The balance of forces driving a zonal jet can be diagnosed using the Eulerian mean equation for the zonal flow (Andrews and McIntyre, 1976;Andrews et al., 1983) (9) to which we have added terms present in the model due to vertical diffusion, bottom drag, and an estimate of the zonal velocity tendency from the Shapiro filter, with = t shap (an approximation as the filter is an adjustment). The final term, in 〈v′∂v′/∂λ〉, is zero in the continuous equations but is small but nonzero due to the vector invariant discretisation used in the MITgcm. It is similar in appearance to the horizontal metric term uv a / , which is explicitly excluded from the continuous form of the vector invariant equations. See Supplementary Data S4 for the full form of each of the momentum terms using this discretisation. We do not interpret our results using the transformed Eulerian mean (Andrews and McIntyre, 1978) or the Eliassen-Palm (EP) flux (Edmon et al., 1980), as ∂〈θ〉/∂p is close to zero in the troposphere, and so the transformed Eulerian mean and EP fluxes are difficult to interpret. Fig. 12 shows the midlatitude zonal jet structure in equilibrium near the end of Run B. The same basic structure is present in all the offequatorial jets in this run. The midlatitude zonal jets are forced by a convergence of eddy angular momentum fluxes with a similar structure to a Ferrel cell, which is characteristic of baroclinic instability, at least Fig. 11. Potential temperature (black, solid) and absolute AM (grey, dashed) contours averaged over the same time periods as in Fig. 4. Note the different latitude ranges in the two runs. R.M.B. Young, et al. Icarus 326 (2019) 225-252 in Earth's midlatitudes. Momentum balance in an eddy-driven Ferrel cell is shown schematically in Vallis (2006, Fig. 11.15), and we may make a direct comparison with the jets here. We examine the cell between 34-38°N associated with a prograde jet. First, the sense of the mean flow is indirect, rising at higher latitudes and falling at lower latitudes (Fig. 12a). Second, on the vertical branches of the cell the vertical velocity and the meridional eddy heat flux convergence ( v T a ( cos )/( cos ), Fig. 12b) have the same sign.
Third, the upper branch of the cell must balance Coriolis accelerations (Fig. 12c) against meridional eddy AM flux convergence (Fig. 12d). At first glance these terms do not appear to balance. However, vertical diffusion acts to cancel the Coriolis acceleration almost completely (see Supplementary Figure S1 for this term and the others in Eq. 9). Hence the dynamically relevant quantity is instead an adjusted Coriolis acceleration that is the sum of these two terms (Fig. 12e). This balances the meridional eddy AM flux convergence on the upper branch of the cell, as required.
Finally, the lower branch of a Ferrel cell balances Coriolis accelerations against surface drag, or in this case our bottom drag (Fig. 12f). At the model base the drag is negative and balances the adjusted Coriolis acceleration (Fig. 12e). These four properties of the cell show that the midlatitude prograde zonal jets are structured like Ferrel cells, once vertical diffusion is taken into account, and that the jets are driven by meridional eddy AM flux convergence.
In the retrograde jets the balance is similar, but with a wider cell and with the sense of each term reversed, i.e. a direct meridional circulation where eddies drive the zonal flow. For the retrograde jet at 30°N, the rising branch is the northern branch of the Ferrel cell to the south (24-26°N), and the descending branch is the southern branch of the Ferrel cell to the north (described above). This can be seen more clearly if we extend our adjusted Coriolis acceleration to an adjusted zonal mean meridional velocity and hence an adjusted meridional mass stream function. This is shown in Fig. 12g, which shows the sense and the width of the overturning circulation associated with both the prograde and retrograde jets. Within the regions of greatest meridional shear there is a weak vertical flow (e.g. at 28°and 32°N in Fig. 12a). These lead to a convergence and divergence of vertical eddy momentum fluxes on the flanks of the prograde jet (Fig. 12h). The prograde jet is accelerated on its equatorward flank and decelerated on its poleward flank. Furthermore, the meridional eddy momentum flux convergence is also biased slightly equatorward of the centre of both the eastward (less clear) and westward (more clear) jets. Both of these asymmetries are reflected in the net acceleration of the zonal flow (Fig. 12i).
In Run A the balance in the prograde jets at high latitudes is similar, although slightly less clear as the jets are not as clean as in Run B.  Figure S1. The heat flux uses T rather than θ to improve the visibility of the cells in the lower atmosphere. R.M.B. Young, et al. Icarus 326 (2019) 225-252 Supplementary Figure S2 shows the same terms as Fig. 12 for this run. Within the prograde jets there is a thermally indirect overturning circulation with the correct balance for a Ferrel cell. Note, however, that in this case vertical diffusion is confined to the equatorial region and hence an appeal to the adjusted Coriolis acceleration is not required. While there is an asymmetry in the vertical eddy momentum fluxes as in Run B, the net acceleration of the jets does not show any obvious bias across the jets.

Time evolution of the flow, and final state at resolution L
The full time evolution of the flow reveals characteristics that will be important later on in the context of equatorial super-rotation. Fig. 13 shows Hovmöller diagrams of the zonal mean zonal velocity over the entire length of Runs A and B at 1 bar. An equatorial jet is established almost immediately during Run A, and it maintains its speed throughout. In Run B the prograde equatorial jet also appears early on (after about 3000 d), but its magnitude then decreases and is weak throughout the resolution L period. In both cases the zonal jets are not fully resolved at resolution L, and become better defined once the resolution doubles. Fig. 14 shows the zonal velocity and temperature at the end of the low resolution part of both runs. Both show less clearly-defined zonal jets at resolution L than at resolution M (Fig. 4a), with easterly jets in Run B only at the bottom of the domain and in the equatorial stratosphere. The temperature in Run A is similar to resolution M, but we notice there is a local minimum in temperature at the equator in Run B. Other differences (not shown) in the low-resolution part of Run B are that the eddy momentum flux and the power from the eddies into the jets are largely confined to the stratosphere. This is not the case in Run A, where it is the same throughout, although both quantities are significantly weaker at resolution L. Supplementary Data S3 shows the evolution of the zonal mean zonal velocity (Fig. 4a) over the whole of Run B.
It is not uncommon in simulations of planetary atmospheres for zonal jets to migrate in latitude (Williams, 2003b;Chan et al., 2007;Chemke and Kaspi, 2015a). Fig. 13 clearly shows that in Run B the zonal jets are migrating towards the equator, at approximately  Fig. 14. As Fig. 4, but averaging over the 100 d at the end of the low resolution part of the run, immediately before the transition from resolution L to M. R.M.B. Young, et al. Icarus 326 (2019) 225-252 1.3 cm s 1 . In Run A any migration is slow and with no obvious pattern. Even at resolution L the under-resolved jets migrate towards the equator throughout, at a rate that does not appear to change once the resolution is doubled. During the resolution L period there is some short-period migration of jets towards the equator on top of the longer period migration, but this disappears once the resolution doubles. The end of Run B occurs during a merger of a southern hemisphere jet into the equatorial jet. As a result, that time is not representative of the equilibrium state of this run, so throughout the paper we have calculated diagnostics for this run at 133000 d, shortly before the merger occurs, rather than at the end, except where stated explicitly.
Williams (2003b, in Jupiter's atmosphere) and Chan et al. (2007, in Earth's oceans) both see equatorward migration of zonal jets, while Chemke and Kaspi (2015a) see poleward migration. Chan et al. (2007) attributed the migration to a residual circulation generated by the EP flux, which was convergent on the poleward flanks of equatoriallypropagating jets, and divergent on the equatorward flank. As there was no obvious asymmetry between the eddy momentum fluxes on the two flanks of the jets, they related this to the vertical component of the EP flux, finding an asymmetry in the baroclinic eddy activity consistent with a poleward bias in the Eady growth rate = fU N / , z where U z is the vertical shear in the zonal flow. They traced the migration to a trend in the static stability, which increased towards the equator, generating an asymmetry in baroclinic eddy activity across the jets.
In our Run B there is a clear asymmetry in the eddy AM flux convergence across the midlatitude jets. Throughout both runs the peaks in the meridional eddy momentum flux convergence coincide with the baroclinic zonal jets, showing they are continuously eddy-driven; this is shown in Supplementary Figure S3. The asymmetries in both the meridional (Fig. 12d, mostly in the westward jets) and vertical (Fig. 12h) eddy AM flux convergence force the jets towards the equator. These asymmetries are either absent or intermittent in Run A (see Supplementary Figure S2). Comparing the positions of the midlatitude jets and the peaks in the Eady growth rate (Fig. 15), Run A has no obvious asymmetry in the Eady growth rate across the jets, while Run B has a clear poleward bias across each jet. Typical Eady growth rates at this level are 0.5-2.0 d 1 for Run A and 0.1-1.0 d 1 for Run B, generally larger than but not unlike the Earth's atmosphere (Vallis, 2006, Ch. 6).
The only difference between the two runs is the interior heat flux, which has a critical effect on the static stability (Fig. 6). Fig. 16 shows the variation in static stability with latitude for both runs at the two pressure levels where there is the largest asymmetry across the jets in Run B. In Run A there is no obvious trend from equator to pole, but in Run B there is a clear trend for the static stability to increase towards the equator. The variation in Eady growth rate across the jets and the upwards trend in static stability towards the equator suggests that Run B behaves consistently with the mechanism described in Chan et al. (2007). Chemke and Kaspi (2015a) also find a poleward bias in the Eady growth rate across the jets, but a poleward migration of the jets, so the mechanism acting must be different. They find the static stability is symmetrical about the jets, and instead find the main trend contributing to the bias in baroclinic growth across the jets is the Coriolis parameter.
To understand why the trend in static stability exists, we plot potential temperature profiles at three latitudes in Fig. 17. In Run A the profile above ∼ 3 bar is similar in form at each latitude, because there is no constraint from below, just shifted depending on the solar flux. Hence the vertical potential temperature gradients will be only weakly latitudinally dependent, if at all. In Run B, however, the profile below ∼ 0.3 bar is similar at all latitudes, as the interior heat flux constrains the tropospheric temperature profile towards a dry adiabat, and so the deep profile is only weakly dependent on the solar flux. This weak dependence can also be seen in equivalent profiles for buoyancy frequency in Fig. 6. Above, where solar flux dominates, the latitudinal variation is similar to Run A. To transition between the two regions, therefore, the vertical potential temperature gradient above the region    R.M.B. Young, et al. Icarus 326 (2019)   of neutral stability must increase towards the equator, and hence the static stability must increase towards the equator, as we see in Fig. 6 and 16. We note that making the interior heat flux a function of latitude, as mentioned in a previous section, should strengthen this trend.

Spin-up and maintenance of the equatorial jets
In this section we first show that longitudinal waves with the properties of Rossby waves exist at the equator in Run B, but not in Run A. Then we relate these waves to equatorial instabilities in the flow during Run B, which leads to the spin-up of a super-rotating equatorial jet, and compare this to how the sub-rotating equatorial jet in Run A is generated. We also diagnose how the equatorial jets in both runs are maintained in the model's equilibrium state.

Equatorial waves
In Run A the most prominent waves are at mid-latitudes (Fig. 3e), while in Run B there are waves throughout the domain. The equatorial region is of particular interest due to the possible role of Rossby waves in driving equatorial super-rotation.
The theory of linear equatorial waves on a β-plane was developed by Matsuno (1966), who showed that longitudinal wave solutions to the 2D barotropic equation on an equatorial β-plane have dispersion relations corresponding to the solutions of for integer = … n 0, 1, 2 , where = c gh is the pure gravity wave speed and h is the equivalent depth of the wave. This has 1-3 physical solutions ω(k) for each n 1. The meridional structure of the wave is proportional to the Hermite polynomial H n for each n, such that waves with odd n are symmetric about the equator, and waves with even n are antisymmetric. Waves decay away from the equator with a length scale given by the equatorial deformation scale = L c/(2 ) eq (Gill, 1982). Solutions correspond to equatorial Rossby waves (n ≥ 1, low ω), inertia-gravity waves (n ≥ 1, high ω), mixed Rossby-gravity waves ( = n 0), and the special case = n 1 with zero meridional velocity everywhere, corresponding to Kelvin waves = ck. Gill (1982, pp. 434-440) covers the various dispersion relation solutions in some detail, and the horizontal structure of these solutions can be seen in, e.g., Matsuno (1966) and Kiladis et al. (2009). Fig. 18 shows horizontal velocity and temperature within 10°of the equator for Runs A and B in equilibrium. The zonal mean has been removed in order to see the equatorial wave structure. Eddy velocities in Run A are much smaller than in Run B (compare the eddy velocity in A with the zonal mean zonal velocity, around 50 m s , 1 while the eddy velocity in B is comparable with the zonal mean zonal velocity, around 20 m s 1 ). Therefore any equatorial waves in Run A will be weak, while waves in Run B will be of comparable magnitude to the zonal mean flow. Second, the latitudinal structure of the eddy fields appears to be broadly symmetrical about the equator in Run A, while in Run B it is less clear, but all fields appear to be more symmetric than antisymmetric.
The propagation of these waves over time is shown as Hovmöller diagrams in Fig. 19. Supplementary Figure S4 shows the same at ± 1.05°latitude. Zonal velocity is shown as it has the clearest pattern, but other quantities are similar. The visible waves in Run A are eastward relative to the mean flow. The dominant mode moves with a phase speed around 40 m s 1 (estimated by eye), and there is a weaker wave moving with a phase speed around 7 m s 1 . In Run B the waves are westward relative to the mean flow. The dominant modes are a stronger wave with phase speed around 8 m s , 1 and a weaker wave with phase speed around 40 m s 1 .
To identify the type of wave(s) present, we computed the frequencywavenumber power spectrum to identify which equatorial wave dispersion relations would fit our data. We followed the procedure described by Wheeler and Kiladis (1999), adapting their code to do so 1 . The analysis was done using the zonal velocity, as above. We began with zonal velocities at = p 1 bar for the same 512 d as Fig. 19, at the four latitude grid points closest to the equator ( ± 1.05°and ± 0.35°). These four latitudes had similar zonal mean zonal velocities and so the effects of shear in the zonal mean (which is not accounted for in the standard theory of equatorial waves) is minimised, although not removed entirely. First we transformed to the frame of reference moving with the mean flow, as in Fig. 19. At each longitude we removed the time mean and the linear trend in time, and tapered the ends of the time series to zero over five days, to satisfy the FFT condition for periodicity in time. We then computed the complex FFT in longitude and time to give the 2D spectrum as a function of zonal wavenumber k (m in non-dimensionalised units) and frequency ω.
The power spectra for the gravest ∼ 70 zonal wavenumbers and wave periods longer than ∼ 2.4 d are shown in Fig. 20. Wheeler and Kiladis (1999) split the 2D power spectrum into its symmetric and antisymmetric parts before normalising the spectrum by a background spectrum to reveal the wave modes present above the noise. However, this was not necessary here, as the differences between the two runs and the types of equatorial waves that are present are clear just from the raw 2D power spectrum.
We fitted (by eye) various dispersion relations to the dominant peaks in the 2D power spectrum (Gill, 1982, pp. 438-9). In Run A, most of the power is in the Kelvin modes. The Kelvin wave dispersion relation = ck is best fit to the data with equivalent depth = h 53 m, corresponding to a phase speed of = v 37 m s p 1 . This dispersion relation passes through peaks in the spectrum at = m 13 and 41, and corresponds to the fast mode in the Hovmöller diagram (Fig. 19a) (Gill, 1982, Eq. 11.6.8). Most of the power is in a broad band that can be fit to the symmetric component = n ( 1) with equivalent depths between = h 15 and 80 m. The peaks in the spectrum are at = m 2 3 and = m 8, with frequencies that can be fit using equivalent depth = h 40 m. Table 3 shows the properties of the waves derived from these fits. The slow mode identified in the Hovmöller diagram corresponds to the waves with phase speeds between −6 and 15 m s , 1 with the wave at = v 10 m s p 1 dominant. The fast mode in the Hovmöller diagram is likely to be the mode with equivalent depth = h 550 m. There is no dominant mode in this case, but = m 8 contains more power than the other modes (while spread over a range of frequencies). The phase speed at that wavenumber is = v 36 m s p 1 , close to the phase speed estimated by eye from the Hovmöller diagram. The equatorial deformation scales derived from these modes are not unreasonable.
Our equivalent depths are considerably less than those deduced from observations, which are typically O(km) (Allison, 1990;Simon-Miller et al., 2012). Our phase speeds are typically a factor five smaller than observed, and the equivalent depth varies as v p 2 (for Kelvin waves, and Rossby waves at low zonal wavenumber). It is not clear what could explain a true difference of this magnitude in the Rossby wave speed, so this is something we could focus on in future work. However, we note that Allison (1990) measured the equivalent depth and hence the phase speed from the vertical wave structure. Instead, we computed the phase speed directly from horizontal wave propagation, basing the equivalent depth on that, with no consideration of vertical structure. Our phase 1 https://github.com/UV-CDAT/wk. R.M.B. Young, et al. Icarus 326 (2019) 225-252 speeds might be more consistent with observations if the vertical structure were accounted for in the same way. For completeness, dispersion relations for other solutions to the wave equation (inertia gravity waves and mixed Rossby-gravity modes in both cases, and Kelvin waves in Run B) are also superimposed on Fig. 20, with equivalent depths used in the other fits, but there was no evidence for these waves. There was some power in the equatorial Rossby wave part of the spectrum in Run A, but the power in the Kelvin modes was much higher.
Finally, in Run A there is a peak in the power spectrum that is fit well by a dispersion relation with phase speed 51 m s 1 . This is almost certainly a remnant of the mean flow, as the mean flow is at this speed, and it has power at all wavenumbers. This suggests that the signal of the zonal mean flow is not perfectly removed, for example due to shear in the mean flow, changes in the mean flow over time, or nonlinear effects. We plotted the same for Run B ( 17 m s ) 1 but it is not clear there is anything there in that case. Schneider and Liu (2009) attribute the source of Rossby waves at the equator to convection from the deep atmosphere. This drives temperature fluctuations that, if not balanced by slow radiative processes, perturb the horizontal divergence, which is a source of vorticity fluctuations and hence Rossby waves. These radiate away and dissipate at higher latitudes, transporting eastward AM into the equatorial region and driving super-rotation.

Spin-up of the super-rotating equatorial jet in Run B
In our simulations, the previous section demonstrated that Rossby waves exist at the equator in Run B, which can provide a source of eastward AM flux there. The source of Rossby waves, however, is unlikely to be the mechanism described by Schneider and Liu (2009). The convective timescale in our model scales as f , 1 which means that dry convection, while present at the equator, is weakest there. It seems unlikely to be sufficient to stimulate enough Rossby wave emission at the equator to drive super-rotation.
Instead, the mechanism is related to the migrating zonal jets described in Section 4.5. At the beginning of the resolution-M part of the simulation, the off-equatorial baroclinic jet structure described in Section 4.4 forms within 100 d (recall the typical Eady growth rate is 0.1 1.0 d 1 ), with a westward jet at the equator, in the absence of any source of angular momentum flux convergence there (Fig. 13b). Due to the dependency of the baroclinic growth rate on latitude across each off-equatorial jet (Fig. 15), these jets immediately begin to migrate towards the equator. Fig. 21 shows a series of equatorial snapshots over the following 1000 d (Supplementary Figure S5 shows the full sequence over 1500 d). As the eastward jets approach the equator (104100-104400 d), maintaining their strength, the westward jet at the equator also maintains its strength, but the gap between the eastward jets on either side of the equator decreases. The flanks of the eastward jets closest to the equator are squeezed into a smaller and smaller space, increasing the latitudinal shear and hence ∂ 2 〈u〉/∂y 2 at the equator (Fig. 21a). At 104500 d it exceeds β in the westward jet at the equator in the upper troposphere, and the westward equatorial jet becomes barotropically unstable (Fig. 21b). Given the turbulent fluctuations in the flow, the growth of the instability is then inevitable.
This change in the sign of u y / 2 2 is accompanied by a sharp increase in the eddy momentum flux convergence at the same locations. This deposits eastward AM at the equator and spins up the super-rotating equatorial jet. The jet remains barotropically unstable until 104800 d (Fig. 21c), but the strongly positive AM flux convergence persists at the equator (Fig. 21d), and by 105500 d a broad super-rotating equatorial jet has established itself at a speed it retains for the remainder of the simulation. The presence of barotropic instability at the same time as a sharp increase in the AM flux convergence indicates that the unstable jet converts zonal to eddy energy in the form of Rossby waves, which then deposit eastward momentum at the equator as they radiate away. In our simulations it is this mechanism, rather than convection, that drives the equatorial super-rotation. Fig. 22 instead shows the spin-up phase at the equator as time series averaged over the pressure and latitude points where the flow becomes barotropically unstable. Ignoring the first couple of 100 d after the change in resolution, when the flow is significantly perturbed, we see the point at which the sub-rotating jet becomes barotropically unstable at 104500 d. Immediately after there is a sharp increase (to 90%) of the fraction of KE held in eddies, implying that the barotropic instability in the mean flow is triggered. The absolute eddy KE increases thereafter, and continues to do so until the region becomes barotropically stable again.
Similarly, when the flow becomes barotropically unstable the eddy AM flux convergence jumps to significantly positive, pumping eastward momentum into the equatorial jet, which is reflected in the subsequent zonal KE and the zonal mean zonal velocity. This convergence is elevated during the period of instability and shortly after, between 104600-105400 d, and then falls back to a level close to zero but, on The zonal mean has been removed at each latitude. Horizontal lines are at the equator and at ± 1.05°. R.M.B. Young, et al. Icarus 326 (2019) 225-252 average, slightly positive. After this period of strong jet forcing the equatorial super-rotating jet reaches the speed it retains throughout the rest of the run. This behaviour at the equator points towards a mechanism involving a secondary barotropic instability on top of the primary baroclinic jets. Barotropic instability converts zonal to eddy KE, generating Rossby waves that radiate away from the equatorial region, depositing eastward AM at the equator.
This mechanism is essentially as described by Williams (2003a). He found that, when the baroclinic zone is sufficiently close to the equator, the equatorward side of baroclinic jets can become barotropically unstable, leading to a strongly super-rotating westerly at the equator. He found this was possible even for Earth conditions, although the faster the rotation, the closer to the equator the baroclinic zone needed to be for super-rotation (7°was sufficient for four times Earth's rotation). While Jupiter's rotation rate is only double that of the Earth, the relevant parameter is rather the thermal Rossby number, where a rotation rate of about eight times Earth's rotation would be equivalent, other things being equal . Based on extrapolating Williams (2003a, Fig. 2) a baroclinic jet at 3 5 would be required, which matches the latitude of the first baroclinic jet when it becomes barotropically unstable (Fig. 21b) reasonably well.

Maintenance of the prograde equatorial jet in Run B
In equilibrium the equatorial super-rotation is maintained by continuous equatorward migration of off-equatorial baroclinic jets, which keep the equatorial zone in a state of near-neutral stability with respect to barotropic instability. Fig. 23 shows ∂ 2 〈u〉/∂y 2 and β near the equator throughout Run B after the increase in resolution. Where the grey line is non-zero, at least one grid point within ± 1.76°(four grid points) of the equator in the upper troposphere breaks the barotropic stability criterion, i.e. ∂ 2 〈u〉/∂y 2 > β somewhere in this region. The initial spike where the flow becomes barotropically unstable during spin-up of the equatorial jet is clear. However, the flow remains close to marginal stability throughout the run, and at various points later on the barotropic stability criterion is broken again.
By comparing this time series with the Hövmoller diagram in Fig. 13b, times when the equatorial flow becomes barotropically unstable again generally coincide with times when off-equatorial jets merge with the equatorial jet. This happens between 121000-126000 d, as a jet merges from the south, and at the very end of the run when another jet merges from the south. This repeated breaking of the barotropic stability criterion implies that once the super-rotating equatorial jet is set up, it is maintained by undergoing further episodes of barotropic instability as off-equatorial jets merge with it. This maintains a source of Rossby waves and hence eastward AM at the equator, maintaining the super-rotation. The time series goes in cycles, becoming barotropically unstable for a short time before fall back to stable conditions and then slowly building up to instability again. In the Hövmoller diagram in Fig. 13b, after 110000 d the equatorial jet takes on a slightly cusped form, which maintains ∂ 2 〈u〉/∂y 2 at the equator closer to β than if the flow at the equator were peaked. Such a peak occurs on the real planet, although the jet is much wider.
As at midlatitudes (Sect. 4.4) we can make use of the Eulerian mean equation (Eq. 9) to diagnose the dominant balance in the equatorial region for the flow in equilibrium. Fig. 24 shows the dominant terms in equilibrium within 30°of the equator as latitude-pressure cross-sections, averaged over a 100 d period finishing at 133000 d, and Fig. 25 shows a latitudinal profile of all the terms at 0.4 bar. Supplementary Figure S6 shows all the terms in the Eulerian mean equation for this case, along with additional latitudinal profiles.
In the stratosphere, the balance is between the Coriolis acceleration ( Fig. 24a) and the meridional eddy momentum flux divergence (Fig. 24b). Within the sub-rotating equatorial jet the Coriolis term is positive and the eddy AM flux is divergent, peaking at the edge of the jet.
In the troposphere the balance is somewhat more complicated. Near the equator vertical diffusion (Fig. 24c) is weak as its timescale is f 1 . Therefore we do not need to appeal to the adjusted Coriolis acceleration in the equatorial balance. Right at the equator (in the cusped part of the jet), the meridional eddy AM flux convergence (Fig. 24b) and the vertical mean AM flux convergence (Fig. 24d) accelerate the eastward jet, and this is balanced by the vertical eddy AM flux divergence (Fig. 24e). The strong vertical eddy AM flux divergence in the equatorial jet in the upper troposphere indicates that as well as Rossby waves radiating away from the equator, there is also a vertical component to the wave action in the equatorial jet, with upwards-propagating waves generated at depth and breaking at the uppermost part of the eastward equatorial jet.
Just off the equator, where the equatorial jet maxima are, the main eastward acceleration comes from Coriolis acceleration (Fig. 24a). This where 〈U〉 is the zonal mean zonal velocity averaged over the four latitudes between ± 1.05°. There is no grid point in u exactly at the equator; the two latitudes closest to it ( ± 0.35°) have been averaged. R.M.B. Young, et al. Icarus 326 (2019) 225-252 is balanced by a combination of the meridional mean AM flux divergence (Fig. 24f) and the vertical eddy AM flux divergence (Fig. 24e). The first carries more zonal momentum poleward from the jet peaks than from the cusped part of the jet at the equator itself, and so decelerates the jet. The second is still relevant, but is not as strong as at the equator. There are also smaller decelerating contributions of O (10 m s ) 6 2 to this balance from the meridional eddy AM flux divergence (Fig. 24b), the vertical mean AM flux divergence (Fig. 24d), vertical diffusion (Fig. 24c), and the Shapiro filter (not shown), which acts to smooth out gradients and hence accelerates the cusped part of the jet and decelerates the jet peaks.
Even averaging over a short 100 d period, the net acceleration due  Table 3 Properties of equatorial Rossby waves identified from the power spectra for Run B in Fig. 20b with symmetric meridional structure ( = n 1). All values are to two significant figures.  Young, et al. Icarus 326 (2019) 225-252 to the sum of all the terms in the Eulerian mean equation is typically (2-4) × 10 m s , 7 2 which is significantly smaller than the dominant individual terms in Fig. 24 and 25, indicating that the jets are well balanced. However, the net acceleration in the equatorial jet fluctuates about zero with a variability significantly larger than the time mean  The black line shows the maximum value of ∂ 2 〈u〉/∂y 2 between 0.3-3 bar and ± 1.76°latitude. The dotted line shows β at the same location (variation over these latitudes indistinguishable by eye). The grey line shows the percentage of grid points within that region that break the barotropic stability criterion, i.e. for which ∂ 2 〈u〉/∂y 2 > β. R.M.B. Young, et al. Icarus 326 (2019) 225-252 acceleration. Fig. 26 shows the net acceleration averaged over 4000 d. We can estimate a "spin-down" timescale associated with the zonal mean flow: corresponding to a spin-down time of 4000-10000 d, depending on pressure. This timescale agrees well with the rate at which changes occur due to jet migration (Fig. 13b), and in Fig. 26a the net acceleration in the offequatorial jets clearly reflects their equatorward migration over this period.
During spin-up of the equatorial jet in Run B most of the terms are similar in form to the equilibrium balance, although the magnitudes of the Coriolis and meridional mean AM flux divergence are larger. The net imbalance at the equator during spin-up is dominated by a particularly strong and deep region of meridional eddy AM flux convergence. The associated net acceleration is at least an order of magnitude larger than the net acceleration in equilibrium. During the spin-up τ spin falls below 10 d in this region, consistent with the time spent for the equatorial jet to spin up. Supplementary Figure S7 shows all the terms in the Eulerian mean equation at one point during the spin-up.
To identify the length scales that contribute to the AM and heat fluxes relevant to the super-rotating equatorial jet, we computed the zonal co-spectrum of the eddy fluxes (Saravanan, 1993). This splits the latitudinal profiles of eddy AM flux and heat flux into their zonal spectral components, and hence identifies which length scales are responsible for the meridional fluxes. The eddy AM flux co-spectrum is and for eddy heat flux it is in terms of zonal wavenumber m and latitude ϕ, where u m, is the (complex) zonal spectrum of the eddy zonal velocity u′ (similarly for v′ and θ′).
We computed the co-spectra of these two quantities and plot them for Run B in equilibrium at 0.4 bar in Fig. 27. Supplementary Figure S10 shows additional levels. In the AM co-spectrum (Fig. 27a), most of the flux is due to wavenumbers 40-80, except right at the equator, consistent with wavelengths around the deformation scale. At   R.M.B. Young, et al. Icarus 326 (2019) 225-252 the equator there is instead power in all modes up to 80. There is divergence of AM flux at low wavenumbers, but convergence (which wins overall) at higher wavenumbers. This is the same as seen by Saravanan (1993) in his case S (superrotating) case. His super-rotating jet was similar to that observed here, with a slightly cusped structure at the equator, but his planet was Earth so the jet itself was wider (some 30°latitude). He found a broad range of higher wavenumbers that converged AM into the equatorial jet, but a small range of low wavenumbers around the edge of the jet that diverged AM away from the equator. In our case the range of latitudes is small where the low-wavenumber flux is divergent. However, this is consistent with some AM being fluxed from the equator itself into the off-equatorial latitudes to create the cusped structure of the equatorial jet, while the majority of the (convergent) AM flux acts, over a wider range of latitudes, to maintain the main body of the super-rotating equatorial jet. For each off-equatorial jet we see corresponding latitudes of AM flux convergence and divergence between wavenumbers 40-80.
The eddy heat flux can be interpreted as a signature of baroclinic instability (Saravanan, 1993). The eddy heat flux co-spectrum (Fig. 27b) shows that most of the flux is poleward and in modes comparable with the relevant length scales for baroclinic instability. This confirms our earlier analysis that showed baroclinic instability is important for generating and maintaining the off-equatorial jets. We also note that in the cusped part of the equatorial jet right at the equator the heat flux at low wavenumbers is reversed, as was the AM flux.
In the stratosphere (well into the retrograde equatorial jet) the eddy heat flux is strongly poleward at almost all wavenumbers, compared with the upper troposphere. Conversely, the AM flux is weak but divergent at the equator at almost all wavenumbers, driving the retrograde equatorial jet in the stratosphere. At 4 bar, the AM flux does not have significant power at low wavenumbers at the equator. Above we linked the cusped jet at the equator to the reversal of equatorial AM and heat fluxes at low wavenumbers, and the flow at 4 bar is not cusped (Fig. 4b), consistent with this interpretation. Supplementary Figure S10 shows both these cases. Fig. 28 shows eddy AM flux co-spectra before and during equatorial jet spin-up. Before spin-up (Fig. 28a) it is noisy up to wavenumber 20, but this dissipates during spin-up. The alternating structure of converging and diverging latitudes at high wavenumbers develops before spinup. During spin-up, however (Fig. 28b), there is an important change at wavenumbers 30-80 near the equator. There is additional AM flux convergence at high wavenumbers within a few degrees of the equator that was not there before spin-up. This is consistent with a new source of Rossby waves that was not present before spin-up, which then spin up the super-rotating equatorial jet. This part of the co-spectrum maintains its power while the flow is in equilibrium (Fig. 27b), showing    R.M.B. Young, et al. Icarus 326 (2019) 225-252 that the same mechanism that spun up the super-rotating equatorial jet also maintains it.

Spin-up and maintenance of the sub-rotating equatorial jet in Run A
The sub-rotating equatorial jet spins-up very quickly in Run A. The main spin-up phase occurs during the first few thousand days (Fig. 13a), and in just 5000-10000 d it reaches the width it retains for the remainder of the run. Doubling the horizontal resolution does not affect the equatorial flow. Fig. 29 shows the dominant terms in the Eulerian mean equation during spin-up. There is a strong mean meridional circulation confined to the upper troposphere, but this is balanced almost completely by vertical diffusion, so is shown as the adjusted Coriolis acceleration (Fig. 29a). The residual is then balanced by the meridional mean AM flux convergence (Fig. 29b). The retrograde equatorial jet is spun up by the meridional eddy AM flux divergence (Fig. 29c). This is strongest at the latitudinal edge of the jet, and only partially balanced by the adjusted Coriolis acceleration there. The other terms in the Eulerian mean equation are small.
The net acceleration (Fig. 30) shows the edges of the jet are accelerated westward, so over time the retrograde jet widens, while within the equatorial jet the Coriolis term stabilises the jet speed by accelerating the jet eastward. Around 2000 d these balance and jet widening ceases. The spin-down time during spin-up (Fig. 30b) shows 10-100 d timescales at the edge of the equatorial jet, consistent with a widening jet over a short period, and a very small net acceleration inside the jet (6000-12000 d timescale).
In equilibrium, the dominant balance in the upper troposphere remains the same. Elsewhere, where vertical diffusion is weak, the Coriolis term is balanced by meridional eddy AM flux divergence. Within the equatorial jet the meridional eddy AM flux is almost everywhere divergent, acting to maintain the sub-rotation, balanced by a very weak meridional circulation via the Coriolis term. This balance is shown in Fig. 31. The net acceleration within the equatorial jet is extremely weak, such that the typical spin-down time within the jet is 0.5-2 million days, significantly longer than the length of the simulation. The momentum budget in Run A is very finely balanced, and hence the zonal jets will be extremely long-lived.
The co-spectra for Run A (shown in Supplementary Figure S11) show that during spin-up the eddy AM flux is poleward within 40°of the equator up to wavenumber 60, with the strongest flux divergence at the edge of the jet. The associated eddy heat flux is consistent with midlatitude baroclinic instability. In equilibrium the peaks are more restricted in wavenumber space than during spin-up, as in Run B. The peaks in latitude and zonal wavenumber indicate that baroclinic instability remains the dominant process driving the circulation in equilibrium.

Conclusion
In this paper we have presented a major update to the Jupiter model initially developed by Yamazaki et al. (2004) and Zuchowski et al. (2009c), and used it to investigate the spin-up of Jupiter's jets with and without an interior heat flux. Other investigations making use of the updated model will be presented elsewhere, including cloud dynamics in Part II, studying the effects of moist convection on the atmosphere, and analysing the model's complex energy cycle in the context of 2D and quasi-geostrophic turbulence.
The most important updates were to convert the model domain from a limited area of the southern hemisphere to a global domain, to add a more realistic two-stream radiative transfer scheme, which represents realistic solar forcing at the top of the atmosphere and interior heat flux from below, and to redesign the model using the MITgcm rather than the Met-Office ExtUM. We have refined other physical parametrisations, in particular vertical diffusion and dry convective adjustment, to use more physically-justified timescales.
In both simulations with and without an interior heat flux a global  Figure S9 shows all the terms. Fig. 26 but for Run A during spin-up, averaged over 1300-1500 d. R.M.B. Young, et al. Icarus 326 (2019) 225-252 banded jet structure emerged, with midlatitude jets driven by baroclinic instability. The response of the atmosphere to heating from below was clear. Without heating there was a broad, sub-rotating, equatorial jet with weak eddy activity at low latitudes and a strong equator-to-pole temperature gradient. With more realistic 5.7 W m −2 heating from below, the equatorial jet became super-rotating, with a significant reduction in its width, a large increase in eddy activity, and a reduction in the equator-pole temperature gradient. The strength of the super-rotating jet, however, remained somewhat less than the observed 100 m s 1 . As the only difference between the two runs, the interior heating must ultimately cause the super-rotation. The following mechanism is proposed. At midlatitudes the flow is baroclinically unstable (Fig. 15,  27b). These instabilities generate Rossby waves that propagate to other latitudes, depositing westward momentum where they break. Baroclinic instability disappears at the equator, as f goes to zero there, so there is a net deposition of westward angular momentum in the equatorial region, which drives a sub-rotating jet by eddy momentum flux divergence (Fig. 24b, 29b). In the absence of other processes, a subrotating equatorial jet develops, as in Run A and in the Run B stratosphere.

Fig. 30. As
When there is interior heating, tropospheric temperatures become neutrally stable with respect to the dry adiabat (Fig. 5). Hence static stability increases towards the equator (Fig. 16) and a poleward bias in the Eady growth rate across the midlatitude jets results (Fig. 15). Via the mechanism in Chan et al. (2007) this creates an asymmetry in eddy angular momentum flux convergence across the jets (Fig. 12d,h), driving the jets equatorward (Fig. 13b). As the eastward zonal jets approach the equator they become barotropically unstable (Fig. 21), generating Rossby waves there (Fig. 20). These waves propagate away from the equator, depositing eastward angular momentum at their source and spinning up the tropospheric super-rotating equatorial jet in Run B.
While a migratory mechanism for equatorial super-rotation is invoked here, we should note that Jupiter's jets do not actually appear to migrate. Our migration timescale is 50-60 years per jet width, and over the 21 year period between Voyager (Limaye, 1986) and Cassini (Porco et al., 2003) the jets do not appear to have moved much at all, except perhaps a 0.5-1.0°equatorward migration poleward of 30°S. Nevertheless, the main use of this kind of model is to explore processes that may be relevant under Jovian conditions, and our simulations demonstrate one possible mechanism for generating equatorial super-rotation in the giant planet context.
Since these simulations have been run, observations from the Juno spacecraft  are challenging many of the assumptions underlying our understanding of Jupiter's atmosphere . Gravity field measurements  should lead to an improved understanding of the deep wind structure , and observations of the poles reveal an environment very different from lower latitudes (Adriani et al., 2018). These insights will be important for informing future model development, in particular to put a physical justification on the velocity boundary condition at the bottom of the model domain, which has long been a problem with giant planet weather layer models.

Accompanying Data
Data from 154760-154860 d in Run A and 132900-133000 d in Run B can be obtained from the Oxford University Research Archive -Data (https://ora.ox.ac.uk) (Young et al., 2018a).   Figure S8 shows all the terms as latitude-pressure crosssections. R.M.B. Young, et al. Icarus 326 (2019) 225-252 National E-Infrastructure. The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work. The authors acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558.

Appendix A. Full model description
This appendix describes in full the various parametrisations that were added to the MITgcm to simulate Jupiter's atmosphere. We forked the MITgcm from the trunk at http://www.mitgcm.org at v1.2085 on 23 October 2014. Table A.1 lists Jupiter's physical and chemical characteristics used by the model.

A.2. Sponge layer
The sponge layer is a standard procedure to avoid gravity waves reflecting off the top of the domain. Counting from the top, the eddy components of the horizontal velocity and temperature at the uppermost three levels are damped with = 0.1, where u jk is the zonal mean zonal velocity at latitude ϕ j and pressure level p k .

A.3. Radiation scheme
At the top-of-atmosphere (TOA) the incoming shortwave radiation is set to the annual mean shortwave incoming solar radiation (ISR) as a function of latitude. The incoming longwave flux is assumed to be zero. For a non-synchronously rotating planet in a circular orbit, the daily average solar insolation at TOA on day j is where S 0 is the solar constant at Jupiter's orbit, A is the bond albedo, μ(ϕ, j) is the daily average of cos θ z , where θ z is the solar zenith angle, and H(ϕ, j) is the half-day length scaled by 2π. Following Vardavas and Taylor (2007, pp. 155-156), μ(ϕ, j) and H(ϕ, j) are  is the solar declination angle, ε is the planetary obliquity, and N year is the year length in local days. When the poles are illuminated continuously or not at all at certain times of year, = H j ( , ) 0 or π as required. The annual mean insolation at latitude ϕ is then year 1 year (A.6) accounting for partial days as necessary (Fig. 1). Radiative transfer within the atmosphere assumes the plane-parallel approximation and is resolved at each horizontal grid point by integrating the downwards SW and LW flux down from the top cell face (at vacuum) to the bottom of the model, and then the upwards LW flux back upwards. We use the two-stream approximation, i.e. that the LW and SW bands can be considered separately, as the black body curves for emission at the Sun's temperature reaching Jupiter and the blackbody curve for emission at Jupiter's effective temperature are separated in wavelength.
Beginning with the ISR as described above, we integrate downwards through the column, calculating the SW and LW fluxes at cell faces. Cell k is bounded by faces k below and + k 1 above. The SW flux is only attenuated, while the LW flux is both attenuated and added to by emission from the cell: represents the effect of integrating over a hemisphere under the plane parallel approximation (Ramanathan et al., 1985). We adopt the simple optical depth profiles used by Schneider and Liu (2009) where τ 0,S , τ 0,L , and p ref are reference values listed in Table A.1. The shortwave dependence on p simulates a well-mixed absorber, and the longwave dependence on p 2 simulates collision-induced absorption (Herzberg, 1952). The downward thermal emission from cell k is derived from Lacis and Oinas (1991, Appendix A) using an exponential dependence of the Planck function on optical depth (Fu and Liou, 1993): We then integrate the LW flux back upwards through the column, calculating the LW flux at each cell face. There is no upwards SW flux, which is assumed to have all been absorbed in the deep atmosphere, without scattering. As with the downwards LW flux, the upwards LW flux is added to by emission from the cell, and attenuated: In both the upwards and downwards paths we evaluate the thermal emission in two parts, above and below the cell pressure level p k . This is more accurate for atmospheres where there is a significant pressure (and hence optical depth) change over the cell. For downwards emission, the emission from the upper part of cell k is  R.M.B. Young, et al. Icarus 326 (2019) 225-252 and similarly for the lower part. At the lower cell face lower ,upper L,˜, (A.20) which combines emission from the lower part with emission from the upper part attenuated by the lower part. There is an analogous formula for upwards emission.
Once the upwards and downwards fluxes are calculated at each cell face, the heating rate within each cell is where ΔF k is the net flux into cell k, ΔA is the horizontal area of the cell face (constant with height for a shallow atmosphere), and m k is the mass of the cell. Assuming hydrostatic balance, the potential temperature tendency is