An understanding of excavation damage zone (EDZ) is critical in safe construction of underground openings, but conventional continuum-based numerical methods such as finite element method or finite difference method have limitations in simulating crack development under coupled thermal-hydro-mechanical (THM) conditions. In this study, a numerical method of predicting EDZ subjected to coupled THM conditions is developed based on discrete element method (DEM), which can simulate propagation of newly formed cracks by adding fictitious joints. A physical model of a high temperature tunnel filled with pressurized water is first used to verify the accuracy and efficiency of the developed numerical method. Based on the case study of the Baoling Mountain tunnel of the Sichuan-Tibet railway, the effects of coupled THM conditions on characteristics of EDZ are investigated. The results show that the two EDZ indices of volume of plastic zones per length and area of cracks per length increase from M-only to THM scenarios. In addition, the shear displacements of the fault in the roof and floor of the tunnel also increase from M-only to THM scenarios. The increasing EDZ size and fault displacement in THM scenario is attributed to the higher ambient rock stress and water pressure in the fictitious joints and the fault. In high temperature regions, thermal stress increases the maximum principal stresses and induces rock damage. Effect of fault dip angle on characteristics of EDZ is also discussed.