An understanding of the subsurface and its associated physical processes is crucial to address challenges such as providing sustainable energy resources. However, these physical processes are high-dimensional and coupled yielding nonlinear hyperbolic partial differential equations that need to be solved for large spatial and temporal domains. This presents major challenges for intrusive reduced order modeling methods. Furthermore, purely data-driven methods are not desirable either since predictions have to be performed and explainable models are required. To overcome these challenges, we use a non-intrusive reduced basis method. We illustrate the potential and open challenges of the methodology for a broad range of geoscientific and geophysical applications. Using the real-case study of a geothermal application, we demonstrate the great advantages that this method offers to enable extensive parameter estimation studies, which in turn, improve our understanding of geothermal systems. Furthermore, we discuss challenges regarding the sampling of the input parameter space for geodynamical applications, where the computation time of a single forward evaluation can require days on a high-performance computing infrastructure.