Abstract The low-enthalpy geothermal systems are commonly deployed in sedimentary geological settings that feature significant levels of deposition-induced heterogeneity. In this paper, realistic levels of heterogeneity in the form of varying porosity variance and spatial correlation lengths are considered for a 3D geothermal system. Using 2600 computationally intensive numerical simulations of two doublets placed in a checkboard pattern, the influence of well and doublet spacings on performance metrics of low-enthalpy geothermal systems are investigated. The simulations strongly support that in varyingly heterogeneous systems, the lifetimes of operation are shorter, and depending on isotropicity or anisotropicity of correlated heterogeneity, the lifetimes vary. Most notably the anisotropically correlated heterogeneity can lead to either positive impact (by diverting the cold water plume) or negative impact (by facilitating an early breakthrough of cold water plume) on the lifetime of the operation compared to isotropically correlated heterogeneity. We also calculate the boundary of the region around the wells designated as the “license area” (where the cold water front reaches to or where a threshold temperature drop of 1 °C occurs). By doing so, it is found that the operator can assume larger extents (of up to 50%) for the license areas of the aquifer than the ones conventionally assumed. To minimize the impact of heterogeneity on operation, the best practice was found to place the doublets in the same spacings as of the wells. Moreover, it is found that the well distance can be significantly shorter than what is commonly realised for heterogeneous geothermal aquifers.