Accounting for interaction of the soil’s constituents due to temperature change in the design of geo-thermal infrastructure requires numerical algorithms capable of reproducing the coupled thermo-hydro-mechanical (THM) behaviour of soils. This paper proposes a fully coupled and robust THM formulation for fully saturated soils, developed and implemented into a bespoke finite element code. The flexibility of the proposed formulation allows the effect of some coupling components, which are often ignored in existing formulations, to be examined. It is further demonstrated that the proposed formulation recovers accurately thermally induced excess pore water pressures observed in undrained heating tests.