We analyse the nonlinear dynamics of the large scale flow in Rayleigh-Benard convection in a two-dimensional, rectangular geometry of aspect ratio $Gamma$. We impose periodic and free-slip boundary conditions in the streamwise and spanwise directions, respectively. As Rayleigh number Ra increases, a large scale zonal flow dominates the dynamics of a moderate Prandtl number fluid. At high Ra, in the turbulent regime, transitions are seen in the probability density function (PDF) of the largest scale mode. For $Gamma = 2$, the PDF first transitions from a Gaussian to a trimodal behaviour, signifying the emergence of reversals of the zonal flow where the flow fluctuates between three distinct turbulent states: two states in which the zonal flow travels in opposite directions and one state with no zonal mean flow. Further increase in Ra leads to a transition from a trimodal to a unimodal PDF which demonstrates the disappearance of the zonal flow reversals. On the other hand, for $Gamma = 1$ the zonal flow reversals are characterised by a bimodal PDF of the largest scale mode, where the flow fluctuates only between two distinct turbulent states with zonal flow travelling in opposite directions.