Real-time management of deep geothermal reservoirs via Digital Twins demands balancing physical fidelity with computational cost. Specifically, we model coupled Darcy flow and heat transport in heterogeneous media that involve solving complex advection-diffusion equations, for which Neural Operators have emerged as a powerful acceleration tool. Although Fourier Neural Operators (FNO) offer mesh-independent acceleration, they suffer from critical spectral bias, acting as low-pass filters that fail to resolve well singularities and sharp thermal fronts in regimes dominated by advection [1]. We present a framework to overcome these limitations without the overhead of complex Transformers [2]. We use the Singularity Enrichment strategy [3], combining analytical decomposition with Sobolev-Regularized Residual Learning. Instead of learning tabula rasa, we incorporate the fundamental analytical solution to the homogeneous problem as an input, training the network solely on the heterogeneity-induced residual. Further- more, a Sobolev (H 1 ) loss minimizes pressure gradient errors, ensuring physically conservative Darcy velocities (u ∝ ∇P ). Validation against high-fidelity FEM simulations shows that enforcing gradient consistency eliminates the Gibbs phenomenon near injection wells and prevents solution oversmoothing, enabling robust forecasting of thermal breakthroughs.