For large-scale multi-physics subsurface applications, digital twins pose several challenges. One of these challenges is providing near real-time predictions for systems governed by partial differential equations. Conventional multi-physics simulations are not real-time compatible, even when faced against state-of-the-art high-performance computing infrastructures, due to their high computational demands. To compensate for the high computational demand, we employ reduced order models. Modeling of geosystems often involves a multitude of different forward solvers. Therefore, non-intrusive methods offer a high potential since they allow for a straightforward transferability of the approach between different solvers. In this work, we employ the non-intrusive reduced basis method, which relies on a combination of proper orthogonal decomposition and various methods from the field of artificial intelligence. They are of high interest because they can capture the complex system behavior, while at the same time offering interpretability. Interpretability is of key importance for the integration in digital twins because it significantly impacts the reliability of the predictions, which is especially important in the context of applications where risks need to be assessed. In this contribution, we use designated case studies of geothermal applications to illustrate the benefits of reduced order modeling for the development of digital twins for subsurface geosystems. Additional challenges arise from the occurrence of scenarios in the prediction phase that have not been considered in the training phase. We demonstrate how the utilization of hierarchical approaches can address this challenge. Furthermore, we illustrate their integration into a Bayesian uncertainty quantification framework, and we discuss possibilities to extend the aforementioned approaches to allow for a continuous integration of observational data.