Natural Fractures are the primary pathways for fluid migration and production in geothermal reservoirs. Reactivation (and possible propagation) of natural fractures is an important component of stimulation, Slip on these natural fractures induces Microseismicity (MEQs) and the locations of MEQ events reveal aspects of stimulation. This work first presents a Geomechanics-Based Stochastic Analysis of Microseismicity (GBSAM) to quantitatively estimate natural fracture orientations from MEQs recorded during reservoir stimulation. This is achieved by combining geomechanics and geostatistics to better constrain uncertainties in natural fracture orientations. The reservoir response to pore pressure changes is modeled using a finite duration line source to simulate water injection. The mechanism for generating MEQs is based on Mohr Coulomb failure criterion while allowing one natural fracture to generate multiple MEQs. Mahalanobis distance, a type of similarity measure technique, is then used to measure the similarity between the simulated MEQ distribution and the field-observed MEQ distribution to find the best natural fracture distribution that fits the field-observed MEQ data. The natural fracture orientations which correspond to best simulated MEQ distribution in GBSAM are considered to be a realistic realization of the subsurface conditions. As an application, the proposed GBSAM is applied to Newberry EGS, and Fenton Hill HDR stimulations to estimate the natural fracture orientations of these two systems. Results from GBSAM show that the natural fractures dip and dip direction inverted from field-observed MEQ distribution are in good agreement with the field results reported from borehole televiewer.