  We analyse the nonlinear dynamics of the large scale flow in
Rayleigh-B\'enard 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.
