The capability of using imperfect statistical reduced-order models to capture crucial statistics in turbulent geophysical systems is investigated. Much simpler and more tractable block-diagonal models are proposed to approximate the complex and high-dimensional turbulent dynamical equations through parametric closure models. New machine learning strategies are proposed to learn the expensive unresolved processes directly from data. A systematic framework of correcting model errors with empirical information theory is introduced, and optimal model parameters under this unbiased information measure can be achieved in a training phase before the prediction. It is demonstrated that crucial principal statistical quantities in the most important large scales can be captured efficiently with accuracy using the reduced-order model in various dynamical regimes of the flow field with distinct statistical structures.