Geothermal systems are increasingly deployed in densely populated regions, where their efficiency and environmental impact depend on interactions with the surrounding subsurface and neighboring installations. These interactions are governed by subsurface permeability structures, that is, geological layers with varying capacities to transmit water, propagate pressure, and transfer heat, which determine thermal and pressure interference and long-term heat accumulation. High-fidelity numerical models for coupled flow and heat transport have been developed to analyze these effects based on borehole data and geological models, including city-scale applications [1]. However, their high computational cost limits their applicability in time-critical planning and permitting. Recent advances in surrogate modeling, particularly neural-network-based approaches, enable fast approximations of such simulations and support scalable, near-real-time assessment of geothermal system interactions [2,3]. A key remaining challenge is the consistent representation and propagation of uncertainty arising from sparse and heterogeneous subsurface data. In geothermal modeling, spatially variable permeability fields must be inferred from limited pointwise measurements, introducing epistemic uncertainty that propagates through flow and heat-transport simulations and limits the reliability of predictive quantities of interest. More broadly, this reflects a common challenge in computational simulation and machine learning, where high-dimensional input fields with unobserved regions are often treated as deterministic, highlighting the need for principled uncertainty quantification. Our contribution develops an uncertainty quantification framework for geothermal groundwater simulations under sparse and heterogeneous subsurface data from [1]. The framework explicitly targets epistemic uncertainty in spatially distributed subsurface parameters, in particular the full two-dimensional permeability field, within high-fidelity and surrogate models. This enables reuse of simulation samples across scenarios with varying well locations and extends state-of-the-art subsurface simulation software such as PFLOTRAN and FEFLOW.