Economic production from shale gas cannot be achieved by natural mechanisms alone; it requires technologies such as hydraulic fracturing in multiple stages along a horizontal wellbore. Developing numerical models for hydraulic fracturing is essential since a successful fracturing job in a shale formation cannot be generalized to another due to different shale characteristics, and restricted access to the field data acquisition. Empirical methods and Linear Elastic Fracture Mechanics (LEFM)-based numerical techniques are still the prevailing design tools in most of the hydraulic fracture applications though they provide a reasonable prediction only for hard (brittle) rocks. The plastic zone and softening effects at the fracture tip have been neglected in modeling hydraulic fracturing, which results in the prediction of conservative fracture geometry and imprecise fracture pressure. These effects can be identified by the Cohesive Zone Model (CZM) but not by LEFM. Moreover, CZM is able to model fracturing interfaces; for instance, natural fractures, which are mechanically weaker than the adjoining materials. Brittle shales are more likely to contain more natural fractures while ductile shales act as good seals for the fractured layers. In this work, we modeled single and double-stage fracturing in a quasibrittle shale layer using an improved CZM for porous media besides including the material softening effect and a novel boundary condition treatment, using infinite elements surrounding the solution domain of interest. We illustrated some superior aspects of CZM compared to LEFM in shales. We demonstrated the significance of Young's modulus, Poisson's ratio, pumping rate, viscosity, and leak-off in the pumping pressure, and fracture aperture. Also, we found that the stress shadow effects on induced fractures significantly influenced the characteristics of subsequent induced fractures. Moreover, we investigated various scenarios in number and sequence of fracturing stages. Our fully coupled analysis provided us with a continuum-based leak-off model eliminating the requirement for an explicit leak-off model such as Carter's leak-off model.
