The temperature of a deep shale reservoir may reach more than 100°C, and the effect of thermal shock on shale hydraulic fracturing has rarely not been considered in previous studies. Based on mesoscopic damage mechanics and the finite element method, a thermo-hydro-mechanical-damage (THMD) coupling model considering temperature, seepage, stress, and damage fields was constructed to investigate the effects of reservoir temperature, convective heat transfer coefficient ( h), in-situ stress difference and bedding plane angle ( αθ) on shale hydraulic fracturing. The results show that multiple hydraulic fractures (HFs) can occur under thermal shock and that HFs control the distribution of seepage, temperature, and stress fields. Reservoir temperature, in-situ stress difference and αθ are primary factors affecting hydraulic fracturing, whereas h is a secondary factor. When the reservoir temperature rises from 50°C to 150°C, the initiation and breakdown pressures decrease by 65.5% and 16.7%, respectively. HFs cross the bedding plane more easily, and fracture complexity is obviously enhanced. A higher h is favourable for slightly reducing the initiation and breakdown pressures, but it has little influence on the fracture complexity. Once the in-situ stress difference is low, there is a high fracture complexity, but HFs are more easily captured by bedding planes to limit the propagation of fracture height. When the in-situ stress difference is high, HFs are more likely to form bi-wing fractures. Whether αθ is too large or small, it is not conducive to improving the fracture complexity. In this study, when αθ is 30°, HFs and bedding planes intersect to form a fracture network. Essentially, thermal shock plays a key role in reducing the initiation pressure and forming multiple HFs during the fracturing process, and fracture propagation mainly depends on the injection pressure. The results can serve as reasonable suggestions for the optimization of shale hydraulic fracturing.