Multistage hydraulic stimulation has the potential to greatly expand the production of geothermal in the United States and worldwide. Zonal isolation and limited-entry completion overcome the problem of flow localization and generate hundreds or thousands of conductive fractures throughout a large volume of rock. In contrast, conventional geothermal stimulation designs are bullheaded as a single stage into a vertical or deviated wellbore, resulting in a small number of dominant flow-pathways. In this study, we perform a modeling study to investigate key physical processes and design considerations for a geothermal system created from multistage hydraulic stimulation. We use a simulator that fully integrates a wellbore simulator, a hydraulic fracturing simulator, and a thermal/compositional reservoir simulator. Thermoelastic and poroelastic stress changes are included, which enables the model the simulate mechanical opening (separation of fracture walls) due to cooling during long-term fluid circulation between wells. The simulator can handle the full life-cycle in a single continuous simulation – multistage fracturing (including crack propagation, proppant, limited-entry, etc.) and long-term circulation. We start by reviewing historic background on the application of hydraulic stimulation to improve geothermal energy production. Next, we discuss key uncertainties regarding stimulation mechanism and fracture geometry. Drawing on this background information, we set up simulations of multistage hydraulic fracturing and long-term fluid circulation through an injector/producer pair. The simulations demonstrate how multistage fracturing enables large flow rates and relatively efficient sweep of heat through large volumes of rock. However, the simulations demonstrate how mechanical opening of fractures due to thermal contraction exacerbates thermal short-circuiting. Produced temperature drops rapidly once mechanical opening reaches the production well. Parameters such as well spacing, fracture spacing, and flow rate can be designed to mitigate thermal breakthrough and maximize discounted return on investment. We integrate the simulator with an optimization algorithm to solve a hypothetical engineering design problem to maximize net present value by optimizing well spacing, fracture spacing, and flow rate. The optimization shows how a balance can be struck between rate acceleration and mitigation of thermal breakthrough.