The coupled behaviors of fluid flow, heat transport and geomechanics in fractured porous media are studied. Emphasis is placed on the coupled hydro-mechanical processes involved in hydraulic fracture propagation, thermal circulation in pre-existing fractures, and reactivation of joints due to fluid injection. A 3D fully coupled hydro-mechanical model is developed to simulate hydraulic fracturing. The model is built on a parallel computation framework. Finite element method is utilized to discretize the governing equations. Hydraulic fracture propagation and fluid flow in the created hydraulic fracture are modeled through a special zero-thickness interface element which is developed based on the cohesive zone model (CZM). The 3D model is verified by considering a penny-shaped hydraulic fracture and a KGD hydraulic fracture (in 3D domain) propagating in both the viscosity- and toughness-dominated regimes. Good agreements have been achieved between numerical results and asymptotic analytical solutions with respect to fluid pressure, fracture height, length and width distributions. The effects of some key CZM parameters and the size of interface element on modeling of hydraulic fracturing are investigated. Using the hydro-mechanical model, height growth of hydraulic fractures in layered formations is simulated. First, the model is validated through a laboratory hydraulic fracturing experiment in the presence of stress contrast. Through explicit modeling of the injection wellbore, the compressibility effects of the wellbore are demonstrated. Comparing the numerical results to those obtained in the laboratory experiment, good agreements in the distribution of fracture aperture, injection pressure and fracture footprint are achieved. Then, numerical analyses are performed to investigate the impacts of in-situ stress contrast, modulus contrast and formation ductility on hydraulic fracture height growth. The interaction between hydraulic fractures and discontinuities is analyzed. The nonlinear mechanical behavior of frictional sliding along interface surfaces is considered. Typical loading paths demonstrated through numerical examples exhibit the capability of the model to simulate different working conditions. Since discontinuities are explicitly simulated through the use of interface element, detailed descriptions of the deformation processes are revealed. For example, information related to aperture opening/sliding and stress distribution along the discontinuities is obtained in the simulations. Numerical simulations of height growth in a layered rock with the existence of a formation interface are performed. The model is first evaluated through the commonly used Renshaw and Pollard’s criterion. Then laboratory experiments on fracture-discontinuity interaction under triaxial-stress conditions are studied. Numerical results match well with those predicted through theoretical formulations and with those observed in laboratory. Typical processes associated with fracture-discontinuity interaction are reveled. With the existence of a horizontal interface, the influence of modulus contrast and stress contrast on hydraulic fracture height growth is analyzed. It is found that the combined effects of material properties of rock, mechanical properties of interfaces, and in-situ stress distribution could effectively inhibit the height growth of hydraulic fractures. Based on the theory of thermo-poroelasticity, together with the equations describing fluid flow and heat transport in fractures, and the equation for fluid flow in wellbore, the coupled processes of fluid flow, heat transport and geomechanics are studied. The thermo-hydro-mechanical behaviors of a fracture/matrix system are modeled through a special zero-thickness interface element. The constitutive law of the interface element is built based on the cohesive zone model (CZM), which is suitable for simulating both tensile and shear failures. Fracture flow is formulated through the commonly used “cubic law”. The fluid flux exchange between the fractures and the surrounding permeable rock matrix is determined by a fluid transfer coefficient. A convective heat transfer coefficient is introduced in the equation governing heat transport in the fractures. Numerical analyses are performed to verify the model, to illustrate some fundamental phenomena, and to provide some applications to laboratory injection and circulation experiments to further validate the model. Numerical simulations reveal the role of mechanical, hydraulic and thermal properties and the coupled processes in the experiments. In addition, the stress and pore pressure distributions during hydraulic stimulation in a heterogeneous poroelastic rock are studied. Results indicate that the stress- and pore pressure distributions are more complex in a heterogeneous reservoir than in a homogeneous one. The spatial extent of stress reorientation during hydraulic stimulations is a function of time and is continuously changing due to the diffusion of pore pressure in the heterogeneous system. In contrast to the stress distributions in homogeneous media, irregular distributions of stresses and pore pressure are observed. Due to the change of material properties, shear stresses and non-uniform deformations are generated. The induced shear stresses in heterogeneous rock cause the initial horizontal principal stresses to rotate out of horizontal planes. At last, an on-site hydraulic stimulation at the Phase I Fenton Hill geothermal reservoir is studied. Four scenarios are proposed to analyze the mechanisms involved in repeated injection-venting experiments. It is found that the stiffness of joint, a key parameter used in aperture calculation, controls the flow-back volume and trapped fluid pressure during venting operations. Considering the size dependent characteristic of joint stiffness and hysteresis behaviors observed during injection and venting, a parameter related to stiffness is gradually changed after each injection-venting treatment. In this way, the numerical simulations yield results close to those observed in the field test.