Successful creations of multiple hydraulic fractures in horizontal wells are critical for economic development of unconventional reservoirs. The recent advances in diagnostic techniques suggest that multi-fracturing stimulation in unconventional reservoirs has often caused complex fracture geometry. The most important factors that might be responsible for the fracture complexity are fracture interaction and the intersection of the hydraulic and natural fracture. The complexity of fracture geometry results in significant uncertainty in fracturing treatment designs and production optimization. Modeling complex fracture propagation can provide a vital link between fracture geometry and stimulation treatments and play a significant role in economically developing unconventional reservoirs. In this research, a novel fracture propagation model was developed to simulate complex hydraulic fracture propagation in unconventional reservoirs. The model coupled rock deformation with fluid flow in the fractures and the horizontal wellbore. A Simplified Three Dimensional Displacement Discontinuity Method (S3D DDM) was proposed to describe rock deformation, calculating fracture opening and shearing as well as fracture interaction. This simplified 3D method is much more accurate than faster pseudo-3D methods for describing multiple fracture propagation but requires significantly less computational effort than fully three-dimensional methods. The mechanical interaction can enhance opening or induce closing of certain crack elements or non-planar propagation. Fluid flow in the fracture and the associated pressure drop were based on the lubrication theory. Fluid flow in the horizontal wellbore was treated as an electrical circuit network to compute the partition of flow rate between multiple fractures and maintain pressure compatibility between the horizontal wellbore and multiple fractures. Iteratively and fully coupled procedures were employed to couple rock deformation and fluid flow by the Newton-Raphson method and the Picard iteration method. The numerical model was applied to understand physical mechanisms of complex fracture geometry and offer insights for operators to design fracturing treatments and optimize the production. Modeling results suggested that non-planar fracture geometry could be generated by an initial fracture with an angle deviating from the direction of the maximum horizontal stress, or by multiple fracture propagation in closed spacing. Stress shadow effects are induced by opening fractures and affect multiple fracture propagation. For closely spaced multiple fractures growing simultaneously, width of the interior fractures are usually significantly restricted, and length of the exterior fractures are much longer than that of the interior fractures. The exterior fractures receive most of fluid and dominate propagation, resulting in immature development of the interior fractures. Natural fractures could further complicate fracture geometry. When a hydraulic fracture encounters a natural fracture and propagates along the pre-existing path of the natural fracture, fracture width on the natural fracture segment will be restricted and injection pressure will increase, as a result of stress shadow effects from hydraulic fracture segments and additional closing stresses from in-situ stress field. When multiple fractures propagate in naturally fracture reservoirs, complex fracture networks could be induced, which are affected by perforation cluster spacing, differential stress and natural fracture patterns. Combination of our numerical model and diagnostic methods (e.g. Microseismicity, DTS and DAS) is an effective approach to accurately characterize the complex fracture geometry. Furthermore, the physics-based complex fracture geometry provided by our model can be imported into reservoir simulation models for production analysis. / text