Multi-stage hydraulically fractured wells are applied widely to produce unconventional resource plays. In naturally fractured reservoirs, hydraulic fracture treatments may induce complex fracture geometries which cannot be modeled accurately and efficiently using Cartesian and corner-point grid systems or standard dual porosity approaches. The interaction of hydraulic and naturally occurring fractures almost certainly plays a role in ultimate well and reservoir performance. Current simulation models are unable to capture the complexity of this interaction. Generally speaking, our ability to detect and characterize fracture systems is far beyond our capability of modeling complex natural fracture systems. In order to evaluate production performance in these complex settings using numerical simulation, fracture networks require advanced meshing and domain discretization techniques. This paper investigates these issues by developing natural fracture networks using fractal-based techniques. Once a fracture network is developed we demonstrate the feasibility of gridding complex natural fracture behavior using optimization-based unstructured meshing algorithms. It can then be demonstrated that natural fracture complexities such as variable aperture, spacing, length and strike can be simulated. This new approach is a significant step beyond the current method of dual porosity simulation which essentially negates the sophisticated level of fracture characterization pursued by many operators.We use currently established code for fractal discrete fracture network (FDFN) models to build realizations of naturally fractured reservoirs in terms of stochastic fracture networks. From outcrop, image log and core analysis, it is possible to extract fracture fractal parameters pertaining to aperture, spacing, length distribution, including center distribution as well as a fracture strike. Then these parameters are used as input variables for the FDFN code to generate multiple realizations of fracture networks mimicking fracture clustering and randomly distributed natural fractures. After incorporating hydraulic fractures, complex fracture networks are obtained for further reservoir domain discretization.In order to discretize the complex fracture networks, a new mesh generation approach is developed to conform to non-orthogonal and low-angle intersections of extensive-clustered discrete fracture networks with non-uniform aperture distribution. Optimization algorithms are adopted to reduce highly skewed cells, and to ensure good mesh quality around fracture tips, intersections and regions of extensive fracture clustering. Moreover, Local Grid Refinement (LGR) is implemented with a pre-defined distance function to control cell sizes and shapes around and far away from fractures. Natural fracture spacing, length, strike and aperture distribution are explicitly gridded thus introducing a new simulation approach that is far superior to dual porosity simulation. Finally initial sensitivity studies are performed to demonstrate both the ca...