Climate-driven adoption of groundwater heat pumps introduces large-scale spatio-temporal thermal interactions that traditionally require expensive physics-based simulations on high-performance computing systems. In this work, we present a machine learning framework that replaces such simulations with a fast, scalable, and physics-aware surrogate model. We propose a Local–Global Convolutional Neural Network (LGCNN) architecture for predicting subsurface temperature fields over domains of arbitrary size and variable numbers of heat sources. The core idea is to decompose the learning problem into local and global components that reflect the underlying physical structure. Local interactions are modeled using convolutional neural networks that learn nonlinear, short-range effects directly from data, while long-range dependencies are handled through an explicit global solver embedded into the pipeline. These components are combined in a sequence of physics-informed stages, enabling stable and interpretable predictions beyond standard end-to-end black-box approaches. Another key contribution is the demonstration of data-efficient training: the model achieves accurate predictions using fewer than five simulated or measured samples, highlighting strong generalization capabilities. The proposed approach is directly applicable to any domain size and number of sources, making it suitable for city-scale optimization tasks. Ongoing work extends the framework toward uncertainty quantification, computational acceleration, real-world benchmarking, and time-dependent modeling, including seasonal dynamics, further positioning LGCNN as a practical surrogate for large-scale physical simulations.