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.