We propose a novel highly efficient, second-order accurate, BDF2-SAV based numerical scheme for a class of nonlinear models that are of importance in geophysical fluid dynamics. The scheme is highly efficient in the sense that only the resolution of a fixed symmetric positive definite linear problem is involved at each time-step. The solutions to the scheme are uniformly bounded for all time provided that the external forcing is uniformly bounded in time. We show that the scheme is able to capture the long-time dynamics of the underlying geophysical model, meaning that the global attractors as well as the invariant measures of the scheme converge to those of the original model under appropriate assumptions. Applications to Lorenz 96 model as well as the incompressible 2D Navier-Stokes system will be discussed.