Within the context for deep geological disposal (DGD) of high-level radioactive waste (HLW), thermo-hydro-mechanical (THM) coupled numerical modeling has become significantly important for studying the safe disposal of HLW. In this work, a 3D mechanical module is incorporated into the thermal–hydraulic (TH) coupled code TOUGH2, thus forming an integrated THM coupled simulator referred to as TOUGH2Biot. The Galerkin finite element method is used to discretize the space for rock mechanical calculation. The mechanical process is sequentially coupled with the fluid and heat flow processes, which further gives feedback to the flow through stress-dependent hydraulic properties (e.g., porosity and permeability). Based on the available geological data at the Meuse/Haute-Marne Underground Research Laboratory (MHM URL) in France, the improved simulator is used to analyze the coupled THM behaviors of the Callovo-Oxfordian claystone (COx) induced by thermal loading. The anisotropy of material parameters (e.g., permeability and thermal conductivity) caused by the bedding and of in-situ stresses are well considered in our model. The numerical simulation can reasonably reproduce the field observations, including changes in temperature and pore pressure at monitoring boreholes during the ALC1604 experiment. The modeling results indicate that the anisotropic effects are remarkable, and temperature, pore pressure, and effective stress along the bedding increase more rapidly than in the vertical direction. Insight into numerical results through the visual model is beneficial for helping us to interpret the field observations and to understand the complex THM problem in the COx claystone formation. The numerical method and the modeling results presented in this work can be effectively used in support of performance assessment studies of HLW disposal sites to build confidence in the safety of future applications of nuclear energy systems.