We introduce a highly efficient, second-order time-marching scheme for nonlinear geophysical fluid models, designed to accurately approximate invariant measures - that is, the stationary statistical properties of the underlying dynamical system. Beyond second-order accuracy in time, the scheme is particularly well suited for long-time simulations due to two key features: (i) it requires solving only a fixed symmetric positive-definite linear system with constant coefficients at each step, and (ii) it guarantees long-time stability, producing uniformly bounded solutions in time for any bounded external forcing, regardless of initial data. We rigorously prove convergence of both global attractors and invariant measures of the discrete system to those of the continuous model in the vanishing time-step limit. Finally we discuss recent progress in the design of ensemble version of the method, and some challenges in the computation of invariant measures.
