Quantitative predictions of the physical state of the Earth’s subsurface are routinely based on numerical solutions of complex coupled partial differential equations together with estimates of the uncertainties in the material parameters. The resulting high-dimensional problems are computationally prohibitive even for state-of-the-art solver solutions. In this study, we introduce a hybrid physics-based machine learning technique, the non-intrusive reduced basis method, to construct reliable, scalable, and interpretable surrogate models. Our approach, to combine physical process models with data-driven machine learning techniques, allows us to overcome limitations specific to each individual component, and it enables us to carry out probabilistic analyses, such as global sensitivity studies and uncertainty quantification for real-case non-linearly coupled physical problems. It additionally provides orders of magnitude computational gain, while maintaining an accuracy higher than measurement errors. Although in this study we use a thermo-hydro-mechanical reservoir application to illustrate these features, all the theory described is equally valid and applicable to a wider range of geoscientific applications.