The paper presents a fully rigorous unsaturated thermo-hydro-mechanical (THM) coupling model based on non-equilibrium thermodynamics and continuum mechanics to investigate the instability mechanisms of boreholes drilled in deep gas-bearing formations. The model incorporates gas-water two-phase seepage, geomechanical loading, and temperature-difference disturbances during deep drilling. Compared to a previous saturated THM model, the new model shows maximum relative errors of 30.2%, 34.6%, 3.4%, and 0.5% for pore pressure, effective radial stress, effective tangential stress, and temperature, respectively. The unsaturated thermal coupling effect significantly influences the effective tangential stress and enhances the time dependence of the failure zone. The cooling condition is favorable for avoiding shear collapse failure, while heating is not. The study also highlights that formations with lower initial water saturation, higher rock permeability, or higher water phase relative permeability are at higher risk of wellbore instability. The thermal conductivity and thermal expansion coefficient of the solid phase have a greater impact on wellbore stability than those of the fluid phases, and the unsaturated THM coupling effect is more significant at lower rock permeabilities. The findings provide valuable insights for wellbore stability analysis during drilling in deep gas-bearing formations.The paper presents a fully rigorous unsaturated thermo-hydro-mechanical (THM) coupling model based on non-equilibrium thermodynamics and continuum mechanics to investigate the instability mechanisms of boreholes drilled in deep gas-bearing formations. The model incorporates gas-water two-phase seepage, geomechanical loading, and temperature-difference disturbances during deep drilling. Compared to a previous saturated THM model, the new model shows maximum relative errors of 30.2%, 34.6%, 3.4%, and 0.5% for pore pressure, effective radial stress, effective tangential stress, and temperature, respectively. The unsaturated thermal coupling effect significantly influences the effective tangential stress and enhances the time dependence of the failure zone. The cooling condition is favorable for avoiding shear collapse failure, while heating is not. The study also highlights that formations with lower initial water saturation, higher rock permeability, or higher water phase relative permeability are at higher risk of wellbore instability. The thermal conductivity and thermal expansion coefficient of the solid phase have a greater impact on wellbore stability than those of the fluid phases, and the unsaturated THM coupling effect is more significant at lower rock permeabilities. The findings provide valuable insights for wellbore stability analysis during drilling in deep gas-bearing formations.