Tight gas and shale oil play an important role in energy security and in meeting an increasing energy demand. Hydraulic fracturing is a widely used technology for recovering these resources. The design and evaluation of hydraulic fracture operation is critical for efficient production from tight gas and shale plays. The efficiency of fracturing jobs depends on the interaction between hydraulic (induced) and naturally occurring discrete fractures. In this work, a coupled reservoir-fracture flow model is described which accounts for varying reservoir geometries and complexities including non-planar fractures. Different flow models such as Darcy flow and Reynold's lubrication equation for fractures and reservoir, respectively are utilized to capture flow physics accurately. Furthermore, the geomechanics effects have been included by considering a multiphase Biot's model. An accurate modeling of solid deformations necessitates a better estimation of fluid pressure inside the fracture. The fractures and reservoir are modeled explicitly allowing accurate representation of contrasting physical descriptions associated with each of the two. The approach presented here is in contrast with existing averaging approaches such as dual and discrete-dual porosity models where the effects of fractures are averaged out. A fracture connected to an injection well shows significant width variations as compared to natural fractures where these changes are negligible. The capillary pressure contrast between the fracture and the reservoir is accounted for by utilizing different capillary pressure curves for the two features. Additionally, a quantitative assessment of hydraulic fracturing jobs relies upon accurate predictions of fracture growth during slick water injection for single and multistage fracturing scenarios. It is also important to consistently model the underlying physical processes from hydraulic fracturing to long-term production. A recently introduced thermodynamically consistent phase-field approach for pressurized fractures in porous medium is utilized which captures several characteristic features of crack propagation such as joining, branching and non-planar propagation in heterogeneous porous media. The phase-field approach captures both the fracture-width evolution and the fracture-length propagation. In this work, the phase-field fracture propagation model is briefly discussed followed by a technique for coupling this to a fractured poroelastic reservoir simulator. We also present a general compositional formulation using multipoint flux mixed finite element (MFMFE) method on general hexahedral grids with a future prospect of treating energized fractures. The mixed finite element framework allows for local mass conservation, accurate flux approximation and a more general treatment of boundary conditions. The multipoint flux inherent in MFMFE scheme allows the usage of a full permeability tensor. An accurate treatment of diffusive/dispersive fluxes owing to additional velocity degrees of freedom is also presented. The applications areas of interest include gas flooding, CO₂ sequestration, contaminant removal and groundwater remediation. / text