In the transport of heat in porous media, diffusion generally dominates over advection due to slow fluid velocities imposed by low permeability. This is the reason why standard Galerkin formulation leading to extra non-symmetric matrix terms may be still used successfully. However, in the presence of fractures the situation may be different. Fractures constitute preferential flow paths where fluid velocities may be significant and advection may become dominant over diffusion (“large advection” with Péclet number >1). This paper focuses on the formulation, numerical implementation and verification of a model to solve the steady-state heat transport problem with large advection along geomechanical discontinuities represented by zero-thickness interface elements. The fluid velocity field is considered as known input data (no hydraulic coupling). The existing SUPG method is modified for its application to zero-thickness interface elements, and the resulting formulation is implemented in an existing FE geomechanical code. An example of application is presented with large advection along a discontinuity crossing a low permeability domain. The results show that the proposed approach leads to stable results, in contrast to standard Galerkin.