On large scales, both oceanic and atmospheric dynamics are well described by geostrophic models. Typically, the Quasi-Geostrophic Equations (QGE) are used to study mid-latitude planetary fluid dynamics by solving the system on a local tangent plane to the sphere using a linear approximation of the Coriolis parameter. However, these models cannot resolve the equatorial region and its consequential effect on mid-latitude dynamics. To address this gap, a fully spherical version of the QGE was recently reintroduced into the scientific literature. Notably, these equations exhibit a Lie-Poisson structure, indicating that the flow is generated by a Hamiltonian, and that an infinite family of constants of motion exists in terms of the prognostic variable, Potential Vorticity (PV). Recently, a numerical scheme was implemented that preserves this Lie-Poisson structure based on an earlier scheme for the 2D rotating Euler equations. This method enables the investigation of geostrophic turbulence on a rotating sphere in the absence of any particular external forcing, as no artificial dissipation is introduced by the numerical method. In this work, we leverage this innovative numerical approach to explore key phenomena in global planetary dynamics such as the critical latitude for zonal jet formation and the anisotropic Rhines barrier.