AbstractEfficient and accurate modeling of the coupled thermal‐hydraulic‐mechanical‐chemical (THMC) processes in various rock formations is indispensable for designing energy geo‐structures such as underground repositories for high‐level nuclear wastes. This work focuses on developing and verifying an implicit finite element solver for generic coupled THMC problems in geological settings. Starting from the mass, momentum, and energy balance laws, a specialized set of governing equations and a thermoporoelastic constitutive model is derived. This system is then solved by an implicit finite element (FE) scheme. Specifically, the residuals and the Jacobians are scripted in a user‐defined element (UEL) subroutine which is then combined with the general‐purpose FE software Abaqus Standard to solve initial‐boundary value problems. Considering the complexity of the system, the UEL development follows a stepwise manner by first solving the coupled hydraulic‐mechanical (HM) and thermal‐hydraulic‐mechanical (THM) equations before moving on to the full THMC problem. Each implementation step consists of at least one verification test by comparing computed results with closed‐form analytical solutions to ensure that the various coupling effects are correctly realized. To demonstrate the robustness of the algorithm and to validate the UEL, a three‐dimensional case study is performed with reference to the in‐situ heating test of ATLAS at Belgium in 1980s. A hypothetical radionuclide leakage event is then simulated by activating the chemical‐concentration degree of freedom and prescribing a constant high concentration at the heater's surface. The model predicts a limited contaminated regime after six years considering both diffusion and advection effects on species transport.