Hydraulic fracturing in unconventional reservoirs exhibits several interesting phenomena including the interaction of hydraulic fractures with multi-scale heterogeneities such as natural fractures, stress/barrier layers, bedding planes, shale laminations, and mineralogy. Moreover, hydraulic fractures originating from different clusters or stages in a multi-stage, multi-cluster treatment interact among themselves. Mathematical models, with various degrees of numerical complexity, are developed for gaining better insights into the physics governing these phenomena. Peridynamics-based hydraulic fracturing model developed by Ouchi (2016) has been demonstrated to capture all of these phenomena. However, its major drawback is that it is computationally expensive. In this dissertation, we have extended the capabilities of the model to multi-phase flow and made it significantly faster by coupling it with the less expensive Finite Volume Method. The single-phase peridynamics flow model for slightly compressible, Newtonian fluids has been generalized for multiphase, multicomponent flow of compressible, non-Newtonian fluids. The generalized flow model has been coupled with the fracturing model and compared with laboratory experiments performed under low confining stresses. The extended model is also applied to simulate the growth of fractures from a new (child) well in the presence of depleted regions created by production from the fractures of an old (parent) well under high confining stresses. The interaction of a hydraulic fracture (HF) with a natural fracture (NF) is investigated. Remote shear failure of the NF due to the pororelastic stress changes caused by the propagating HF are considered. Consistent with the experiments, the remote shear failure is shown to result in the bending of the HF towards the NF before intersecting with it. Accounting for the effects of remote shear failure and poroelasticity, numerical crossing criteria for the HF-NF interaction are developed. The hydraulic fracturing model based on peridynamics (PD) theory is coupled with the less expensive Finite Volume Method (FVM), following the PD-FEM coupling method proposed by Galvanetto et al. (2016). Significant improvements in computational performance are achieved by the coupled model relative to the pure PD-based model, without compromising the unique original capabilities. By monitoring material damage in remote heterogeneous regions, a workflow for estimating the extent of the Stimulated Reservoir Volume (SRV) around a primary hydraulic fracture is developed. A sensitivity study for the effects of elastic properties of the formation, injection rate, and the reservoir fluid type on SRV extent is presentedPetroleum and Geosystems Engineerin
