  In order to characterize the global circulation of the subsurface ocean of
Jovian and Saturnian moons, we analyze the properties of 21 three-dimensional
simulations of Boussinesq thermal convection in a rapidly rotating spherical
shell. Flow is driven by an adverse temperature contrast imposed across the
domain, and is subjected to no-slip boundary conditions. We cover a region of
parameter space previously unexplored by global simulations, both in terms of
rapid rotation and vigor of convective forcing, closer to, yet still admittedly
far from, the conditions appropriate for the subsurface ocean of Ganymede,
Europa, Enceladus, and Titan. Our most extreme simulations exhibit a dynamic
global circulation that combines powerful east-west zonal jets, planetary
waves, and vortices. A spectral analysis of the kinetic energy distribution
performed in cylindrical geometry reveals a high degree of anisotropy of the
simulated flows. Specifically, the axisymmetric zonal energy spectra follow a
steep $-5$ slope in wavenumber space, with the energy amplitude exclusively
controlled by the rotation rate. In contrast, the non-axisymmetric residual
spectra display a gentle $-5/3$ slope, with the energy amplitude controlled by
the thermal buoyancy input power. This spectral behavior conforms with the
theory of zonostrophic turbulence and allows us to propose tentative
extrapolations of these findings to the more extreme conditions of icy
satellites. By assuming that kinetic energy dissipates via Ekman friction we
predict an upper bound for the zonal velocity ranging from a few centimeters
per second for Enceladus to about one meter per second for Ganymede, with
residual velocities smaller than the zonal velocity by an order of magnitude on
each moon. These predictions yield typical jets size approaching the ocean
depth of Titan, Ganymede and Europa and $10$ to $40\%$ of the ocean depth on
Enceladus.
