Enhanced Geothermal Systems (EGS) are those geothermal reservoirs artificially fractured to create paths for injected low-temperature fluid which is then heated up along the flow path until production for electricity generation. This heat recovery involves three tightly coupled processes: thermal, hydraulic and mechanical which interacts with each other and in turn affects the energy production. The local temperature field would be disturbed by injected fluid, resulting in thermal/poroelastic responses near the hydraulic fractured area which are the dominant factors of fluid flow. In this paper, the three-dimensional (3D) Embedded Discrete Fracture Model (EDFM) was adopted to describe the geometry of the fracture and simulate fluid flow and heat transfer between fractures and the matrix, while mechanics, including displacement of the strong discontinuity (fractures), was solved by the 3D eXtended Finite Element Method (XFEM). With the capability of modeling fractures of arbitrary shapes within a 3D reservoir domain using 3D EDFM-XFEM, a coupled THM model was developed based on the unconditionally stable fixed-stress split sequential-implicit method, where the fluid flow/heat transfer module and mechanics module are solved iteratively until convergence within a time step. Fluid flow/heat transfer and XFEM with internal/external tractions are both validated by comparison with existing simulators. We conducted simulations for two synthetic geothermal reservoir heat recovery cases to investigate the effects of the injection temperature and boundary traction condition on the production temperature and fracture deformation. The results indicate that the fracture aperture and permeability are sensitive to temperature variation and hence impact the production rate/temperature. Thermal strain might be the dominant factor of rock deformation, especially in the shallow depth where geostress is at a low level.