Groundwater Heat Pumps (GWHPs) are increasingly adopted for air conditioning in urban areas, thus reducing CO2 emissions, and this growth needs to be managed to ensure the sustainability of the thermal alteration of aquifers. However, few studies have addressed the propagation of thermal plumes from open-loop geothermal systems from a long-term perspective. We provide a comprehensive sensitivity analysis, performed with numerical finite-element simulations, to assess how the size of the thermally affected zone is driven by hydrodynamic and thermal subsurface properties, the vadose zone and aquifer thickness, and plant setup. In particular, we focus the analysis on the length and width of thermal plumes, and on their time evolution. Numerical simulations are compared with two simplified methods, namely (i) replacing the time-varying thermal load with its yearly average and (ii) analytical formulae for advective heat transport in the aquifer. The former proves acceptable for the assessment of plume length, while the latter can be used to estimate the width of the thermally affected zone. The results highlight the strong influence of groundwater velocity on the plume size and, especially for its long-term evolution, of ground thermal properties and of subsurface geometrical parameters.