Building on the previous hydraulic model, a thermally coupled 2D model was created. The upper boundary was fixed at 10 °C to represent atmospheric conditions, while a constant heat flux of 0.065 Wm^-2 was imposed at the lower boundary. Matrix permeability (1×10^-20 m²) and porosity were set to exceedingly low values (1×10^-4) so that heat transfer occurs solely by conduction, with no convective contribution.

To assess the effect of fractures, a second 2D model incorporated three fracture models as 1D elements and distributed loosely across the domain (see figure on the left). The model depth was extended to 3,000 m below sea level. The matrix was assigned values of 1×10^-16 m² for permeability and 1×10^-4 for porosity, whereas the fractures were assigned much higher permeability (1×10^-8 m²) and porosity (1×10^-3). The resulting temperature contour plot shows that the fractures significantly disrupt the flow field: water enters the fractures and diffuses along them. An increased flow velocity profile is observed, particularly in the central fracture (see figure on the right).

Finally, the approach was extended to a preliminary 3D model that includes hydraulic movement and heat transfer but no fractures. The same boundary conditions were applied (10 °C at the top surface and a 0.065 Wm^-2 heat flux at the bottom). Matrix permeability and porosity were again set to 1×10^-16 m² and 1×10^-4. The simulation provides the 3D temperature field, as well as vector fields for heat flux and velocity magnitude. These results indicate that even without fractures, the interaction of fluid flow and conduction generates complex temperature and flow patterns, providing a solid foundation for subsequent fully fractured 3D analyses.