A cohesive zone model (CZM) has been developed to couple fluid flow with elastic, plastic and damage behavior of rock during hydraulic fracturing in naturally fractured formations. In addition to inelastic deformations, this model incorporates rock anisotropics. Fracture mechanics of microcrack and micro-void nucleation and their coalescence are incorporated into the formulation of the CZM models to accurately capture different failure modes of rocks. The performance of the developed elastoplastic and CZM models are compared with the available data of a shale play, and then the models are introduced into a commercial finite element package through user-defined subroutines. A workflow to derive the required model parameters for both intact rock and cemented natural fractures is presented through inverse modeling of field data. The hydraulic fractures' growth in the reservoir scale is then simulated, in which the effect of fluid viscosity, natural fracture characteristics and differential stresses on induced fracture network is studied. The simulation results are compared with the available solutions in the literature. The developed CZM model outperforms the traditional fracture mechanics approaches by removing stress singularities at the fracture tips, and simulation of progressive fractures without any essential need for remeshing. This model would provide a robust tool for modeling hydraulic fracture growth using conventional elements of FEA.