In order to characterize the global circulation of the subsurface ocean of Jovian and Saturnian moons, we analyse 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 vigour 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 behaviour conforms with the theory of zonostrophic turbulence, as coined by Sukoriansky et al. (2002), 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 at the ice–ocean boundary, we predict an upper bound for the geostrophic zonal velocity ranging from a few centimetres per second for Enceladus to about one metre 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.