Cassini observations reveal a regime of zonostrophic macroturbulence on Jupiter

a College of Marine Science, University of South Florida, St. Petersburg, FL 33701, USA b Atmospheric, Oceanic & Planetary Physics, Clarendon Laboratory, Department of Physics, University of Oxford, UK c Department of Mechanical Engineering, Ben-Gurion University of the Negev, Beer-Sheva, Israel d Perlstone Center for Aeronautical Engineering Studies, Ben-Gurion University of the Negev, Beer-Sheva, Israel e John Adams Institute for Accelerator Science, Department of Physics, University of Oxford, UK f Department of Physics, University of Warwick, UK


a b s t r a c t
In December 2000, the Cassini fly-by near Jupiter delivered high-resolution images of Jupiter's clouds over the entire planet in a band between 50°N and 50°S. Three daily-averaged two-dimensional velocity snapshots extracted from these images are used to perform spectral analysis of jovian atmospheric macroturbulence. A similar analysis is also performed on alternative data documented by Choi and Showman (Choi, D., Showman, A. [2011]. Icarus 216, 597-609), based on a different method of image processing. The inter-comparison of the products of both analyses ensures a better constraint of the spectral estimates. Both analyses reveal strong anisotropy of the kinetic energy spectrum. The zonal spectrum is very steep and most of the kinetic energy resides in slowly evolving, alternating zonal (west-east) jets, while the non-zonal, or residual spectrum obeys the Kolmogorov-Kraichnan law specific to two-dimensional turbulence in the range of the inverse energy cascade. The spectral data is used to estimate the inverse cascade rate and the zonostrophy index R b for the first time. Although both datasets yield somewhat different values of , it is estimated to be in the range 0.5-1.0 Â 10 À5 m 2 s À3 . The ensuing values of R b J 5 belong well in the range of zonostrophic turbulence whose threshold corresponds to R b ' 2.5. We infer that the large-scale circulation is maintained by an anisotropic inverse energy cascade. The removal of the Great Red Spot from both datasets has no significant effect upon either the spectra or the inverse cascade rate. The spectral data are used to compute the rate of the energy exchange, W, between the non-zonal structures and the large-scale zonal flow. It is found that instantaneous values of W may exceed by an order of magnitude. Previous numerical simulations with a barotropic model suggest that W and attain comparable values only after averaging of W over a sufficiently long time.
Near-instantaneous values of W that have been routinely used to infer the rate of the kinetic energy supply to Jupiter's zonal flow may therefore significantly overestimate . This disparity between W and may resolve the long-standing conundrum of an unrealistically high rate of energy transfer to the zonal flow. The meridional diffusivity K / in the regime of zonostrophic turbulence is given by an expression that depends on . The value of K / estimated from the spectra is compared against data from the dispersion of stratospheric gases and debris resulting from the Shoemaker-Levy 9 comet and Wesley asteroid impacts in 1994 and 2009 respectively. Not only is K / found to be consistent with estimates for both impacts, but the eddy diffusivity found from observations appears to be scale-independent. This behaviour could be a consequence of the interaction between anisotropic turbulence and Rossby waves specific to the regime of zonostrophic macroturbulence.
Ó 2013 The Authors. Published by Elsevier Inc.

Introduction
Explaining the latitudinal dependence of zonal velocity profiles and, more generally, the large-scale dynamics of giant planet atmospheres are classical long-standing problems in geophysical fluid dynamics (Golitsyn, 1973;Beebe, 1994;Dowling, 1995). A wealth of observational information accumulated over the last 0019 several decades from NASA space missions and ground-based observations, combined with numerous theoretical and numerical studies, have led to progress in this area of science (Ingersoll et al., 2004;Vasavada and Showman, 2005;Del Genio et al., 2009). Some model studies, such as those by Heimpel and Aurnou (2007) and Schneider and Liu (2009), have yielded realistic zonal profiles. These profiles alone, however, provide only limited information on the physics of the general circulation. Assessing scaledependent flow energetics, energy transfers, and diffusion characteristics relies upon higher-order statistical moments, and the lack of observational data necessary to compute these moments has impaired our ability to understand turbulent processes on the giant planets despite ever-improving computing resources.
The Reynolds numbers characteristic of jovian atmospheric flows are very large and so these flows are highly turbulent. In this paper we focus our attention on various aspects of planetary turbulence, which will henceforth be referred to as macroturbulence (see, e.g., Held, 1999). The spatial and temporal distributions of turbulent energy sources and sinks are poorly known, but it can be expected that over certain ranges of length scales, flows develop spectral ranges characterized by a small number of ''local'' parameters while being largely independent of other details. The emergence of such ranges can be quantified in terms of properly chosen non-dimensional numbers. As an illustration we recall that inertial ranges emerge in both three-dimensional (3D) (e.g., Tennekes and Lumley, 1972) and two-dimensional (2D) (e.g., Lesieur, 1997) turbulence. In these two cases, the kinetic energy spectra obey the classical Kolmogorov and Kolmogorov-Kraichnan (KK hereafter) laws respectively.
Jovian macroturbulence exhibits some peculiarities. For instance, even though the high variability of small-scale flow features is well documented, large-scale structures, such as the zonal jets on all giant planets and the Great Red Spot on Jupiter, have preserved their main characteristics and appearance, apparently for centuries (Hooke, 1665). This seemingly paradoxical situation is typical of anisotropic turbulence with dispersive waves that facilitate the formation of slow manifolds with very long characteristic time scales (e.g., Vanneste, 2013). These slow manifolds emerge in the low-dimensional sub-spaces orthogonal to the directions in which the waves do not propagate (Sagaut and Cambon, 2008, p. 143).
Slow manifolds reveal their presence via the steepening of 1D kinetic energy spectra in the wave's non-propagating directions. For the jovian planets, the large-scale slow manifold can be associated with the zonal direction, which is orthogonal to the non-propagating (meridional) direction for barotropic Rossby-Haurwitz waves. The corresponding physical-space flow configuration is represented by a system of slow-varying alternating zonal jets with a steep meridional spectrum. Since this spectrum pertains to the zonal flows, we shall refer to it hereafter as the zonal spectrum. The spectrum formed by all other (non-zonal) modes will be referred to as the residual spectrum. Computer simulations indeed detected such configurations in barotropic 2D turbulence with a b-effect (meridional variation of the Coriolis parameter) both on a b-plane (Chekhlov et al., 1996) and on the surface of a rotating sphere . Galperin et al. (2006) coined the associated flow regime zonostrophic turbulence. Section 2 provides a brief survey of this regime.
The analysis by Galperin et al. (2006) suggests that due to a strong tendency towards barotropization, the barotropic modes of circulation on giant planets may contain significant fractions of the kinetic energy. In this case, as discussed by Del Genio et al. (2009), a barotropic theory of zonostrophic turbulence can be applied to explain and quantify some features of jovian macroturbulence.
Even in a simplified barotropic case, the presence of the zonostrophic regime can only be established by analyzing both zonal and residual spectra. Only the zonal and total kinetic energy spectra have been measured until recently (see Section 2). The Cassini spacecraft fly-by near Jupiter in 2000 provided a longawaited breakthrough for comprehensive studies of the spectral properties of Jupiter's atmosphere. It yielded two-dimensional cloud images with resolution and coverage sufficient to reconstruct daily-averaged 2D snapshots of the velocity field at Jupiter's cloud tops. These velocity fields can be used to compute fully-2D energy spectra and other diagnostics.
Accordingly, in the first part of this paper we describe our algorithm for evaluating fluid velocities from images of moving clouds. This is done by tracking clouds and other brightness features between images separated in time, under the assumption that their motions reflect the underlying wind velocities. For Jupiter's atmosphere, cloud tracking techniques have a long history going back to the Voyager missions. Early data processing methods required a human operator to identify individual features moving between images (Ingersoll et al., 1981). The utility of this approach was limited by the time necessary to inspect numerous cloud features and the size of features that could be reasonably resolved by the human eye. The seminal paper by Limaye (1986) introduced an automated method using 1D correlations along latitude circles, which was capable of retrieving zonal velocities at high latitudinal resolution. Through the 1990s and 2000s this approach was extended and refined to identify velocities in two dimensions using 2D correlation techniques (e.g., Read et al., 2005). Some of the most recent methods incorporate an advection equation into the procedure (Asay- Davis et al., 2009), or take a different approach using an 'optical flow' technique (Liu et al., 2012). Asay-Davis et al. (2009) provide a good review of the methods used in the past. Similar methods have also been applied to other planets such as Venus  and Saturn (Sánchez-Lavega et al., 2003). Cloud images have also been used to estimate winds on Earth (Héas et al., 2012).
We describe our method in detail in Section 3 along with the observations from the Cassini flyby at closest approach to the planet in December 2000, while Section 4 describes the basic velocity fields and associated products. This dataset has been analyzed by a few groups previously, with the interest placed either on zonal diagnostics (Salyk et al., 2006) or the total kinetic energy spectrum in various image filters (Choi and Showman, 2011). Our method complements these other groups' approaches, and the visual comparison of our fields with these earlier results is encouraging. There is an extended error analysis associated with the cloud tracking section, which is described in the Appendix.
The second part of the paper describes our use of the ensuing velocity fields to analyze spectral properties of Jupiter's macroturbulence. First, in Section 5 we deduce the probable presence of the zonostrophic regime. We show that the residual spectrum possesses a relatively robust inertial range consistent with an inverse energy cascade. For the first time, the rate of this energy cascade, , is estimated from the data. We then evaluate the zonostrophy index, R b , as defined in Section 2. The estimated value of R b places Jupiter's macroturbulence squarely in the zonostrophic regime.
This regime has specific physical-space manifestations. First, there is a strong disparity between the energies of zonal jets and vortical and wave motions, which can be quantified by estimating the zonal and non-zonal spectra. A second pertains to the energy exchange between zonal flows and eddies. In Section 6 we show that the energy exchange between the large-scale vortices and zonal jets fluctuates in time, and may sporadically exceed by an order of magnitude. Finally, a third manifestation is traced to the laws of meridional diffusion, as the meridional diffusivity in this regime becomes scale-independent and is determined solely by and b. Section 7 presents an analysis of the meridional diffusion using data from the Shoemaker-Levy 9 comet and Wesley asteroid impacts. Section 8 provides concluding remarks.

The regime of zonostrophic turbulence
As mentioned earlier, computer simulations of 2D barotropic turbulence on a b-plane and on the surface of a rotating sphere gave rise to a classification of flow regimes possible in these systems . One of these is the regime of zonostrophic turbulence, which is distinguished by an anisotropic inverse energy cascade, strongly anisotropic spectra, and slowly varying systems of strong alternating zonal jets. Invoking the barotropic framework is appropriate for flows with a small Burger number, Bu = (L d /L) 2 , where L d is the first baroclinic Rossby radius and L is the maximum length scale in the system (Galperin et al., 2006). This condition stipulates funneling of the kinetic energy of large-scale circulations into the barotropic mode (Charney, 1971;Rhines, 1977;Salmon, 1998;Vallis, 2006). The regime of zonostrophic turbulence develops in the barotropic mode if the zonostrophic inertial range is sufficiently wide .
The notion of a zonostrophic inertial range extends Kraichnan's classical definition of the energy transfer range (ER) (Kraichnan, 1971) to 2D turbulent flows with a b-effect. For homogeneous, isotropic, 2D turbulence in a statistically steady state on the surface of a non-rotating sphere, the ER is confined between the wave numbers of the small-scale forcing, n n , and the large-scale friction, n fr , and features an isotropic inverse energy cascade with a constant rate . The energy spectrum in this range attains the celebrated KK form (e.g., Huang et al., 2001).
On a rotating sphere, a b-effect gives rise to a transitional length scale, L b = p (/b 3 ) 1/5 , at which the turbulent and Rossby wave frequencies are equal (Pelinovsky, 1978;Vallis and Maltrud, 1993;Huang et al., 2001).
and R are the angular velocity and the radius of the sphere, respectively. This length scale is related to the wave number 1=5 which has the units of inverse length.
When spherical harmonic decomposition is invoked, the basis functions, Y m n ðk; /Þ; k and / being the longitude and the latitude, respectively, depend on the nondimensional total and zonal indices, n and m, respectively (see, e.g., Boer, 1983;Boer and Shepherd, 1983;Burgess et al., 2013). Since / 2 [0, p], a length scale L corresponding to the total index n is given by n = pR/L. Clearly, the wave numbern and the total index n are related by n ¼ Rn. Since this relationship is quite straightforward, we shall not differentiate between the wave numbers and indices henceforth and omit the hat in the notations for the wave numbers.
The smallness of Jupiter's forcing scale is supported by the relatively small size of observed convective structures (Ingersoll et al., 2000;Li et al., 2006), and the magnitude of the first baroclinic Rossby radius Achterberg and Ingersoll, 1989), either (or both) of which characterize possible energy sources for the large-scale circulation. That is constant is of course an idealization, but is plausible in the statistical sense. If n b < n n , then n b splits the ER into two subranges, (n b , n n ) and (n fr , n b ) (Galperin and Sukoriansky, 2010, Fig. 5). In the former, the influence of a b-effect (and, hence, the dependence on b), is weak and the flow behaves largely like classical non-rotating KK (2D) turbulence. In the latter, which is the zonostrophic inertial range, the b-effect dominates and the flow undergoes increasing anisotropization .
Following Boer (1983) and Boer and Shepherd (1983), the energy spectrum in spherical polar coordinates is given by where n and m are the total and zonal indices, respectively, the modal spectrum Eðn; mÞ is the spectral energy density per mode (n, m), and Dm = 1 preserves correct dimensions when the wave numbers are used instead of the spectral indices, which indeed will be the case in most of the presentation henceforth. The spectrum E(n) can be decomposed into a sum of zonal and nonzonal (or residual) components, E(n) = E Z (n) + E R (n), where the zonal spectrum can also be written as E Z ðnÞ ¼ Eðn; 0ÞDm.
In the regime of zonostrophic turbulence, the zonal and the residual spectra are Sukoriansky et al., 2002) i.e., the residual spectrum preserves the KK form . Using these expressions, one can estimate the value of n at which E Z (n) and E R (n) intersect, which is the transitional wavenumber mentioned above: In their simulations Galperin and Sukoriansky (2008) and  found that (C Z /C K ) 3/10 ' 0.5. For large n, the residual spectrum is expected to be isotropic. Hence, for each m the modal residual spectrum can be written as Eðn; mÞ % E R ðnÞ=ð2nÞ ¼ ð1=2ÞC K 2=3 n À8=3 . Equating this to Eq. (2b), the transition between À5 and À8/3 scaling exponents in the modal spectrum takes place at Both Eqs. (3) and (4) will be used later as benchmarks in the data analysis.
The steep zonal spectrum given by Eq. (2a) sets zonostrophic turbulence apart from conventional 2D turbulence, where energy moves to large scales while enstrophy cascades to small scales. According to Eq. (2a), both the kinetic energy and the enstrophy, n 2 E Z (n) / n À3 , accumulate on large scales. A recent study by Venaille et al. (2012) demonstrated that a b-effect increases layer-wise enstrophy and enhances the effectiveness of the barotropization process. This finding highlights another aspect of the symbiosis between the barotropization and a strong b-effect.
An additional characteristic scale of 2D turbulence on a rotating sphere is defined by the Rhines wavenumber, n R = (b/2U) 1/2 , where U is the RMS velocity (Hide, 1966;Rhines, 1975). In flows with a strong b-effect, n R ' n fr (Sukoriansky et al., 2007). The importance of the b-effect can be measured using the zonostrophy index, R b = n b /n R , which characterizes the width of the zonostrophic inertial range (Sukoriansky et al., 2007;. Zonostrophic turbulence emerges for R b J 2.5 (Galperin et al., 2006;Galperin and Sukoriansky, 2008;Sukoriansky et al., 2007).
The spectra defined by Eqs. (2a) and (2b) are verifiable predictions. The zonal spectrum (Eq. (2a)) has been compared with available data for the jovian planets previously, and found to be in good agreement Sukoriansky et al., 2002;Del Genio et al., 2009;Barrado-Izagirre et al., 2009). The length of the observational record was clearly insufficient to compute meaningful statistical averages, but the agreement between theory and observations for both the spectral slope and the amplitude has been encouraging.
The data from the 2000 Cassini fly-by made it possible to test both the zonal and the residual spectra against observations. Choi and Showman (2011, hereafter CS) showed nearly identical zonal and total spectra on large scales. Further analysis of the spectra is one of the main purposes of this paper. In particular, if the measured residual spectrum agrees with (2b), then one could diagnose values of and R b directly from the observational data. A large R b (2.5 or larger) would point to the regime of zonostrophic macroturbulence. This regime has already been hinted at by the strong flow anisotropy observed in Jupiter's atmosphere.
Meridional diffusion of a conservative tracer in the regime of zonostrophic turbulence was investigated by Sukoriansky et al. (2009) and . They found that the meridional diffusivity K / in this regime becomes scale-independent and is determined solely by and b; see Eq. (18) below. The evaluation of using Eq. (2b) allows one to compute the magnitude of K / and compare it with observational data. Such a comparison could provide an independent test for the presence of the zonostrophic regime. The scale independence of Eq. (18) is an important consequence of the different roles played by turbulence and Rossby waves in transport processes (Sukoriansky et al., 2009;. If confirmed, Eq. (18) would provide a new insight into the physics of the general circulation and meridional diffusivity on Jupiter.

Cassini observations and cloud tracking
In this section we present the method we used to extract velocity fields from observations taken during the Cassini flyby of Jupiter.

Cassini observations
In December 2000 the NASA spacecraft Cassini passed Jupiter, and among other observations imaged the northern and southern hemispheres of the planet using the narrow-angle camera within the Imaging Science Subsystem (ISS) instrument. A small subset (304) of the many raw images taken at its closest approach to the planet (at around 140 Jupiter radii) have been pre-processed, projected onto regularly-spaced System III west longitude (k)/planetocentric latitude (/) coordinates centred over the equator, and made publicly available for analysis via the NASA Planetary Data System (Vasavada et al., 2008). The images in this dataset span almost four complete Jupiter rotation periods, beginning at 13:32:01 UTC on 11 December 2000 and ending at 04:37:06 UTC on 13 December 2000. Images were taken approximately every 63 min as the planet rotated below the spacecraft.
Images are available using four filters: a blue (BL1) filter at 455 nm, a near-IR (CB2) continuum band at 750 nm, and two methane bands (MT2 and MT3) at 727 nm and 889 nm. To calculate velocities from this dataset we selected images in the CB2 filter. This filter has the highest native resolution at 0.05°pixel À1 , compared with 0.1°pixel À1 for the others. This allows us to compute kinetic energy spectra from the velocity fields up to the highest possible wavenumber, and because these images show light reflected from Jupiter's cloud tops, they show the clouds' visible appearance best.

Cloud tracking
We calculated the horizontal wind fields u = (u, v) using an automated ''cloud tracking'' method based on Correlation Imaging Velocimetry (CIV) (Fincham and Spedding, 1997;Fincham and Delerce, 2000); here, u and v are the zonal (eastward) and the meridional (northward) velocities, respectively. The approach is similar to the methods used in the work cited in the introduction. The particular software we used was initially developed to analyze laboratory data from the 13 m rotating tank platform at LEGI, Grenoble (Coriolis Platform/ LEGI, 2013).
The principle of the cloud tracking method is to take a pair of images of the same place but separated in time, and correlate the brightness of pixel patches as they move from the first image to the second. At wavelengths where the brightness is due to reflection from clouds, moving patches correspond to cloud movement from which the underlying wind velocity can be estimated. The algorithm works on the assumption that brightness features will retain their general appearance as they evolve over the time between images. This puts a constraint on the time between images in a pair: too short, and the error in the estimated velocities will be comparable to the velocities themselves; too long, and the cloud features will deform too much to be identifiable in the second image.
From the CB2 filter data we identified all image pairs separated by 63 min, for northern and southern hemispheres separately. Each image pair was aligned using the sub-spacecraft longitude contained within the image metadata, and each image in the pair was then cropped to leave only the overlapping pixels. Early tests showed it was not possible to obtain good velocities near the poles, due to the low contrast and deformation of the images onto the latitude-longitude grid at those latitudes, so we cropped all images from 54°to the pole in each hemisphere. An example cropped image pair is shown in Fig. 1. After this procedure we were left with 70 aligned image pairs, 35 in each hemisphere, which are listed in the Supplementary Material. Initially we had used images separated by one jovian rotation period ($10 h), but this did not work because the cloud features had deformed too much over that time to extract good velocities. It also required an unacceptable amount of human intervention because latitudes with strong zonal jets had to be processed separately using different tracking parameters.
For each image pair the CIV algorithm contained two stages. In the first stage, ''CIV1'' (Fincham and Spedding, 1997), a first estimate of the wind vectors was made based only on translation of brightness features. A square pixel patch or ''correlation box'' in the first image with 23 pixels to a side was compared with equivalent patches within a larger ''search box'' (75 Â 60 pixels) centred on the same point in the second image. The conversion from pixels to distance on the planet is 1 pixel 0.05° 61 cos / km using R J = 69,911 km (NASA, 2010), where / is the planetocentric latitude, and the factor cos / applies to the longitudinal distance only. Hence the correlation box had side 1.15°at the equator, and the search box had sides 3.75 Â 3°there. The translation of the correlation box that gave the highest 2D correlation coefficient between the brightness of the matched patches in the two images defined the velocity of that patch (with the position of the velocity vector being the mid-point of the displacement). A fitting and interpolation procedure to the location with maximum correlation coefficient allowed sub-pixel accuracy. This procedure was repeated for all correlation boxes spaced by 15 pixels (0.75°at the equator) in both longitude and latitude. The conversion from pixel displacement between images to velocities in m s À1 is 1 pixel 16.1 cos / m s À1 , with the cos / term present only in the zonal (u) velocity.
There are several ways to generate false positives using this method, so a number of automated filters were then run to remove erroneously-identified and poor quality vectors. Some of these filters were quite harsh, as we preferred to lose good vectors rather than keep bad ones. The thresholds were selected empirically through a visual inspection of the pre-and post-filtered velocity fields. First, all vectors with correlation coefficient below 0.2 were removed. Vectors where the correlation maximum was within two pixels of the search box edge were removed, as were vectors corresponding to correlation boxes with a low brightness variation, and those where the fit to the correlation function was unstable. To account for any local or global spurious values, we also removed vectors whose longitudinal displacement between images exceeded 20 pixels, or whose latitudinal displacement exceeded 10 pixels.
Vectors with speeds further than 1.5 standard deviations from the global mean velocity magnitude were removed (note ''global'' here means all vectors within the image pair, not over the whole planet). Finally, vectors with either velocity component further than 3.3 standard deviations from the local mean or median (local here meant within 15 pixels) were also removed. The global filter preferentially removed most vectors close to the planet's limb, where low contrast resulted in many spurious vectors. Finally, the vectors were interpolated to a regular grid and smoothed, using a thin plate spline (Coriolis Platform, 2011). The parameters for this step are listed in the Supplementary Material. In stage two, ''CIV2'', the correlation analysis was repeated, but now with rotation and deformation (shearing and stretching) also taken into account (Fincham and Delerce, 2000). The correlation box was reduced to 19 pixels (0.95°a t the equator) to capture smaller scales in the flow. This box size largely defines the ''true'' resolution of all the subsequent data products (see Section 3.4). The first stage velocities were used as an initial estimate, and the analysis was done within an area a few pixels around this, including sub-pixel translations. This procedure was repeated at points spaced by 10 pixels (0.5°at the equator) in longitude and latitude. Finally the same set of filters were applied, with the same parameters except the local filters used a 20 pixel averaging area, and the global filter removed vectors over three standard deviations from the mean velocity magnitude, rather than 1.5 standard deviations as was used in the first stage. Because in the next stage we need to combine the image pairs together, the field was not interpolated onto a regular grid here, as was done at the end of the first stage.
The optimum parameters for each step of this procedure were determined empirically using a few image pairs, selecting parameters primarily on the basis of a visual survey of the output wind field, as well as correlation statistics and the number of filtered vectors. The output was particularly sensitive to the correlation and search box sizes, with poor choices yielding obviously degraded velocity fields. Large vortices such as the Great Red Spot were particularly useful for determining the optimal parameters. This feature has a well-defined flow field that has been extensively studied using similar methods and observations in the past (e.g., Read et al., 2005;Choi et al., 2007), and hence when suboptimal parameters were used this was clear.
This procedure was run on all 70 image pairs, giving 70 vector velocity fields at different times and at various locations over the planet.

Post-processing of vector fields from individual image pairs
From the 70 vector fields corresponding to the different image pairs, we created one global mosaic for each jovian day. Table 1 lists the start and end times for the four days. Each image pair was assigned to the day corresponding to the timestamp of the first image in the pair. Most of the image pairs overlap, particularly near the equator, so we require a procedure to combine the vectors appropriately. Each image pair produces 60-70°of usable data in longitude. Before combining the fields, we removed vectors near the planet's limb where the original images are most distorted by the projection onto the latitude-longitude grid. For each pair this meant removing all vectors within 2.5°of the latitudinal edge nearest the equator, and at longitudes (for all latitudes) less than 5°f rom the limb at the equator. Second, we removed vectors at latitudes where the quality of the wind field suffers due to insufficient contrast in the raw images, which corresponded to all vectors poleward of ±50°latitude. These crops are fairly harsh, but we decided it was better to remove some good vectors rather than keep bad ones. Each image pair then contained vectors covering 60-70°in both longitude and latitude. We also removed vectors close to Io and its shadow. This left 1,123,505 vectors covering 360°in longitude and ±50°in latitude over the four days. We shall refer to this set as the ''filtered vectors''.
For each day we combined the filtered vectors from individual image pairs into mosaics on a regular 0.5°latitude-longitude grid, with grid points at (k i , / j ) and grid spacing dk in longitude and d/ in latitude (dk = d/ = 0.5°). We discuss the choice of 0.5°grid spacing in Appendix A.3. To combine the images, at each grid point we first identified all the vectors within x m° ( Fig. 2), setting x m equal to the grid spacing 0.5°. The velocity assigned at the grid point was the mean of these k vectors, with the jth vector weighted by the inverse of the product of the following three quantities: (a) dg j , the distance between the grid point and the jth velocity vector, (b) do j , the distance between the position of the jth velocity vector and the central overlap pixel for the image pair containing it, and (c) dt j , the total time difference between vector j and all the other vectors in the set:  Table 1 Times and dates corresponding to the start and end of each jovian day. We defined a day using the System III rotation period of 9.92425 h 9 h 55 min 27.3 s (Weiss, 2004, Table A2.3). The start of day one corresponds to the first image in the first image pair we processed (i.e., not the first image in the raw dataset).

Day
Start End If all the t i are the same (i.e., all vectors are from the same image pair) then dt j was set to one. The weight assigned to the jth vector is then Normalising each x j with W, the weighted velocity at the grid point is then Weighting (a) ensures that vectors closer to the grid point are assigned higher weight. Weighting (b) favours vectors closer to the centre of their image, i.e., those where the original image is least distorted. Finally, weighting (c) makes sure that vectors separated in time by a long interval from the other vectors in the set are assigned less weight. This is necessary because where the images have wrapped around the planet there may be overlap between two image pairs a full jovian rotation apart, and there the flow will have distorted considerably. This ensures that if there is a group of vectors, a few of which are a full rotation separated from the others, the vectors a full rotation out will be assigned less weight. This weighting is not particularly harsh, however. For example, if there are three vectors at time zero and one at time ten hours, they will be assigned relative time weights of 3:3:3:1.

Resolution of the mosaics
The true resolution of the final 0.5°grid mosaics was estimated from the CIV method and the mosaicing method. We estimated it by starting with two d-functions separated by x res . This was then convolved with a top-hat of width 0.05°, corresponding to the pixel width, followed by a top-hat of width 0.95°, corresponding to the width of the CIV correlation function (at the equator), and finally with 1/x, representing the function used to weight the filtered vectors by their distance from the grid point in the mosaicing step.
Because the convolution with the 1/x function is so steep, the separation x res at which the two initial d-functions can be resolved is not very sensitive to the resolution criterion used. The CIV correlation box width is the main determinant of the true resolution. For example, using Rayleigh's criterion, 1 i.e., 26.3% contrast along the line joining the two points to be resolved, the separation required is only 0.97°. Even if the criterion is 50% contrast the separation need only be 0.99°. So given this is an estimate, the true resolution is around 1°. Fig. 3 shows the velocity, vorticity, total and eddy kinetic energy fields. Day 2 is shown as it contained the most complete coverage. Equivalent figures for the other three days are shown in the Supplementary Material.

Fields produced by cloud tracking
The velocity field ( Fig. 3(a)) shows that the main features of Jupiter's atmosphere are successfully identified using our method. The Great Red Spot and Oval BA are obvious, as are the major jet features. These features are even clearer in the vorticity field ( Fig. 3(b)), as we can pick out much finer detail because vorticity is a differentiated quantity. The band of anticyclonic vorticity around the Great Red Spot is picked out, as is the more quiescent region at the centre of the Spot, and the meandering westward jet to the north. The zonal jets are easily identified in the vorticity field by virtue of the strong shear they impart on the flow, which appears as regions of strongly negative or positive vorticity. Many of the smaller-scale features are clearly real as well. This can be seen better by cycling through the vorticity fields for the four days. For example, two small cyclonic features around 265°W 40°N and 264°W 46°N appear in all four days, at very nearly the same position and with very similar magnitude, showing the method can identify persistent features with angular diameter less than 1.5°across in this case, just larger than the resolution limit. It is also possible to pick out smallscale features moving with the mean flow, for example in the jet at 15-20°N, or the feature around 150°W 5°N, which is advected eastward.
There are some artifacts of the mosaic processing, however, and probably some from the CIV step as well. In the day 2 vorticity field the clearest artifact is at 260°W 5°N (notice how obvious this is in the eddy kinetic energy field in particular). Such artifacts tend to occur primarily at the joins between image pairs, and particularly in the equatorial region, where image pairs from both hemispheres are combined together as well as pairs separated by an hour or two. The point density map (Supplementary Material,Fig. S5b) shows how many filtered vectors are within each dk by d/ grid box before the mosaic step. Near the equator the number of points is larger, and in particular the region from 200°to 240°W is the longitude band where the earliest image pairs for this day overlap with the latest image pairs in the day, some 10 h later. Hence in this region there will be grid boxes containing vectors separated by 10 h, so perhaps it is not surprising that we find some artifacts in this region, as the flow will have evolved considerably over that time.
The random error in individual filtered velocity vectors associated with the CIV procedure is around 7.4 m s À1 , and in general the error in the mosaiced velocities is greater than the error in individual filtered vectors. Appendix A.1 contains an error calculation for both the CIV and the mosaic steps, and the Supplementary Material (Figs. S6 and S7) contains maps of the estimated random error at each point in the mosaiced velocity field. While these errors are sometimes large, the median errors for day 2 are 11.7 m s À1 for u and 9.3 m s À1 for v, near the low end of the range. It is straightforward to see by eye that the regions where artifacts are present are exactly the same as those with the highest random error. This makes sense as the highest error will be in grid boxes with a large spread of velocity values within the box, and hence it may not be possible to define a representative velocity vector for that grid point. It is somewhat counterintuitive that the grid Fig. 2. Schematic of the grid used when combining image pairs into a mosaic. We averaged over a circle of radius x m degrees around the grid point, with vectors within the circle weighted as described in the text. Fig. 3. Velocity, relative vorticity, total kinetic energy, and eddy kinetic energy for day 2, where u is the zonal mean velocity, averaging over all the points in the mosaic at that latitude.

B. Galperin et al. / Icarus 229 (2014) 295-320
points surrounded by the largest number of points are most likely to lead to poor velocities in the mosaics. Fig. 2 of CS, which shows equivalent plots from their somewhat different cloud tracking analysis. The similarities between the two are encouraging, and give us confidence that both groups' cloud tracking methods are sound and that the features identified are real. We see a broadly similar range of values in both the total and eddy kinetic energy fields (our Fig. 3(c) and (d)). For example, apart from the obvious large jets and vortices, which are in general agreement, both analyses reveal strongly energetic features within the mean flow around 30°W, 140°W, 160°W, 225°W, and 290°W in the jet at 5°N, and a gap in the zonal jet around 220°W/5°S. Both analyses also pick out the stronger kinetic energy signature in the northern part of the Great Red Spot compared with the southern part, and the quiescent regions in the Great Red Spot, Oval BA, and the anticyclone at 175°W/37°N.

Fig. 3 should be compared with
The methodology used in CS to produce their velocity fields is discussed in full detail there and initially in Choi et al. (2007), but because we shall compare our results with theirs again later, we summarize it here for the reader. Their image selection is essentially the same as ours, with image pairs separated by 63 min. Before the tracking procedure they use a Minnaert correction to compensate for illumination differences between images, and apply a high-pass filter to enhance image contrast.
The first stage of their cloud tracking is very similar. They use a 2°correlation box to estimate the translation of pixel patches between images, the estimate being where the cross-correlation coefficient is maximized within the search area. A second stage then Fig. 4. Zonal mean diagnostics using data from all four days. Each profile uses all the filtered velocity vectors from all four days in latitude bins 0.5°wide, without any weighting, rather than taking the average of the mosaiced fields, so any errors introduced by the mosaicing step are avoided. refines this estimate using a 0.5°pixel patch to search around the first estimate. As far as we can tell, this second stage does not account for rotation and deformation of the patch. They repeat this analysis with correlation boxes centred on each pixel in the original images, resulting in a vector field at resolution 0.05°. It is unclear whether the additional resolution over and above the resolution determined by the correlation box size, i.e., 0.5°, is justified.
To combine the image pairs into global mosaics they removed vectors poleward of 50°, as we did, and combined overlapping image pairs by weighting individual vectors by their distance to the overlapping edges of the relevant images. Within 1°of the equator they did not use the tracked winds but interpolated across the equator using weighting nearest-neighbour averaging. They did not use a correlation coefficient cutoff during this procedure. Instead they removed bad vectors with a method similar to those used to remove cosmic ray noise, and also manually removed vectors identified as spurious by eye, filling in the gaps with nearest-neighbour averaging. Note that the start and end times for their jovian days (CS Table 1) are slightly different from ours.
Our method and that in CS generally complement each other, and it is encouraging that they converge on similar results. Our method accounts for rotation and deformation in the correlation analysis and includes a few more parameters in the weighting for the image-combining step. Theirs, on the other hand, relies on pre-processing the cloud brightnesses and manually removing spurious vectors. The automated filters employed in our approach did not remove the spurious vectors and simply left them where they were.
Both methods have their strengths and weaknesses in recovering details of the Jupiter's circulation and neither of them can be expected to perform better in all situations. Since the methods are mutually complementary, our strategy will be to consider them in tandem rather than in competition. Such an approach will better constrain the emerging flow parameters and distill the physical laws that govern Jupiter's atmospheric dynamics.

Zonal mean profiles
Fig. 4 displays the zonal mean profiles using the filtered (not mosaiced) data from all four days. Fig. 4(a) shows our zonal mean zonal velocity, plus the Porco et al. (2003) and Limaye (1986) profiles from the Cassini and Voyager missions respectively. Our profile generally agrees well with these earlier profiles, which calculated the zonal mean zonal velocity using the same principle for cloud tracking but using a 1D correlation technique that could only resolve the zonal velocities. There are a few differences of O (10 m s À1 ) at the peaks of the strongest eastward jets, which may be explained by the larger latitudinal averaging box that we used. The standard error in this profile is shown in Fig. 4(c), and is thinner than the solid line in Fig. 4(a).
The zonal mean meridional velocity profile is shown in Fig. 4(b); it oscillates about zero. The standard error in v is generally smaller than the standard error in u, but the percentage error will be larger because, in general, j uj ) j vj. The mean meridional velocity from this profile is 0.76 ± 0.09 m s À1 , which is nonzero within the error bars. Salyk et al. (2006) also found a nonzero v, although theirs was negative. This suggests some systematic error(s) are unaccounted for, but the lack of data poleward of ±50°may be to blame as well.
We include in Fig. 4(d) the eddy momentum flux u 0 v 0 and latitudinal gradient of the zonal mean zonal velocity d u=dy, primarily for comparison with previous results such as Salyk et al. (2006, Fig. 5). A positive correlation between the two indicates that the energy flux is, on average, from the eddies to the zonal flow. For this analysis the correlation coefficient R, between the two is between 0.4 and 0.5. Salyk et al. (2006) estimated their corresponding correlation coefficient at about 0.86. We find that our correlation is poorest near the equator. Fig. 5 is akin to Fig. 4(d) and presents the zonally-averaged, latitude-dependent correlation coefficient, and the meridional shear for both our and the CS data sets. Since the CS data has a smaller grid spacing, the ensuing c 12 (/) shows more detail. Although the behaviour of c 12 (/) is generally similar for both data sets, large differences are noticeable near the equator. Poor correlation with the mean meridional shear in this region for our data has already been mentioned earlier. The correlation looks no better for the CS data. One can also detect the prevalence of negative c 12 in the northern hemisphere and marginally positive c 12 in the southern hemisphere. The mean meridional shear, however, reveals no such tendency. The latitudinally averaged values of the correlation coefficient are estimated at À0.025 and À0.09 for our and CS data, respectively. When the latitudinal dependence is taken into account, the variation of c 12 (/) is between 0.4 and À0.5, typical of well developed turbulence. For comparison, recall that in the region where 'the law of the wall' holds in a classical 3D shear layer, the averaged value of c 12 is about À0.4 (Tennekes and Lumley, 1972).

Zonal and residual spectra
In this section we present the kinetic energy spectra derived from our and the CS datasets, and compare them against the model spectra defined by Eq. (2).

Computing the spectra
The 2D velocity fields were used to compute the zonal spectrum E Z (n) and non-zonal (residual) spectrum E R (n) using a well established procedure described in Boer (1983), Boer and Shepherd (1983), and Huang et al. (2001). We processed both our and the CS data in the same way.

Filling the data
The spectral analysis requires a value at each grid point over the globe, but there are gaps in our mosaics (the CS data was provided to us complete). There is no data poleward of ±50°, but there are also gaps where there are no filtered vectors within the required distance from the grid point. To fill in the spaces, for each missing point we used the equivalent data point in the mosaic from one of the other days. For day 1 the points in day 2 were used, for day 2 the points in day 3 were used, for day 3 day 4 was used, and for day 4 days 3 and 2 were used. At a few points where data was missing in both the filling and filled days, a random data point was selected from the same day and the same latitude, but from a different longitude. Experiments with the model problem showed this to be the most non-distorting strategy. Outside ±50°l atitude, the velocities were filled in with zeros. Studies with a model problem revealed that these high latitude points only affect the spectra at relatively low n, which is of no consequence to the residual spectrum whose Kolmogorov part extends to n > 100 as discussed below.
Only the first three days were used for the spectral analysis because there is a gap in the fourth day's data that significantly degrades the quality of the spectrum compared with the other days, even when filled using this method (see Supplementary Material).

Estimating the error bars
We used two independent methods to estimate the uncertainties in the spectra. They produced similar results, giving us confidence that our estimates are reasonable.
First, we estimated the errors using the standard deviation of the spectral amplitude at each wavenumber over the three days analyzed. Second, we used a bootstrap technique (Efron and Tibshirani, 1993, p. 99). Beginning with the filled velocity field for day 2, we reconstructed the velocity field by resampling the data longitudinally in blocks of pre-defined size (we tried 5°Â 5°, 10°Â 10°, and 20°Â 20°). We filled each latitude band with blocks from the same latitude band by randomly sampling blocks from that band centred on any longitudinal grid point, with replacement. Twenty-five resampled velocity fields were produced in this way, and the spectrum for each calculated, with the standard deviation of the power at each wavenumber over the 25 instances giving us an estimate of the 1r error in the spectrum.

Results
Fig. 6 summarizes the resulting spectra averaged over jovian days 1, 2 and 3, and compares our spectra obtained from the data with 0.5°resolution with their counterparts computed using the CS data. In order to account for the impact of the absence of data at high latitudes on spectral amplitudes, all the spectra were multiplied by a factor of 1.3, which is the ratio of the total surface area of a sphere to the area between ±50°latitude, where the data is available. The accuracy of this correction factor was tested in computer simulations with a model problem. Fig. 6(a) and (c) indicate that, on large scales, the zonal spectrum E Z (n) exceeds E R (n) by about two orders of magnitude such that the total (i.e., E(n)) and zonal spectra are practically indistinguishable. This result provides evidence for the strong anisotropy inherent to zonostrophic turbulence, and emphasizes the difference between the total 2D and 1D spectra characteristic of anisotropic quasi-2D turbulence with a slow manifold. As mentioned earlier, slow manifolds in flows with anisotropic dispersive waves are low-dimensional flow subsets in the directions normal to the directions of non-propagation of linear waves. For flows on a rotating sphere featuring Rossby-Haurwitz waves, the non-propagation direction corresponds to m = 0. Similar differences between the total 3D and 1D spectra, energy accumulation in slow manifolds, and steepening of the corresponding 1D spectra are discussed in e.g., Smith and Waleffe (2002), Sagaut and Cambon (2008), Sukoriansky et al. (2005), and  for flows with stable stratification and/or rotation. Zonal flows (m = 0) constitute a 1D slow manifold in turbulence with a b-effect. Retention of the kinetic energy in zonal flows and steepening of the zonal spectrum towards the slope given by Eq. (2a) are attributes of this slow manifold. The physics of flow anisotropization is reflected more clearly in 1D rather than 2D spectra. By emphasizing this aspect, we diverge from CS where the total 2D spectrum was emphasized. Fig. 6(b) demonstrates that the scaling of the smoothed zonal spectrum with Eq. (2a) is robust. Our value of the coefficient, C Z ' 2, differs from previous estimates. Although the slope of E Z (n) seems to become shallower for n J 200 in the CS data, our data preserves the n À5 slope for all n. There is further discussion of this spectrum below. Fig. 6(c) shows that our and CS's estimates of the residual spectrum are close to each other but differ quantitatively. To better quantify the differences, Fig. 6(d) presents the compensated residual spectrum C R (n) / E R (n)n 5/3 . In such a representation, the KK spectrum Eq. (2b) should appear as a horizontal straight line. This scaling emerges for n J 180 and 200 in our and CS's data, respectively.
By comparing ours and CS's compensated residual spectra, we notice that the smaller grid spacing CS spectrum also reports the smallest energy. This inconsistency presumably stems from the differences in the data processing. Both our and CS's data show well-defined KK ranges which, nevertheless, exhibit different energy levels. These levels will be used to compute high-and lowend estimates of the inverse energy cascade rate, , using Eq. (2b).  (b) show the data points used for the least-squares determination of the spectral slopes for our and CS's data, respectively. The large spread in the CS data is a consequence of large differences between their day 1 and days 2 and 3. (c) and (d) show the corresponding compensated spectra, C R (n) / E R (n)n 5/3 , scaled around unity, with ±r error bars. 5.3. Estimation of the rate of the inverse energy cascade, length scales, and the zonostrophy index Fig. 7 shows the data points used for the least-squares estimates of the spectral slopes as well as the three-day averaged compensated residual spectra computed from our and CS's data with corresponding ±r error bars.
The best fit power law, n a , to E R (n) using the least-squares estimate yields a = À1.84 with coefficient of determination R 2 = 0.644 in the range n 2 (180,235) for our data, and a = À1.72 with R 2 = 0.737 in the range n 2 (210,395) for the CS data. These values of a are within 10% and 3% of the theoretical power law n À5/3 in Eq.
(2b), respectively. Using Eq. (2b) and assuming that C K = 6, from Fig. 7, we estimate to be about 10 À5 m 2 s À3 and 0.5 Â 10 À5 m 2 s À3 for our and CS's data, respectively. These are the first ever estimates of inferred directly from Jupiter observations. Using Eq. (3) with b = 2.5 Â 10 À12 m À1 s À1 and C Z = 2, we can now compute n b ' 57 and 67 for our and CS's data. From Eq. (4) the respective wavenumbers of the intersection of the modal spectra are then n z = 430 and 550. The wavenumbers n b and n z are shown in Fig. 8 for both our and the CS datasets. Since the values of n z are large they are not fully resolved in either dataset, and even less so are the scales on which the residual spectrum is supposed to become isotropic.
Spectral estimates of the mean total, zonal, and residual kinetic energies per unit mass are b E tot ¼ 1236:4 m 2 s À2 ; b E Z ¼ 1148:6 m 2 s À2 , and b E R ¼ 87:8 m 2 s À2 , respectively. Using b E tot , the RMS velocity U and the Rhines' wavenumber are found to be about 50 m s À1 and 11, respectively. The resulting values of the zonostrophy index are R b ' 5.2 and 6.1 for our and CS's data.
In Section 2, we stated the criterion for the establishment of the zonostrophic turbulence regime is R b J 2.5. This criterion was developed from the analysis of numerical simulations where C Z in Eq. (2a) was around 0.5, and so the prefactor in Eq. (3) was (C Z /C K ) 3/10 ' 0.5. A more precise estimate of C Z based on Jupiter's zonal spectrum is 2, so the prefactor in Eq. (3) increases to 0.73. This increase in n b increases the criterion for defining the zonostrophic regime in terms of R b from 2.5 to 3.65. This value is only hypothetical, however, because it has not been substantiated in numerical simulations or in other data. In any case, the values of R b estimated from Jupiter observations using either our or CS processing significantly exceed this value, confirming that the jovian atmosphere indeed conforms to the regime of zonostrophic macroturbulence. The various diagnostic quantities calculated from our and the CS data are summarized in Table 2.  (3) and (4). From our measured values of b E tot , b E Z , and b E R we can calculate the zonal and non-zonal mean flow indices, zmf b E Z = b E tot and nzmf b E R = b E tot , akin to those introduced by Srinivasan and Young (2012); the values for our data are zmf % 0.93 and nzmf % 0.07. Furthermore, we can approximate the zonal and residual spectra shown in Fig. 8 by simple models consisting of the lines given by Eqs. (2a) and (2b) up to the frictional wavenumber n fr , and straight horizontal lines for n 6 n fr . By integrating we find analytical expressions for b E Z and b E R valid in the zonostrophic regime: Some implications of Eq. (9) were discussed by Galperin et al. (2001) and Sukoriansky et al. (2002). In this regime, n fr ' (10C Z ) 1/4 n R , hence the anisotropy index, i.e., it is a rapidly growing function of R b . Eqs. (9) and (10) can be used to derive an analytical dependence of zmf and nzmf on R b in the zonostrophic regime, zmf ¼ c=ð1 þ cÞ; nzmf ¼ 1=ð1 þ cÞ: For R b ' 5, Eqs. (12) and (13) yield zmf ' 0.93 and nzmf ' 0.07, in good agreement with the preceding estimates derived directly from the spectra. The thresholds of the zonostrophic regime, Fig. 9. Vorticity and velocities around the Great Red Spot from our CIV analysis. The top two rows show the data for day 3, and the bottom row shows velocity profiles for all four days. Day 3 is used because in the day 2 data (Fig. 3(b)) there is a probable data processing artifact to the south of the GRS, which spoils the velocity profiles (see (e)). The grey lines in (b)  In conclusion, we re-iterate that the regime of zonostrophic macroturbulence is well-established in the reported observations through multiple diagnostic attributes: spectra Eqs. (2a) and (2b), strong anisotropy measured in terms of c or zmf, and characteristic wavenumbers n b and n z . The physical backbone of this regime, the anisotropic inverse energy transfer from eddies to zonal flows, will be discussed in Section 6.

Sensitivity of the spectra to the Great Red Spot
The Great Red Spot is the largest and longest-lived vortex on Jupiter, and so it is important to understand its role in the dynamics and energetics of Jupiter's atmospheric circulation. Fig. 9 provides a closer look at this spectacular vortex. The top panels show the relative vorticity and velocity fields in the region containing the Spot from the day 3 mosaic, along with velocity profiles through the centre of the Spot, and the bottom panels show the equivalent velocity profiles for all four days.
The main features of the GRS are well known and have been widely discussed over the years (e.g., Golitsyn, 1973;Dowling and Ingersoll, 1989;Marcus, 1993). There is a high-speed anticyclonic collar of approximately 100-150 m s À1 W-E / 50-100 m s À1 N-S around a quiescent central region. There is a westward jet meandering around the north side of the Spot, and an eastward jet meandering around the south side, which contains within it a smaller cyclonic vortex of diameter approximately 7°/7500 km in longitude and 4°/5000 km in latitude. The velocities in the Spot are constant over the four days, within the error bars.
In the quiescent central region there is evidence for very weak cyclonic flow. The signal/noise ratio is small here, but the feature is consistently identified for all four days in the zonal velocity (Fig. 9, low left Fig. 19b). Liu et al. (2012) in particular examine this phenomenon in some detail. In almost all cases the velocity reversal is seen only in the zonal velocity, while the meridional velocity across the Spot is flat. The exception is Asay-Davis et al. (2009, Fig. 6a), who see a weak reversal in the meridional velocity in one of their Galileo datasets, but not in the other.
Here we are particularly interested in how the spectral characteristics of the flow depend on the presence of the GRS. In particular, how do the spectral moments change with the GRS removed? This will enable us to diagnose and quantify various aspects of the interaction.
We created a GRS-free dataset for day 2 of our data by replacing the data window containing the GRS by a randomly selected, structureless window of the same size located inside the same latitude band, which is shown in the top panel of Fig. 10. The zonal and residual spectra calculated for both the original and GRS-free versions of our day 2 data are presented in the bottom panels of Fig. 10. The GRS-free mean residual kinetic energy per unit mass is b E GRS ¼ 75:0 m 2 s À2 . Hence the total kinetic energy per unit mass of the GRS is b E R À b E GRS = 12.8 m 2 s À2 , based on the difference between these spectra. This is only about 1% of b E tot , but $15% of the total residual kinetic energy. On the other hand, Golitsyn (1973) estimated that nearly 90% of b E Z resides in the equatorial belt, which is therefore by far the most energetic structure in Jupiter's atmosphere.
As expected, since the GRS is a non-zonal structure and its spherical harmonic representation does not contain zonal modes, the zonal spectrum is practically unaffected by removing the GRS.
The GRS has no visible impact on the residual spectrum on scales smaller than n ' 30, or about 7000 km. These scales are in the inverse cascade range. Evidently, the GRS does not disturb the inverse cascade and may possibly absorb a fraction of the upscale energy flux. As could be expected, removing the GRS causes significant decrease of the eddy energy around n $ 10-20 as the Spot is by far the largest and most energetic non-zonal structure on these length scales.
The next section addresses the interactions and the energy exchange between the zonostrophic turbulence, the GRS, and other large eddies.

Kinetic energy transfer and energy conversion conundrums
To investigate the three-way interaction between the smallscale forcing, large-scale vortical structures, and zonal flow, we consider the kinetic energy transfer equation in the Lorenz (1955) and Peixoto and Oort (1992) energy cycle. In particular, we are interested in the term that represents the rate of transfer of eddy kinetic energy, K 0 , to the kinetic energy of the mean zonal flow, K. In Holton's (2004) notation, this term is where hi is a global average. The globally-averaged production term u 0 v 0 d u=dy gives the barotropic energy flux per unit mass from all non-zonal fluctuations to the zonal mean flow, which will be denoted by W. We discretize the spherical coordinates by using j = 1 Á Á Á N to denote the latitude bins. In this notation, W is Using Cassini data, Salyk et al. (2006) estimate W in the range (7.1-12.3) Â 10 À5 W kg À1 , and cite other values found in the literature, including from variants of their own analysis procedure, as being in the range (6-30) Â 10 À5 W kg À1 . Using some assumptions regarding the depth of mass distribution in Jupiter's atmosphere, they estimated that about 5% of Jupiter's total thermal energy release is converted to the kinetic energy of the zonal flows. Sánchez-Lavega (2011) estimates the conversion rate at 15% compared with only around 0.1% for the Earth and suggests that ''the kinetic energy conversion from eddies to the mean flow is a process much more efficient on the two giant planets (i.e., Jupiter and Saturn) than on Earth''. This large difference is difficult to explain within the framework of conventional thermodynamics and so it has long been considered an unresolved 'puzzle' (e.g., Salyk et al., 2006). Note that this conundrum bears on two aspects of the global energy budget: (1) conversion of thermal energy to the kinetic energy of eddies and (2) conversion of eddy kinetic energy (EKE) to the kinetic energy of the zonal flow. Conversion of thermal energy to EKE results from baroclinic instability or convection. From thermodynamic considerations, it is expected to be equally efficient (or, more precisely, inefficient) on all planets although, if understood in terms of the Carnot cycle, the efficiency should increase with decreasing temperature. The problem arises when the rate of the conversion (per unit mass in our analysis) associated with part (2) is equated to the instantaneous conversion rate, W. This obscures important physics.
Recall that the classical Lorenz kinetic energy transfer equation pertains to a statistically steady state. Along with the term (14), this equation involves other (baroclinic) terms that cannot be determined without detailed knowledge of the vertical structure of planetary circulation (Read, 1986;Irwin, 2009). However, as follows from the present study, there still remains an uncertainty regarding the interpretation of the instantaneous value of ½K 0 K or, equivalently, the conversion rate per unit mass, W, even in a barotropic framework (where the baroclinic terms disappear). Do they provide adequate characterizations of the mean energy conversion rate from EKE to the mean zonal flow? The Lorenz diagram takes no consideration of the nature of the flow dynamics and, more specifically, the direction of energy transfer between different length scales. Even in the case of purely 2D flows with inverse transfer, the rate of this transfer per unit mass is time-dependent (see, e.g., Sukoriansky et al., 2012). It characterizes the eddy-mean flow interactions and its mean rate is quantified by a fraction of . A single instantaneous parameter, W, may be time dependent and differ from by a large margin, in which case it would be inadequate to quantify the mean rate of the energy exchange. A warning in this direction for real flows comes from Section 5.3 in which we found that for the anisotropic turbulence of Jupiter, is in the range between 0.5 Â 10 À5 m 2 s À3 and 10 À5 m 2 s À3 , which is smaller than Salyk et al.'s (2006) estimate of W by an order of magnitude. In the light of this finding, the original conundrum can be re-formulated in terms of clarifying a relationship between and W. This step requires introducing a scale-dependent energy flux and will be addressed in the next subsection.
It can be argued that the conundrum only exists assuming that all jovian dynamics is constrained within a thin weather layer. Ingersoll and Pollard (1982) addressed this issue by assuming that jovian dynamics is driven by the mechanism of deep convection. Their analysis was based almost entirely on scaling in a Boussinesq incompressible fluid. They hypothesized that convection (of an electrically insulating fluid) in a deep spherical shell will tend to form quasi-barotropic eddy structures capable of transferring zonal momentum into jets. Although it suggests that deep convection may offer a qualitative explanation for deep jet formation (via quasi-horizontal Reynolds stresses permeating the entire deep interior), it says little about the quantitative scales and amplitudes of motion expected to emerge.
More recent numerical model studies, such as those by Heimpel et al. (2005) based upon Boussinesq deep convection, and anelastic models by Jones and Kuzanyan (2009) and Gastine and Wicht (2012), put the deep convection hypothesis on a more quantitative basis. They showed that the Reynolds stresses induced by deep convection, whether in a Boussinesq context or in the presence of strong fluid compression, can possibly drive jovian-strength zonal flows. However, for pragmatic reasons these models (even when dealing with a realistically compressible fluid) cannot match the key similarity parameters (e.g., Rayleigh and Ekman numbers) for the deep interiors of gas giants by many orders of magnitude. Thus, even though such model simulations appear to be able to produce strong zonal flows that penetrate the entire molecular envelopes of these planets, the heat fluxes necessary to drive these strong flows are much larger than observed.
The recent study by Showman et al. (2011) addresses this issue head on and concludes that, if heat fluxes and values for Rayleigh and Ekman numbers closer to the real values in Jupiter and Saturn are applied, the resulting zonal jets would be at least an order of magnitude too weak to account for the jets observed at the cloud tops. Some deep convection models that exhibit strong zonal jets tend to ignore the braking effects of hydromagnetic forces that may well damp and inhibit deep jet formation (e.g., Liu et al., 2008). Several models in which these forces are accounted for, i.e., dynamo models, produce crude prograde equatorial jets but their high latitude jets are unrealistic (e.g., Heimpel and Gómez Pérez, 2011;Duarte et al., 2013).
The most recent attempts to determine the depth of penetration of cloud-level winds on Uranus and Neptune from observations seem to indicate that they do not extend far into the deep interior -at least at anywhere near the strength observed at the cloud tops (Kaspi et al., 2013). For Jupiter, this issue is open. The Juno mission (Matousek, 2007) launched by NASA in 2011 will hopefully deliver answers.

A scale-dependent energy flux
We now introduce a scale-dependent energy transfer function that approaches W on large scales. Following Sukoriansky et al. (2012) we define a dynamic wavenumber cutoff, n c , and compute the energy flux, W < , from all modes with n < n c to the mean flow. is calculated from its spherical harmonics with the modes n > n c discarded. The flux from the small scale modes with n > n c to the mean flow is then W > = W À W < , with W < (0) = 0 and W > (0) = W. Fig. 11 shows W > as a function of n c /n b for days 1, 2, and 3 for our and CS velocity fields. The positiveness of W > for all n c clearly indicates that the energy flux is directed from all scales to the zonal flow. The flux strength increases with increasing length scale, and on the largest scales exceeds by a factor of 5-20, depending on the dataset.
These results can be compared with the barotropic simulations by Sukoriansky et al. (2012) under the assumption that the barotropic model is applicable to 3D flows where the barotropic mode contains a significant fraction, if not most, of the energy. In these simulations R b was about 2.6, while the ratio W > / fluctuated with a standard deviation of about 2. If the temporal variability of this ratio increases with R b , which is a plausible assumption, then it may be much greater for Jupiter, where R b J 5, than in the simulations. The barotropic model attributes the large-scale variability of W > to complicated interactions between zonal jets, eddies and non-dispersive packets of nonlinear waves called zonons (see, e.g., Sukoriansky et al., 2008). Furthermore, these interactions cause W > to undergo large-amplitude oscillations and even become negative. In other words, the simulations indicate that the energy exchange between the zonal flows and eddies may proceed both ways. Only a long-term averaged value of W is positive with a mean value of W/ around 0.5 for relatively low values of R b . This ratio may approach 1 with increasing R b . It is precisely the averaged value of W that provides a true characterization of the mean energy flux from all non-zonal structures to the zonal flows.
The predictions of the barotropic model could be far reaching. Unfortunately, the variability of the large-scale fluxes cannot be fully assessed using only three days worth of data, but Fig. 11 demonstrates that W may change by a factor of 2 even in a short span of one day. The large-amplitude oscillations of W, as well as the confinement of the averaged value of the ratio W/ between 0 and 1, are all important predictions that need verification in long-term observations.

The Great Red Spot -zonal flow energy exchange
As described by Sukoriansky et al. (2012), when n c is small W > has a dip associated with direct energy flux from the zonal flow to non-zonal structures. This dip is clearly visible in both our and CS's data (Fig. 11), although it is more pronounced in the latter. There exists a local maximum of W > at n c ' 30 and a local minimum at n c ' 20. The non-zonal structures in this range receive energy from the zonal flow as the latter receives less energy from all scales with n > 20 than from all scales with n > 30. The scales with n [ 20 may therefore receive energy from the scales in the range 20-30 and transfer it back to the zonal flow. This behaviour of W > reflects a complicated energy exchange between large eddies and the zonal flow.
One may ask what role the GRS plays in this exchange. We address this question by comparing W > with f W > , the equivalent flux computed with the GRS removed (see Section 5.4). Fig. 12 displays the difference DW > ¼ W > À f W > , for all three days, using both datasets. Positive DW represents a decrease in energy flux to the zonal flow due to removing the GRS, and negative DW the opposite. On the largest scales, energy flux is from the GRS to the zonal flow while it is the opposite on smaller scales. One discerns high variability of DW > as over a span of a few days it varies by a factor of two and may even change sign; the variability is significantly less in the CS data. Nevertheless, the energy flux averaged over the three days is about the same in both datasets, and may exceed by a factor of 2. Apparently, the interaction between the GRS and the zonal flows alone already accounts for a significant part of the disparity between W and . Other coherent eddies may further enhance this disparity.
Another aspect of the two-way energy exchange between the GRS and the zonal flow taking place over a wide range of spatial and temporal scales concerns the possibility that the GRS is a Rossby solitary wave embedded in a strong zonal shear flow (e.g., Maxworthy and Redekopp, 1976a,b;Read, 1987;Nezlin and Sutyrin, 1994;Williams, 1996). The same possibility exists for Jupiter's other coherent vortices as most of them reside in areas of strong meridional shear (Legarreta and Sánchez-Lavega, 2008). If this indeed is the case, then the multi-scale zonal flow -eddy interaction is a product of an ongoing process of focusing the GRS and other coherent vortices by the zonal shear.

Energy conversion conundrums
The interaction between vortices and zonal flows offers a possible resolution to the eddy-zonal energy conversion conundrum mentioned earlier in the case of shallow water dynamics. If large-scale oscillations of W are indeed sustained by the inverse energy cascade from smaller scales, as in the barotropic case, then the cascade rate would be the proper variable for assessing the energy conversion rate per unit mass. Barotropic models suggest that this can be much smaller than instantaneous values of W. A long-term, relatively high resolution monitoring of Jupiter's fully 2D velocity fields would therefore greatly help to finally sort out this puzzle.
The energetics pertinent to the zonostrophic regime suggests an explanation for another energy transfer conundrum as formulated by Ingersoll et al. (2004) for Jupiter and Suomi et al. (1991) for Neptune. The absorbed sunlight (power per unit area) on Jupiter and Neptune is, respectively, only 3.3% and 1/900th of that on Earth, but Jupiter's and Neptune's winds are, respectively, 3-4 and 6-7 times stronger, facts which seem to contradict each other. This contradiction is further accentuated by the arguments advanced in Vallis (2006, p. 144): ''Arguably, the magnitude of the velocity in the atmosphere and ocean is ultimately given by the strength of the forcing, and so ultimately by the differential heating between pole and equator.'' It remains unmitigated by the subsequent qualifying remark by Vallis, ''although even this argument is not satisfactory, since the forcing mainly determines the energy throughput, not directly the energy itself, and the forcing is itself dependent on the atmosphere's response.'' Similar ideas can be traced in the scaling analysis of the deep convection theory by Showman et al. (2011) where the rate of the convective forcing is explicitly related to the strength of zonal jets.
In sharp contrast to the Earth's atmosphere, the presented data demonstrates that Jupiter's troposphere conforms to the regime of zonostrophic macroturbulence. As elaborated in , this regime is distinguished by strong barotropization and kinetic energy concentration in the zonal jets. Both of these features have been confirmed in recent fully 3D simulations. 2 In fact, these simulations showed that the kinetic energy of zonal flows far exceeds not only the barotropic EKE but also the baroclinic EKE. The zonal spectrum computed from the Cassini data appears to be broadly consistent with Eq. (2a) indicating that it may be independent of the forcing. This would further imply that the strength of Jupiter's winds depends not on the absorbed sunlight, intensity of the internal energy sources, or the temperature difference between the equator and pole (which is generally small on all giant planets), but on the global value of b (=X/R) and the magnitude of large-scale 'friction'. The precise nature of the large-scale 'frictional' processes is still poorly known, and is likely non-universal (e.g., radiative cooling, Scott and Polvani (2008), and Ohmic dissipation, Liu et al. (2008), have been proposed), though they appear to be weak on Jupiter and progressively weaker on the other giant planets (Sukoriansky et al., 2002;Ingersoll et al., 2004), at least at the altitudes of their weather layers. The magnitude of this friction can be roughly judged by the number of zonal jets, which are fewest on Neptune . As a consequence, despite being farthest away from the Sun, Neptune nevertheless possesses the most powerful winds in the Solar System.
Summarizing the importance of the inverse cascade in giant planets' circulations (for which R b is sufficiently large), we reemphasize that however small the cascade rate is, it is crucial for maintaining the large-scale circulation in the zonostrophic regime. But the strength of the zonal winds may then be independent of and set mainly by b and the level of large-scale 'friction'. This conclusion is in agreement with Vallis's remark cited earlier that ''the forcing mainly determines the energy throughput, not directly the energy itself.'' The potentially important role of shallow water quasi-geostrophic dynamics on Neptune is further underscored in the recent study by Kaspi et al. (2013) which asserts the comparative shallowness of its zonal jets.

Meridional diffusion on Jupiter
Along with its unique dynamics, the regime of zonostrophic turbulence is characterized by specific laws of meridional and zonal diffusion. These laws provide an independent diagnostic tool for regime identification from observations and may help to understand and quantify large-scale atmospheric diffusion.
In the zonostrophic regime, the diffusion processes are determined by the interaction between turbulence and Rossby-Haurwitz waves. If the source of the diffusing tracer is relatively small, a tracer cloud's spread is at first ballistic, following the dominant wind velocity. When the cloud expands to the size of turbulence-dominated scales, its diffusion follows the Richardson law with a scale-dependent eddy diffusivity (Sukoriansky et al., 2009), where n characterizes the maximum length scale of the diffusing tracer cloud. It is assumed that all scales smaller than n À1 contribute to K / . At later times, when the tracer cloud spreads to length scales Oðn À1 b Þ which are dominated by waves, the coefficient of meridional eddy diffusivity becomes scale-independent and is given by (Sukoriansky et al., 2009; By substituting the expression for n b from Eq. (3) into Eq. (17) we find Since the value of C Z found from observations differs from the one obtained in numerical simulations (e.g., Huang et al., 2001), the coefficient in Eq. (18) may differ from 0.5. Therefore expressions in this section should be understood only as an order-of-magnitude analysis. For Jupiter's troposphere, substituting b ' 2.5 Â 10 À12m À1 s À1 and our earlier estimate of % 10 À5 m 2 s À3 into (18), we find that K / % 2 Â 10 10 cm 2 s À1 . Eqs. (16) and (17) can be construed in terms of mixing length theory, where the meridional diffusivity at a scale l $ n À1 is given by where v = E 1/2 is a measure of the velocity fluctuation corresponding to the length scale l as obtained by integrating the residual spectrum, Eq. (2b), from n to 1. The calculation yields K / $ 1/3 n À4/3 , in agreement with Eq. (16). The largest scale containing turbulent eddies and contributing to the meridional scalar diffusion is Oðn À1 b Þ (Sukoriansky et al., 2009;. Processes on larger scales are dominated by Rossby-Haurwitz waves, which do not contribute to diffusion, and so in Eq. (19) we should take l $ n À1 b and compute v by integrating the spectrum (2b) from n b to 1. This yields Eq. (17). Lacorata and Espa (2012) confirmed this estimate in laboratory experiments.
A scale limitation in terms of n À1 b reflects the effect of Rossby wave elasticity imposed by the wave's restoring force introduced by the meridional variability of the Coriolis parameter, i.e., the b-effect (Baldwin et al., 2007). This force limits the meridional deviation of fluid particles.  explored an analogy between Rossby wave elasticity and internal gravity wave elasticity which manifests via the Ellison-Britter-Osborn (EBO) model for density-stratified fluids. In that case, buoyancy is a restoring force that limits the amplitude of the vertical oscillations of fluid particles and leads to the scale-independence of the vertical (diapycnal) eddy diffusivity on large scales dominated by internal waves.
Such scale-independence of the vertical diffusivity in strongly stratified flows has been confirmed in observations. Ledwell et al. (1998) showed that the diapycnal diffusivity associated with the diapycnal mixing in the oceanic pycnocline remains approximately scale-independent on time scales from 6 to 24 months, despite considerable vertical spread of the tracer billow during that time. In that case, the diapycnal diffusion is carried out within turbulent patches dominated by 3D overturning turbulence. The vertical scale of these patches is characterized by the Ozmidov wavenumber, which depends on the Brunt-Väisälä frequency, measuring the strength of the ambient stratification, and the turbulence dissipation rate. The Ozmidov wavenumber is the stable stratification analogue of the transitional wavenumber n b , and is the relevant parameter for the diapycnal eddy diffusivity. The evolution of an initial tracer concentration profile after one year's integration using a one-dimensional diffusion equation with constant, scaleindependent eddy diffusivity estimated from the EBO model was in good agreement with measurements reported by Ledwell et al. (2011, their Fig. 2).
We note that a barotropic model cannot account for the intricacies of the vertical distribution of spreading material. Instead, it can describe the evolution of the horizontal extent e.g., of impact clouds, under the assumption that they represent advection by the mean flow plus horizontal diffusion in a vertically-integrated sense. In this situation, the meridional mean particle separation (see, e.g., LaCasce, 2008) is given by where t is the time elapsed and K / is given by Eq. (18). This can be compared directly with observations.

Diffusion of acetylene and ethane
The scale-independence of meridional diffusivity in Jupiter's stratosphere emerges by fitting Cassini observations of acetylene (C 2 H 2 ) and ethane (C 2 H 6 ) meridional mixing to a simple diffusion model, which was done by Liang et al. (2005). They found that, for the lower stratosphere below 5 mbar, a good fit to the tracer diffusion can be obtained with K / = 2 Â 10 10 cm 2 s À1 . This value is in good agreement with our own earlier estimate. Above 5 mbar the diffusivity needs to be decreased to about 2 Â 10 9 cm 2 s À1 to fit the data. As mentioned earlier, our relatively simple barotropic model does not explain changes in the vertical. Note that Eq. (18) was derived for the weather layer, which is confined to the upper troposphere and lower stratosphere, so comparisons with the upper stratosphere are, perhaps, overstretched. Even for the lower stratosphere the dynamics are not well understood, and the extent of the zonal flow's penetration into the stratosphere is poorly constrained (although see Flasar et al., 2004). Overall, the scaleindependence of the meridional diffusivity and its numerical agreement with Eq. (18) may be taken as a hint that zonostrophic macroturbulence extends well into the stratosphere. Verification of this assumption as well as the more fundamental question about the nature of the stratospheric energy cascade still remain largely open, however.

Debris diffusion after the Shoemaker-Levy 9 comet impact
Further testing of Eq. (18) as a model for the meridional diffusivity in Jupiter's stratosphere can be done using observations of the dispersion of debris and gases after the Shoemaker-Levy 9 (SL9) comet impact (see, e.g., Harrington et al., 2004) in 1994 and the Wesley asteroid impact on 19 July 2009 Sánchez-Lavega et al., 2011;de Pater et al., 2010;Hammel et al., 2010). Friedson et al. (1999) concluded that, in the case of SL9, advection due to the residual circulation was insufficient to describe the temporal dispersion of the impact cloud fully. To do so, they needed to include meridional diffusion processes and invoke a scale-independent meridional diffusivity with magnitude (1-10) Â 10 10 cm 2 s À1 in the region between 10 and 100 mbar. Our estimate of K / , Eq. (18), agrees well with these results and produces no negative eddy diffusivities (as the Friedson et al. (1999) calculations did). Griffith et al. (2004) traced the meridional transport of hydrogen cyanide (HCN) in Jupiter's stratosphere for six years following the SL9 impact around 44°S. The dispersion after the first ten months was consistent with a scale-independent meridional eddy diffusivity of (2-3) Â 10 10 cm 2 s À1 , spreading to between 30°S and 65°S. This figure is again in good agreement with our estimate using Eq. (18). During the following years, HCN remained confined near 44°S, with close-to-Gaussian diffusion, but also diffused into the northern hemisphere stratosphere with diffusivity (2-5) Â 10 11 cm 2 s À1 . Such a diffusion pattern does not indicate the presence of any mixing barriers.

Debris diffusion after the 2009 Wesley asteroid impact
The 2009 impact event provided further detailed information on short-term meridional dispersion. The impact took place at an approximate planetocentric latitude of 305°W 55°S. The size of the impactor was estimated at about 1 km Hammel et al., 2010). The planet was subsequently observed from July 20 through December 31, until the cloud became indistinguishable from the background. Visible-wavelength images from the Hubble Space Telescope taken four days after the impact reported the meridional and zonal extents of the impact cloud to be 2800 km and 6000 km, respectively. Multi-wavelength studies of the debris and aerosol spread after the impact revealed the cloud's vertical structure de Pater et al., 2010). As mentioned earlier, our analysis is restricted to spatial boundaries of the impact cloud in a vertically-integrated sense. The latitudinal drift of some features observed during the first two weeks following the impact were found to be 2 m s À1 equatorward at 55°S (the middle of the cloud) and 0.8 m s À1 poleward at 60°S . The cloud was confined between latitudes of 53.5°± 0.5°S and 61.4 ± 0.5°S afterwards.
These observations suggest that the final meridional spread of the cloud, L, was about 8°/9000 km after five months. The transitional wave number n b was estimated in Section 5.3 to be around 60, which corresponds to a length scale of L b = pR/n b % 3500 km.
The initial two-week spread can be reasonably estimated at about 3° Sánchez-Lavega et al., 2011). This is comparable with L b , so we can estimate the subsequent meridional mean particle separation using Eq. (20), replacing D with L. Using the value of K / calculated earlier (Eq. (18)) and t = 3 Â 10 7 s (five months), we find D % 5000 km. Added to the impact cloud spread, the total debris dispersion over the five month period is about 8500 km, close to the observed spread.
After the five months had elapsed, the impact cloud seemed to become trapped between two eastward jets . Hence the question arises as to whether or not the eastward jets are impenetrable mixing barriers. We recall that long-time observations by Griffith et al. (2004) of HCN following the SL9 impact revealed no mixing barriers. This result suggests that a fivemonths record is insufficient to assess the barriers' existence.
An interesting final observation is that even though the impactors deposited significant amounts of energy in Jupiter's atmosphere in both the SL9 (de Pater et al., 1995) and the Wesley impacts (Pond et al., 2012), the meridional dispersion seems to have been determined by the background energy of the residual circulation rather than the energy produced by the impacts.
Summarising, the meridional diffusion of material derived from hydrocarbon concentrations and from both the SL9 and Wesley impacts is consistent with a scale-independent diffusivity (Eq. (18)) inherent to a zonostrophic macroturbulence framework. The value of K / agrees well with various data collected over many years. This is intriguing evidence that the zonostrophic regime on Jupiter may possibly be extending into the lower stratosphere. Such a possibility should be investigated in future observations.

Concluding remarks
This paper collates extensive evidence in support of the zonostrophic turbulence regime in Jupiter's large-scale atmospheric circulation. This major result is obtained using two independent methods for extracting fluid velocities from cloud images. Both were based on tracking patches of cloud brightness between images spaced in time to estimate the underlying wind speeds and directions. Both methods were limited in different ways, so the general agreement between the wind fields and the spectral diagnostics obtained using the two methods is very encouraging.
By combining earlier diagnostics with the new ones established in this study, there is clearly a case to be made not only for the existence of inverse energy transfer to large scales, but also for the presence of the KK residual spectrum and the regime of zonostrophic turbulence in general. The relevant diagnostic results are: (1) the flow field is strongly anisotropic in both spectral and physical space, (2) the spectra conform to the spectral laws given by Eqs. (2a) and (2b), (3) the zonostrophy index, around R b % 5À6, is significantly larger than the zonostrophic regime threshold around 2.5, (4) a strong disparity exists between the magnitudes of zonal mean zonal velocities and departures from the zonal mean, (5) the total energy flux W is from nonzonal structures to zonal flows, (6) the scale-wise energy flux is from the eddies to the zonal flow at almost all length scales down to the resolution limit, (7) there is a hint at large temporal fluctuations in W, and (8) the diffusivity given by Eq. (18) not only leads to scale-independence of the meridional eddy diffusivity, but also predicts a value in good agreement with multiple observations. One of the difficulties with the data is that the KK spectral segment with the À5/3 slope only appears for n > 180 or so up to the grid spacing limit of both databases. So, even with the backing of a carefully considered error budget, it is hard to make completely conclusive claims about the existence and interpretation of this part of the spectrum. This also means analyzing images from the other filters is of marginal value, since they are at only half the spatial resolution.
All these results suggest that, although important, the zonal velocity profiles on Jupiter and, probably, other giant planets provide limited information on atmospheric dynamics. Being firstorder statistical moments of a turbulent flow field, they only partially characterize the slow mode of circulation. A more complete understanding of the circulation requires observational information of a high enough quality so that higher-order moments can be retrieved, starting with the second-order moments, the 1D spectra. These spectra illuminate the anisotropic nature of Jupiter's circulation and kinetic energy distribution over different scales. Nevertheless, even the second-order moments are insufficient to assess the energy exchanges between different scales and structures as well as the circulation variability on a variety of spatial and temporal scales. These processes begin to manifest in the third-order moments, i.e., the energy fluxes. Future observational programs should run for periods of weeks, months, and possibly years in order to enable the observers to measure the higher-order moments and their variability. In addition to measuring the dynamic characteristics, these programs should ideally collect simultaneous dynamic and thermodynamic information such that the velocity and temperature covariances and higher-order moments could also be computed.
How important is the inverse energy cascade for planetary dynamics? Some recent studies (e.g., Ioannou, 2003, 2007;Farrell and Ioannou, 2009;Marston, 2012;Srinivasan and Young, 2012) suggest that energy transfer to the zonal flow need not necessarily involve an inverse cascade. This research is only at its early stages and it has not yet been used to produce realistic spectra, analyze large fluctuations of W such as those discussed in Section 6, verify the emergence of zonons, or derive a model for the meridional diffusion. The concept of the inverse energy cascade, on the other hand, is one of the corner stones of the framework of zonostrophic turbulence which is supported by the evidence of the KK range in the residual spectrum. It is well suited to deal with the above phenomena.
The large-amplitude fluctuation of W discussed in Section 6 is a result whose verification is important for further progress in our understanding of Jupiter's atmospheric dynamics. Some light could be shed on this using data already available and waiting to be processed. The task may be somewhat easier than the one undertaken in this paper because only the large-scale behaviour of W needs to be investigated. But such data is only available in high quality for a period of months. Continuous monitoring of Jupiter's atmosphere is necessary for more precise estimates of energy spectra, energy transfers, and other spectral characteristics such as higher-order velocity correlations. The zonostrophic turbulence framework could possibly be extended into the stratosphere, as the zonostrophic meridional diffusivity (Eq. (18)) seems to be applicable there. This result suggests that tropospheric jets may extend further into the stratosphere than previously thought (not inconsistent with observations by Flasar et al. (2004), at least away from the equator). This suggestion should be tested by future observations.
Our results point to a difficulty with present numerical simulations of Jupiter's atmospheric circulation. As shown in Section 5.3, the range of isotropic turbulence starts at wavenumbers n J n z ' 500. This range includes the energetically important isotropic inverse energy cascade that determines the value of . To properly resolve this range, a numerical model would need to possess some 1000 active meridional modes and the correspondingly large number of the zonal modes, not mentioning the vertical resolution. In addition, long-time integrations would be required, due to the slowness of large-scale circulation processes. At present such demands are impractical based on existing computer resources, particularly if fully 3D simulations are required. Typical resolutions for state-of-the-art numerical simulations of Jupiter's atmosphere are around 0.7° (Lian and Showman, 2010).
We need to emphasize that one of the important features setting apart atmospheric circulations on terrestrial planets and gas giants is the direction of energy transfer. While it is direct in the Earth atmosphere (as a result of the large Rossby radius where the baroclinic instability takes place), it is inverse on giant planets (due to small-scale convective forcing).
One may hypothesize that the regime of zonostrophic macroturbulence is established on any rapidly-rotating planet with a system of slowly evolving zonal jets and a sufficiently large value of the zonostrophy index, R b . The number of jets, which can be estimated from visual observations, would give an approximate value for the frictional wavenumber, n fr , and hence the Rhines' wavenumber n R (see, e.g., Sukoriansky et al., 2007). Along with the other two observable parameters, X and R, this yields the total zonal kinetic energy (see, e.g., Galperin et al., 2001;Sukoriansky et al., 2002), which is close to the total kinetic energy (per unit mass) associated with the tropospheric circulation. If the inverse energy cascade rate can be approximated by 1% of the thermal energy received from the planet's star(s), then one could reproduce the entire kinetic energy spectrum of the planet from these observations alone! The meridional diffusivity could be estimated as well. Conversely, if the diffusion characteristics are observable, measurable and indicative of a scale-independent value of K / , they can be used to estimate the rate of the inverse cascade . In summary, the regime of zonostrophic macroturbulence provides a new conceptual framework for understanding and quantification of large-scale circulation on the giant planets in our Solar System.

A.1. Sources of random error in the velocity vectors
There are three main sources of random error in the individual CIV vectors, and further error is introduced in the mosaicing step.
First, in the individual wind vectors produced by CIV, the navigation error in NASA's limb-fitting algorithm is 0.2 pixel in the position of each 0.05°resolution image (Salyk et al., 2006, p. 431), which corresponds to a velocity error of 4.6 cos / m s À1 in the u velocity for a 63 min image separation, and an error of 4.6 m s À1 for the v velocity. The factor of cos / comes from the conversion from pixel displacement to velocity. Second, the position of the feature tracked in the correlation box introduces a 'shear error' (Choi et al., 2007) of up to half the correlation box size multiplied by the local horizontal velocity shear. The horizontal shear is approximately 10 À5 s À1 (Salyk et al., 2006, Fig. 5), leading to a velocity error of about 5.8 cos / m s À1 in u and 5.8 m s À1 in v. Finally, there is a small but non-negligible error of about 0.05 pixel associated with the formal tracking error in the CIV method itself (Fincham and Delerce, 2000), corresponding to an additional 0.8 cos / m s À1 in u and 0.8 m s À1 in v. Combining these errors in quadrature gives a total random error in the u velocity component of about 7.4 cos / m s À1 , and 7.4 m s À1 in v, or 7:4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ cos 2 / p m s À1 in the wind speed. This value is comparable with other groups' analyses (e.g., Choi et al., 2007).
The mosaicing step introduces further error because, at a particular grid point, we combine filtered vectors that have a spread of velocity values as well as individual errors associated with them. Furthermore, each of these vectors is weighted as described above. The combined error is calculated using a method based on Brandt (1970, pp. 97-99) and Upton and Cook (1996, p. 72). Using this, we estimate the u-velocity component error variance as follows. Say we have k vectors, j = 1 Á Á Á k, being combined to define the mosaiced velocity u at a particular point, with u-velocities u j , random errors r u j , and normalized weights x j , as defined in Section 3.3. Then the error variance in the mosaiced velocity is x j r 2 u j þ ðu j À uÞ 2 h i ðA:1Þ and similarly for v. This expression is independent of the method used to combine the image pairs-the method determines the x j and the set of vectors used, but the error variance expression is the same.
A.2. Sources of systematic error in the CIV procedure The main sources of systematic error in the CIV procedure come from applying the method to cloud images. We must assume that (1) there exists negligible vertical motion, (2) the images are horizontal sections through the fluid, and (3) the motion of clouds actually reflects the underlying winds.
We can minimize errors associated with the first of these assumptions by using image pairs separated by 63 min rather than a full Jupiter rotation period. Originally we used the jovian rotation period as the separation, but obtained very poor results, as the flow distortion over this time prevented the algorithm from identifying features accurately. In the raw images there is sufficient overlap between the northern and southern hemisphere images for velocities to be retrieved from these pairs, which are separated by just seven minutes. But as we expected the patch displacement between these 7-min pairs to be swamped by measurement noise and errors, we did not explore this option. Earlier work suggests that in general, the vertical velocity w is very small. Ingersoll et al. (2004) estimate it at around 0.1 cm s À1 from observations, and Zuchowski et al. (2009) put it at 0.5 cm s À1 using a model. In both cases, w 2 ( u 2 + v 2 , and so the contribution to the kinetic energy from the vertical flow is negligible and within the error bars. Errors due to the other assumptions are unavoidable in the jovian context. In a laboratory setting, a horizontal section through the fluid is guaranteed by an imaging method that uses a horizontal laser sheet to illuminate passive tracers. How well cloud motions reflect the underlying winds is still an open question. An ongoing project led by Kunio Sayanagi (Sayanagi et al., 2010) is addressing this question by using a variety of cloud tracking methods within a general circulation model (GCM). The true wind is then known exactly and so it is possible to compare wind fields obtained by cloud tracking directly against the true winds computed in the GCM.

A.3. Varying the mosaic grid spacing
We combined the filtered vectors from the image pairs into mosaics using four different grid spacings: 0.3°, 0.5°, 0.75°, and 1.00°, each using the standard method described in Section 3.3. This was done to see how the quality of the fields, zonal diagnostics, and spectra varied with grid spacing.
The fields (Fig. 3) have some obvious differences. As the grid point separation increases, the mosaics become smoother, contain fewer small-scale features, and the range of vorticity values decreases. At the 0.5°and 0.3°separations, however, some grid points did not have any filtered vectors within range, so there are a handful of gaps in the mosaiced data at 0.5°and a large number of gaps at 0.3°. The large number of gaps makes the 0.3°grid separation unsatisfactory for our purposes.
There are also some obvious differences in the zonal diagnostics (Fig. 4) for different latitudinal grid separations. The extremes of all the profiles are reduced at the larger grid separations, because the averaging is over a larger sample of vectors. At 0.3°grid separation there is a similar problem as there was in the fields; some latitudes contain little or no velocity vectors. When this happens, there are gaps in the profiles, and when only a very few vectors are available, the mean value has a large standard error associated with it.
In the meridional velocity profiles there is an additional sampling effect favouring grid separations that are a multiple of 0.5°, the sampling distance in the CIV procedure. The image pairs in a given hemisphere are aligned in latitude, so the CIV analysis centres the search box at the same latitude points for each image pair. Hence the output set of CIV vectors are at very similar latitudes for each image pair, with a small deviation for individual vectors based on the calculated meridional velocity. When the latitudinal grid spacing is smaller than 0.5°some latitude bands will contain many velocity vectors, but others will only contain vectors where the meridional velocity is very large. The mean meridional velocity will then be very large (with correspondingly large standard errors) at every other latitude point. A similar effect occurs at 0.75°spacing. When the grid points are aligned with the CIV points every other grid box will contain all the vectors from one latitude band but only around half of the vectors in the next. One latitudinal box will be biased towards southward winds (as those vectors will be positioned slightly south of the latitude required to enter the southward box), and the next will be biased towards northward winds, giving a zig-zag pattern in the meridional velocity profile. Both the zonal mean zonal and meridional velocity profiles are included in the Supplementary Material, for the grid spacings missing from Fig. 4. Fig. A.13 shows how the zonal and residual spectra depend on the mosaic grid spacing. The three grid spacings 0.5°, 0.75°, and 1.00°are presented. The differences at large scales are negligible. Fig. A.13(a) demonstrates that the scaling of the smoothed zonal spectrum with Eq. (2a) is robust and emerges at all grid spacings. Fig. A.13(d) presents the compensated residual spectrum C R (n) / E R (n)n 5/3 . In such a representation, the KK spectrum Eq. (2b) should appear as a horizontal straight line. One can see that our lower grid spacing data (1°and 0.75°) do not attain a clear KK scaling. This scaling emerges for n J 180 at the highest grid Table A.3 Various options used to test the sensitivity of the results to the mosaicing method. Variant 7 is the default used to calculate the fields and spectra shown earlier in the paper, and is described fully in Section 3.3. The x m in variant 6 is the smallest possible x m that ensures no space in the individual images is left unsearched.

Method variant
Averaging area dg j weight?
x m do j weight? dt j weight? Within radius x m of grid point Within radius x m of grid point 1/dg j dk 1/do j 1/dt j 8 Within radius x m of grid point 1/dg j 2.dk 1/do j 1/dt j Fig. A.14. Schematic of the alternative grid used when combining image pairs into a mosaic (Table A.3, variants 1-4). We average over a grid box of side x m = dk (=d/).  Table A.3 used to create the mosaic, using the day 2 data only. The axes are the same as in Fig. 6. Variant 7 is the default used for the other figures in this paper. Only the spectra for n P 30 are plotted-there is no visible difference between the spectra for n < 30. For clarity the key lists the lines as they appear from top to bottom in the figure.
spacing of 0.5°. At small scales, the larger grid spacings (and hence larger averaging kernels in the mosaic step) produce less energy, because the larger kernel averages out more of the small-scale features. The CS data (which has a grid spacing of 0.3°) produces the least energy of all, but presumably because of differences in the image processing. As with the fields themselves, our spectrum for 0.3°is qualitatively different from the other grid spacings and so we omit it from the figure.
Our comparison indicates that the 0.5°grid separation is the optimal grid spacing for this data. This oversamples the true resolution of the mosaics by about a factor of two, but minimizes sampling effects and the number of empty grid points, while retaining as much of the small-scale structure as possible.

A.4. Varying the mosaicing method
We tested the sensitivity of the fields and spectra to the method used to combine the filtered vectors in the image pairs into global mosaics. Several different variations on the method described in Section 3.3 were used, which are summarized in Table A.3. Variants 1-4 used a grid box for averaging (Fig. A.14), and variants 5-8 used the grid point averaging mentioned earlier in Fig. 2. The different variants were done with or without time weighting, using either 1/do j or 1/(do j ) 2 weighting, and with different search radii x m in the grid point averaging case.
There were fewer differences between the mosaiced fields for the different variants than for the different grid spacings, and the robustness of the mosaics to the variant used is very encouraging. The main difference between the variants was the number of resulting points with missing data, which was quite sensitive to the variant used. For a given mosaic grid spacing, the grid box variants generally contained many more missing points than the grid point variants, because the area searched for vectors in each grid box is x 2 m while the area searched for each grid point is px 2 m . Similarly the number of missing points decreased for variants 6 through 8, with variant 8 having no missing data within the ±50°l atitude bands.
One trend that was clear was as x m increases from dk= ffiffiffi 2 p to 2Ádk the mosaics become smoother and less detailed, with the range of vorticity values decreasing as extreme values are averaged out. The number of missing points decreases, however, so there is a tradeoff to be made.
In terms of the mosaics themselves, all the other differences between the variants were smaller than this. The fields were always improved slightly by including time weighting. In those cases some-but not all-of the artifacts in regions where image pairs overlap disappeared. There was also a marginally smaller overall error in the velocity components with time weighting included, particularly visible in the overlap regions.
There was slightly more detail visible for the grid box variants than for the grid point variants, because the averaging area is smaller, but this must be traded off against the increased number of missing points. The u and v velocity errors were also marginally smaller for the grid box variants, which is consistent with the larger averaging area for the grid point variants. The larger the averaging area the more vectors combined into one, and hence the larger the spread in velocity values.
Another point against the grid box variants is that they have a slightly poorer true resolution. Using Rayleigh's criterion to estimate the resolution (26.3% drop in the brightness along the line connecting the two objects to be resolved), the grid point variant's resolution is only 0.97°for a grid spacing of 0.5°, while the resolution of the grid box variants is 1.08°. Because the 1/x averaging function is very steep the grid point method's true resolution is only ever marginally larger than the CIV correlation box size, but the grid box method's resolution is sensitive to the mosaic grid spacing.
Overall, there is little to choose between the variants based on the fields produced. Time weighting should be included, and the large averaging radius in variant 8 should be avoided. Some features of the grid box variants give better results, but there are many more missing data points. This could be circumvented somewhat by using a larger grid spacing, but as we have discussed already the larger grid spacing should be avoided for other reasons. Fig. A.15 shows the corresponding eddy kinetic energy spectra for the eight variants. Again, the robustness of the spectrum to the different variants is encouraging, but there are a few points to note.
First, the inclusion of time weighting makes very little difference to the spectra: lines 1, 3, and 5 are virtually the same as lines 2, 4, and 7 respectively. Second, variant 8 (grid point with x m = 2Ádk) is the anomaly, with the eddy energy in that case falling off much faster with n than the others. This makes sense because the averaging radius (and hence the smoothing length) is twice the grid spacing in this case, so small-scale variability is smoothed out much more by this method than by any of the others. Combined with the visible differences in the mosaics for variant 8, a smoothing radius so large should be avoided.
Finally, the grid box variants look like they may be levelling off at the smallest scales, while the grid point variants do not, with slightly more energy in the smallest scales for the grid-box variants. As discussed above, the grid box variants average over a smaller area so may preserve some of the energy on the smallest scales that is smoothed out by the grid point averaging. The levelling off can be linked to the missing data that was discussed above. When there is missing data, values from other days are copied in (in this case from day 3). The copied in values will be spatially uncorrelated with respect to the surrounding values, which will introduce white noise on the grid scale and hence both increase the amount of energy at the smallest scales (from the noise) and cause it to level off (because it is white noise). Hence we should avoid variants that include a large number of missing data points (blocks of missing data copied in should have a smaller effect as these velocities will be correlated on the grid scale). The grid point variants definitely produce fewer missing points, with the larger x m producing fewer missing points.
So on this basis and from the discussion of the different mosaics, variant 7, as presented originally in Section 3.3, seems justified. Residual spectra for the tests introducing uncertainty into the image pointing, using our day 2 data only. The axes are the same as in Fig. 6. The three shades of grey correspond to the different uncertainties. The shaded areas correspond to the range of spectral coefficients calculated over the three instances run in each case.
This uses the grid point method, including time weighting, but with a search radius that is not too large to damp out small-scale features.

A.5. Uncertainty in the image pointing
Uncertainty in the relative position of the image pairs may introduce error into the spectra such that spectral features on length scales larger than individual image footprints could be artifacts. 3 Each jovian day's velocity field is constructed from 10 to 12 image pairs, so one might expect this uncertainty to result in artifacts on scales with total wave number <10.
We tried to determine whether this uncertainty is important by running a test that introduced uncertainty into the positions of the image pairs before combining them into global mosaics. The navigation error according to Salyk et al. (2006, p. 431) is 0.01°in the CB2 filter images. We introduced errors of d = 0.1°, 0.5°, and 1.0°in the positions of the images before they were combined, i.e., much larger than the quoted error. A random number uniformly distributed between ±d was added to the latitude and longitude positions of each image before they were combined. Each case was repeated three times to measure the spread in the spectrum introduced by this uncertainty. This was done for day 2 only. Fig. A.16 shows the residual spectra for the three uncertainties, and one of the vorticity fields from the 1.0°case is included in the Supplementary Material (Fig. S10). From the visible appearance of the vorticity field it is clear that a 1.0°uncertainty is too largeseveral jets do not line up, which would be visible in the original fields (Fig. 3) if the uncertainty were this large. This indicates that any pointing uncertainty is smaller than 1.0°.
On this basis, the spectra in Fig. A.16 give us confidence that any pointing uncertainty that does exist has a small effect on the spectrum. We can ignore the red areas corresponding to 1.0°uncertainty, and the black and blue areas (0.1°and 0.5°uncertainty) have a small spread almost everywhere (except for n % 3 in the 0.5°case). The spread from this uncertainty is either smaller than or comparable to the random error discussed in Fig. 7. Even if the pointing uncertainty is 50 times the quoted pointing error in the position of individual image pairs (i.e., the 0.5°case), the effect this has on the final spectrum is small.

A.6. Wavelength-dependent errors
It is probable that the random error varies with length scale. The spatial resolution of the velocity measurements was discussed in Section 3.4, and is around 1 Â 1°. So on scales close to 1-2°we may expect some degree of inter-correlation between adjacent gridded values. This is slightly complicated by the interpolations that take filtered velocity vectors from a slightly irregular grid, where each point is located at the mean location of the initial and final displacement of each image patch, onto a genuinely regular grid.
On the largest scales, there is scope for generating artifacts due to the mosaicing of data from different image pairs, which we examined in Section A.5 above. However, the kinetic energy and vorticity fields (Fig. 3) do not contain any obvious large-scale artifacts identifiable by eye. There is a possible artifact around the individual image size (around 10-12 longitudinal wavenumbers), but this is not prominent. The human eye is extremely good at detecting weak patterns in noise, so if no obvious large-scale artifacts appear in the maps then it is unlikely they will affect the spectra very much. Between these two limits, i.e., scales between around 2°and 40°, or total wavenumbers between around 10 and 200, there is no reason to anticipate significant intercorrelations between velocities other than those representing physical structures.