Thermal convection in a rotating spherical shell with free-slip boundaries can excite a dominant mean zonal flow in the form of differentially rotating cylinders concentric to the principal rotation axis. This process is studied numerically for Prandtl numbers of order 1, Ekman numbers in the range E = 3 × 10−4−10−5, and Rayleigh numbers up to 100× critical. Small-scale convection transfers kinetic energy into the mean zonal flow via Reynolds stresses. For low Ekman number and high Rayleigh number, the force balance is predominantly among the Coriolis, inertial and buoyancy forces, and viscosity plays a minor role. A modified Rayleigh number Ra* is introduced, which does not depend on viscosity or thermal diffusivity, and asymptotic scaling laws for the dependence of various properties on Ra* in the limit of negligible viscosity (E → 0) are estimated from the numerical results. The ratio of kinetic energy in the zonal flow to that in the non-zonal (convective) flow increases strongly with Ra* at low supercritical Rayleigh number, but drops at high values of Ra*. This is probably caused by the gradual loss of geostrophy of the convective columns and a corresponding decorrelation of Reynolds stresses. Applying the scaling laws to convection in the molecular hydrogen envelopes of the large gas planets predicts the observed magnitude of the zonal winds at their surfaces.