[Truncated abstract] Rocks, as the host to natural chains of coupled thermal, hydraulic and mechanical processes, are heterogeneous at a variety of length scales, and in their mechanical properties, as well as in the hydraulic and thermal transport properties. Rock heterogeneity affects the ultimate hydro-carbon recovery or geothermal energy production. This heterogeneity has been considered one important and difficult problem that needs to be taken into account for its effect on the coupled processes. The aim of this thesis is to investigate the effect of rock heterogeneity on multi-physical processes. A fully coupled finite element model, hereinafter referred to as a porosity-based model (PBM) was developed to characterise the thermal-hydraulic-mechanical (THM) coupling processes. The development of the PBM consists of a two-staged workflow. First, based on poromechanics, porosity, one of the inherent rock properties, was derived as a variant function of the thermal, hydraulic and mechanical effects. Then, empirical relations or experimental results, correlating porosity with the mechanical, hydraulic and thermal properties, were incorporated as the coupling effects. In the PBM, the bulk volume of the model is assumed to be changeable. The rate of the volumetric strain was derived as the difference of two parts: the first part is the change in volume per unit of volume and per unit of time (this part was traditionally considered the rate of volumetric strain); and the second is the product of the first part and the volumetric strain. The second part makes the PBM a significant advancement of the models reported in the literature. ... impact of the rock heterogeneity on the hydro-mechanical responses because of the requirement of large memory and long central processing unit (CPU) time for the 3D applications. In the 2D PBM applications, as the thermal boundary condition applied to the rock samples containing some fractures, the pore pressure is generated by the thermal gradient. Some pore pressure islands can be generated as the statistical model and the digital image model are applied to characterise the initial porosity distribution. However, by using the homogeneous model, this phenomenon cannot be produced. In the 3D PBM applications, the existing fractures become the preferential paths for the fluid flowing inside the numerical model. The numerical results show that the PBM is sufficiently reliable to account for the rock mineral distribution in the hydro-mechanical coupling processes. The applications of the statistical method and the digital image processing technique make it possible to visualise the rock heterogeneity effect on the pore pressure distribution and the heat dissipation inside the rock model. Monitoring the fluid flux demonstrates the impact of the rock heterogeneity on the fluid product, which concerns petroleum engineering. The overall fluid flux (OFF) is mostly overestimated when the rock and fluid properties are assumed to be homogeneous. The 3D PBM application is an example. As the rock is heterogeneous, the OFF by the digital core is almost the same as that by the homogeneous model (this is due to that some fractures running through the digital core become the preferential path for the fluid flow), and around 1.5 times of that by the statistical model.