We investigate the role of higher order moments in ensemble prediction studies using bred initial perturbations, ensemble averaged direct numerical simulations (dns), and an inhomogeneous statistical closure model for 2D turbulence interacting with Rossby waves, mean fields and topography. Our closure methodology extends the previous work of Epstein [1969, Tellus 21] and Fleming [1971, Mon. Wea. Rev. 99], based on a cumulant discard hypothesis, to include all moments higher than second order. Unlike previous quasi-normal closures the quasi-diagonal direct interaction approximation (qdia) model does not need to be Markovianized to guarantee positivity of energy and enstrophy and is valid for both strongly non-Gaussian and strongly inhomogeneous flows. In addition the qdia does not become unbounded for long integration periods as is the case for the eddy damped quasi-normal method [Fleming, 1971, Mon. Wea. Rev. 99]. We compare the qdia to dns for Northern Hemisphere flow during early November 1979 when a large and persistent high-low blocking dipole formed over the Gulf of Alaska. The formation and decay of atmospheric blocking events has typically been associated with a loss of predictability induced by large scale flow instabilities. Results from the qdia statistical closure and dns are shown to be in close agreement even throughout periods of block formation. We also examine distinct regimes of error growth arising from second and higher order effects.