The aim of the present study is to help understand hydraulic fracture propagation in shale formations through numerical simulations. The hydraulic fracture propagation regime in shale is analyzed considering the anisotropic nature of the shale rock formation and slick-water fracturing fluid. It is determined that the dominant mechanism of hydraulic fracture propagation is the so-called transitional regime that is characterized by a negligibly small fluid lag region and zero fluid front pressure. For the modeling of the hydraulic fracture evolution over time, we assume orthotropic linear-elastic rock media and that the flow of the fracturing fluid is governed by the Reynold's lubrication equation. For the discretization of the coupled solid-fluid equations within the 2D plane-strain context we use the extended finite element method for the rock media and the finite volume method for the lubrication equation. The problem of the hydraulic fracture evolution over time is modeled as stable quasi-static crack growth where time is the result of upholding the mass conservation principle between the fluid inflow and the crack volume. The Picard iterative approach is used to solve the discrete non-linearly coupled solid-fluid equations. Our model is verified against several analytical solutions. Subsequently, a five-stage hydraulic fracturing problem is simulated to study the interactions between the different fractures. Results show that the on-going hydraulic fractures are attracted by the pre-existing hydraulic fractures as a result of the change of the local stress field relative to the initial in-situ stress field. For the cases considered, fracture deflections are found to be most extensive with decreasing fracture spacing and in-situ stress difference, but insensitive with the increasing the ratio of Young's moduli.