Simulating Thermo-Hydro-Mechanical (THM) problems in variably saturated porous media presents significant challenges due to the high nonlinearity of the coupled processes, and the numerical instabilities that may lead to nonphysical oscillations in pressure, stress, and temperature solutions. In this work, a robust Mixed Finite Element (MFE) scheme based on the lowest order Raviart-Thomas space is developed for the THM problem in deformable unsaturated porous media. The MFE method for fluid flow and heat transport is integrated with the Crouzeix-Raviart (CR) finite element method for the displacement field. To prevent nonphysical oscillations induced by the hyperbolic convection term in the heat transport equation, the MFE method incorporates an upwind scheme. The Method of Lines (MOL) is employed to transform the partial differential equations into a system of nonlinear ordinary differential equations integrated in time with high-order methods using the DASPK time solver. The coupled THM equations are solved simultaneously using a monolithic scheme to avoid splitting errors. The primary unknowns of the developed formulation are the hydraulic head, the temperature, and the displacement vectors that are assigned at the mesh edges. The proposed MFE models are implemented using Fortran code and developed for both poroelasticity and thermo-poroelasticity problems. The developed models are validated by comparisons against analytical and standard finite element solutions for classical benchmarks. Several numerical experiments are performed under saturated and unsaturated conditions to show the efficiency and accuracy of the developed model. The robustness of the proposed MFE model is demonstrated by solutions free of oscillations. Furthermore, a two-phase THM model is established with COMSOL Multiphysics to investigate the CO2 plume evolution and the reservoir stability during the CO2 injection into the sealed Unayzah formation in Saudi Arabia and highlight the application of the THM model in CO2 sequestration field.