Commercial development of enhanced geothermal systems in low-permeability rocks relies on fracturing treatments to create complex-fracture networks and an appropriate circulation strategy to maintain high flow rates at sufficiently high temperatures. However, it remains challenging to model complex-fracture propagation and heat energy extraction as a whole. This paper develops a fully coupled thermo-hydro-mechanical model to simulate reservoir stimulation and heat production in naturally fractured geothermal reservoirs. The proposed model is validated against a widely used model, TOUGH2, concerning heat sweep in a vertical fracture. This model is then applied to study multi-staged fracturing and geothermal extraction related to a doublet of horizontal wells. The hydro-geomechanical properties are chosen from the Soultz geothermal reservoir at a depth of approximately 3600 m. Numerical results demonstrate that: (1) mixed tensile and shear fracturing can constitute an important stimulation mechanism for naturally fractured geothermal reservoirs; (2) well interlinked, zigzag artificial fractures between injection and production wells readily lead to channeling flow; (3) keeping a segment of horizontal wells open and placing them further apart are beneficial to the formation of sufficiently diffuse flow pathways; (4) increasing well spacing tends to improve thermal performance; however, for the case of a one stage opening, the improvement of heat sweep efficiency is not significant; and (5) an alternating circulation scheme could achieve superior thermal performance. This study establishes an effective modeling workflow for the design and optimization of naturally fractured geothermal reservoirs, and provides an integrated modeling framework for evaluating recoverable energy potential from geothermal reservoirs.