The most effective method for stimulating shale gas reservoirs is horizontal drilling with successful multi-stage hydraulic fracture treatments. Recent fracture diagnostic technologies have shown that complex fracture networks are commonly created in the field. The interaction between preexisting natural fractures and the propagating hydraulic fracture is a critical factor affecting the complex fracture network. However, many existing numerical models simulate only planar hydraulic fractures without considering the pre-existing fractures in the formation. The shale formations already contain a large number of natural fractures, so an accurate fracture propagation model needs to be developed to optimize the fracturing process. In this paper, we first understood the interaction between hydraulic and natural fractures. We then developed a new, coupled numerical model that integrates dynamic fracture propagation, reservoir flow simulation, and the interactions between hydraulic fractures and pre-existing natural fractures. By using the developed model, we conducted parametric studies to quantify the effects of rock toughness, stress anisotropy, and natural fracture spacing on the geometry and conductivities of the hydraulic fracture network. Lastly, we introduced new parmeters Fracture Network Index (FNI) and Width Anistropy (Wani) which may describe the characteristics of the fracture network in shale gas reservoirs. This new knowledge helps one understand and optimize the stimulation of shale gas reservoirs.