This study describes the geothermal response of the Phlegraean Fields as well as the impact of changes in its thermal and hydrodynamic properties brought on by a deep borehole heat exchanger (DBHE). For this purpose, we have developed a specialized model based on the Galerkin Method (GM) and the iterative Newton–Raphson algorithm to perform a transient simulation of heat transfer with fluid flow in porous media by solving the related system of coupled non-linear differential equations. A two-dimensional domain characterized with an anisotropic saturated porous media and a non-uniform grid is simulated. Extreme characteristics, such as non-uniformity in the distribution of the thermal source, are implemented as well as the fluid flow boundary conditions. While simulating the undisturbed geothermal reservoir and reaching the steady temperature, stream function, and velocity components, a DBHE is placed into the domain to evaluate its impact on the thermal and fluid flow fields. This research aims to identify and investigate the variables involved in the Phlegraean Fields and provide a numerical approach to accurately simulate the thermodynamic and hydrodynamic effects induced in a reservoir by a DBHE. The results show a maximum temperature change of 107.3°C in 200 years of service in the study area and a 65-year time limit is set for sustainable geothermal energy production.