Abstract: To investigate the forming mechanism of complex fracture network during the process of hydraulic fracturing in shale gas reservoirs, numerical simulation of fracture network propagation has been carried out. Taking the e ects of stress shadowing and natural fractures into account, coupled mathematical model of fluid flow in the wellbore and fractures is established. The stress and displacement discontinuity are solved by using displacement discontinuity method. New iterative algorithm of pressure and fracture widths has been constructed and solved with Newton iteration method. The mathematical model and numerical algorithm are validated by comparing numerical solutions to analytical solutions with consideration of leak-o e ect. Parameters sensitivity analysis is performed to study the influencing factors of hydraulic fractures propagation. As to simultaneous propagation of multiple cluster fractures: the smaller fracture spacing is, the more unevenly fracturing fluid is distributed to each fracture; the fracture near the heel of wellbore receives more fluid than other fractures and becomes wider; When taking the e ect of natural fractures into account, the smaller approaching angle is or the weaker stress anisotropy is, the more possibly hydraulic fractures change propagation direction and the more complex fracture networks become.