This study presents a coupled thermo-hydro-mechanical-phase field model, grounded in a thermodynamic framework, to elucidate the process of hydraulic fracture evolution in thermo-poroelastic media. The thermo-poroelastic characteristics of the rock exhibit variations contingent upon the phase field value. Consequently, the energy stored within the fluid and rock becomes a function of the phase field values, leading to the formulation of novel phase field evolution equations. For discretizing the thermal-hydrodynamic-mechanical-phase field model, the finite element method is employed. An interleaved scheme is utilized to concurrently solve for displacement, pressure, phase field, and temperature. The corresponding iterative structure is established using the Newton–Raphson method. Several hydraulic fracturing scenarios are investigated to showcase the potential and versatility of the presented model. These scenarios include: a) investigating the impact of different reservoir temperatures on hydraulic fracture propagation, b) analyzing the effects of intersection angle and reservoir temperature on fracture propagation in naturally fractured reservoirs, and c) examining the influence of reservoir temperature on hydraulic fracture propagation in grain-scale polymineral rock.