This thesis focuses on the evolution of zonal flows in self-gravitating accretion discs and their resulting effect on disc stability; it also studies the process of disc gravito-turbulence, with particular emphasis given to the way the turbulent state is able to extract energy from the background flow and sustain itself by means of a feedback. Chapters 1 and 2 provide an overview of systems involving accretion discs and a theoretical introduction to the theory of accretion discs, along with potential methods of angular momentum transport to explain the observed accretion rates. To address the issue of the gravito-turbulence self-sustenance, a compressible non-linear spectral code (dubbed CASPER) was developed from scratch in C; its equations and specifications are laid out in Chapter 3. In Chapter 4 an ideal (no viscosities or cooling) linear stability analysis to non-axisymmetric perturbations is carried out when a zonal flow is present in the flow. This yields two instabilities: a Kelvin-Helmholtz instability (active only if the zonal flow wavelength is sufficiently small) and one driven by self-gravity. A stability analysis of the zonal flow itself is carried out in Chapter 5 by means of an axisymmetric linear analysis, using non-ideal conditions. This considers instability due to both density wave modes (which give rise to overstability) and slow modes (which result in thermal or viscous instability) and, thanks a different perturbation wavelength regime, represents an extension to the classical theory of thermal and viscous instabilities. The slow mode instability is found to be aided by high Prandtl numbers and adiabatic index γ values, while quenched by fast cooling. The overstability is likewise stabilised by fast cooling, and occurs in a non-self-gravitational regime only if γ ≲ 1.305. Lastly, Chapter 6 illustrates the results of the non-linear simulations carried out using the CASPER code. Here the system settles into a state of gravito-turbulence, which appears to be linked to a spontaneously-developing zonal flow. Results show that this zonal flow is driven by the slow mode instability discussed in Chapter 5, and that the presence of zonal flows triggers a non-axisymmetric instability, as seen in Chapter 4. The role of the latter is to constrain the zonal flow amplitude, with the resulting zonal flow disruption providing a generation of shearing waves which permits the self-sustenance of the turbulent state.