In this work, a 2D pore pressure cohesive zone model is presented to simulate the hydraulic fracture propagation in naturally fractured formation, in which the fracturing process is governed by a bilinear cohesive zone model equipped with Coulomb's friction law and the fluid flow within the hydraulic fracture is described by the lubrication equation. In this model, fluid leak‐off into rock matrix is ignored and a fully explicit temporal integration scheme is adopted to overcome the convergence issue of conventional implicit scheme (eg, Newton‐Raphson method). The advantage of the proposed model is that it does not require any special crossing criterion to determine the interaction behavior when the hydraulic fracture hits the natural fracture. Implementation of the model was described in detail, and then, the model was verified with well‐known analytical solution of KGD problem and criteria of hydraulic fracture crossing natural fracture. The capability of the proposed model to capture the interaction behavior between hydraulic fracture and natural fracture was demonstrated by the good agreement of the modeling result and analytical solution. Several numerical cases were performed to investigate the impact of key factors on fracture network evolution during hydraulic fracturing treatment. Results show that fracture network is not only governed by the spatial distribution of natural fracture, but also greatly affected by the initial hydraulic aperture of natural fracture and in situ stress contrast.