Natural fractures and fracture networks characteristics significantly influence the success of reservoir stimulations. Subsurface and surface investigations consider various geophysical techniques to probe and measure natural fracture characteristics at depth. However, the relationship between geophysical properties of a rock mass and its hydro-mechanical parameters is not unique. As a result, numerical modeling can be called upon to advance the characterization of fractured reservoir stimulation. To simulate fractured reservoir stimulation (e.g., Enhanced Geothermal System), a fully coupled thermo-hydro-mechanical-seismic (THMS) finite element model with a dynamic 3D fracture network capability is further developed which incorporates the processes of fracture flow, poroelastic deformation, shear dilation, fracture network propagation and induced seismicity. A dual-scale semi-deterministic fracture network is generated using data from image logs, foliations/micro-fractures, intact and fractured core. Fracture network propagation under injection is considered using a rapid analytical approach. The induced stress field (stress shadow) caused by fracture pressurized is quantified using the Eshelby solution. The numerical model allows for multiple seismic events to occur on and around a fracture. The location, number and magnitude of seismic events are constrained in an effective physics-based, and conceptually rigorous manner. Several verification and illustrative simulations are carried out. The model is also applied to a field case at intermediate scale (EGS Collab Testbed Experiment 1). The results demonstrate the effectiveness of the coupled thermo-hydro-mechanical-seismic model with a dynamic 3D fracture network in fractured reservoir stimulation and shed light on stimulation mechanisms and potential permeability creation and induced seismicity. The spatial distribution of natural fracture roughness which plays a significant role in the patterns of seismic and aseismic slip on the fractures is considered in this dissertation. A novel approach is used to account for the spatial distribution of roughness on natural fractures. The approach utilizes the concept of Joint Roughness Coefficient (JRC) to obtain a shear stiffness distribution on the fracture surface. The spatial variability of the JRC is obtained through the concept of stochastic partitioning of the fracture plane which can be calibrated using laboratory and field data. The spatial distribution of JRC is then used to map the spatial distribution of fracture shear stiffness. A critical velocity as a function of fracture shear stiffness is defined to distinguish between the seismic/aseismic states. The technique is employed to simulate acoustic emission (AE) distribution during shearing of a rough fracture at laboratory scale. To represent a large-scale real fracture roughness, synthetic anisotropic roughness is used. The synthetic anisotropic roughness is described by three parameters, including two fractal dimensions in perpendicular directions along the fracture, and amplitude elevation of scanning points on the fracture surface. These parameters are used as input for estimating JRC of a synthetic anisotropic roughness. A quadratic polynomial equation is derived to link the Hurst exponent with JRC covering its full spectrum (0-20). Such JRC- Hurst exponent relation is incorporated into the coupled thermo-hydro-mechanical-seismic model with a dynamic 3D fracture network to simulate a field case at a large scale (Utah FORGE). Also, the impact of cyclic injection on the stimulation is studied. Finally, a Geomechanics-Based Stochastic Analysis of Microseismicity (GBSAM) is developed and applied to quantitatively estimate the natural fracture orientations in a reservoir from MEQs (induced seismicity) recorded during reservoir stimulation. This is achieved by combining geomechanics and geostatistics to better constrain uncertainties in natural fracture orientations. The GBSAM is applied estimate the natural fracture orientations (dip and dip direction) in two case studies namely, the Newberry EGS and the Fenton Hill HDR. Results from the GBSAM are in very good agreement with results from boreholes image logs or previous studies in those case studies.