Scientific machine learning methods that incorporate physics knowledge on a data-driven learning have become quite promising for the representation and prediction of nonlinear fluid flow systems with multiple scales in space and time. This work addresses one of these methods, the Operator Inference (OpInf), in the context of model order reduction. By solving a multivariable regression problem in latent space, whose basis is computed through a proper orthogonal decomposition (POD) of the respective high-fidelity dataset, the OpInf seeks for optimal low-dimensional operators that represent the system dynamics. However, this method still requires improvements in its regularization strategy and reliability in complex scenarios, as well as in the robustness of the obtained reduced models for long-term extrapolation when trained with limited data. For that, a recent and efficient algorithm for hyperparameters search, a sequential operator inference and an ensemble learning strategy were successfully implemented in the present work. Other modifications to the standard OpInf were also investigated, such as an incremental data reduction and POD-based forcing terms. To test them, different physical systems were considered: transient heat conduction; oscillating lid-driven cavity flow; nonlinear wave propagation; natural convection; atmospheric CO2 dispersion and sea surface height elevation due to tidal surges. Overall, it was demonstrated OpInf-based models may have very good predictive capabilities for highly turbulent flows and parameter-dependent systems. Furthermore, it was shown these models may be employed for fast response climate-related predictions as they are capable to handle noisy geospatial measurements. Finally, the results suggest the OpInf may be a reliable alternative to deep learning neural networks for model order reduction due to its lower computational costs and good performance beyond the training horizon.