<p>Defining the mechanism of hydraulic fracture (HF) propagation in geothermal reservoirs, which is influenced by thermal stress and natural fractures (NFs), is crucial for improving reservoir fracturing stimulation performance and improving the effectiveness of geothermal energy development. In this paper, we establish a field-scale 3D thermos-hydro-mechanical coupled fracture propagation simulation method based on the continuum‒discontinuum element method (CDEM). The accuracy of this method is verified through a theoretical model and high-temperature hydraulic fracturing experiments. We conduct numerical research on hydraulic fracture propagation in geothermal reservoirs containing NFs under various control factors. The results indicate that hydraulic fracturing in geothermal reservoirs containing NFs forms a belt-shaped fracture zone perpendicular to the direction of the minimum geostress. Improving thermal stress and increasing the NF density can significantly improve the effectiveness of hydraulic fracturing, potentially increasing the stimulated reservoir volume (SRV) by 17.2% and 63.8%, respectively. As the thermal expansion coefficient increases and the horizontal stress difference, injection displacement, fracturing fluid viscosity, and NF friction coefficient decrease, hydraulic fracturing in geothermal reservoirs becomes more conducive to the formation of complex fracture networks and the improvement of the SRV. These research findings provide valuable theoretical guidance for hydraulic fracturing construction in geothermal reservoirs.</p>