Due to the increasing availability of large datasets provided by high-fidelity numerical simulations (HFNSs) of physical phenomena, data-driven reduced-order models (ROMs) have become a reliable alternative to those described by partial differential equations, especially when dealing with optimization problems that require multiple evaluations of parameter dependent solutions. These models may be built through different methodologies, among which are the physics-informed neural networks (PINNs) and operator learning frameworks (OLFs). While some recent studies suggest that PINNs have some difficulties in handling convective-dominant physical systems, others show OLFs struggle to accurately predict state variables of interest along the testing interval. In the present work, a ROM for the two-dimensional Rayleigh-Bénard convection flow was built through an incremental proper orthogonal decomposition (IPOD) of a high-fidelity dataset followed by a non-intrusive operator inference (OpInf) approach. Then, the predictive capabilities of the model were assessed for two test cases with different high Rayleigh numbers (Ra=109 and Ra=1011). Very good agreements were achieved between the OpInf-based ROMs and the HFNSs for temperature and velocity fields along training and testing intervals. For instance, time-averaged normalized ROM errors along the testing interval for all the flow variables were below 15% at a given probe location when Ra=109. For each case, while the computational fluid dynamics simulation took about 2 h, the reduced model required less than 3 s to predict the system’s states given its initial condition. This fact illustrates the substantial computational time reduction provided by the ROM. Regarding the case in which the Rayleigh number is higher, the ROM had more difficulties in approximating velocity fields due to the existence of more complex turbulence patterns. Finally, the quality of the data decomposition and of the selected regularization technique, often seen in artificial intelligence algorithms, was investigated. The results suggest the OpInf technique may be successfully applied on the surrogate modeling of turbulent buoyancy-driven flows.