Efficient and accurate hydrothermal and mechanical mathematical models in porous media constitute a fundamental tool for improving the understanding of the subsurface dynamics in volcanic areas. We propose a finite-difference ghost-point method for the numerical solution of thermo-poroelastic and gravity change equations. The main aim of this work is to study how the thermo-poroelastic solutions vary in a realistic description of a specific volcanic region, focusing on the topography and the heterogeneous structure of Campi Flegrei (CF) caldera (Italy). Our numerical approach provides the opportunity to explore different model configurations that cannot be taken into account using standard analytical models. Since the physics of the investigated hydrothermal system is similar to any saturated reservoir, such as oil fields or CO2 reservoirs produced by sequestration, the model is generally applicable to the monitoring and interpretation of both deformation and gravity changes induced by other geophysical hazards that pose a risk to human activity.