Fracture Network Model Research Articles

Implementation of manifold coverage for 3D rock fracture network modeling and its application in rock permeability prediction

Fracture network modeling is an important prerequisite for studies on rock mechanical behaviors and permeability properties. At present, fracture patterns in most discrete fracture network (DFN) models are expressed by Euclidean geometries, which leads to distortion in describing complex fractured rocks. In order to solve this problem, a block texture synthesis method is proposed in this paper to generate multi-resolution two-dimensional (2D) DFN models for fractured coal. The results show that the generated multi-resolution 2D DFN models largely retain the fracture texture structure in the original exemplar. In addition, the concept of manifold coverage is adopted to optimize the 2D fracture image, and the three-dimensional (3D) texture-based DFN models are generated by the solid texture synthesis method from the optimized image. The results indicate that the 3D texture-based DFN models have good fidelity of natural fractures in fracture sets, orientation, density, roughness, fluctuation, aperture and connectivity from the perspective of geometric appearance. Comparison between the 3D DFN model with texture coverage and random coverage shows good consistency in fracture proportion, geometric pattern, and permeability properties. Therefore, the 3D DFN models based on random coverage provide a new idea for modeling fractured rocks and lattice models in numerical simulation of rock mechanical behaviors and permeability properties by various methods, such as RFPA2D/3D, FLAC2D/3D, and Lattice Boltzmann Method (LBM).

Identification of the Potential Critical Slip Surface for Fractured Rock Slope Using the Floyd Algorithm

A rock slope can be characterized by tens of persistent discontinuities. A slope can be massive. The slip surface of the slope is usually easier to expand along with the discontinuities because the shear strength of the discontinuities is substantially lower than that of the rock blocks. Based on this idea, this paper takes a jointed rock slope in Hengqin Island, Zhuhai as an example, and establishes a three-dimensional (3D) model of the studied slope by digital close-range photogrammetry to rapidly interpret 222 fracture parameters. Meanwhile, a new Floyd algorithm for finding the shortest path is developed to realize the critical slip surface identification of the studied slope. Within the 3D fracture network model created using the Monte Carlo method, a sequence of cross-sections is placed. These cross-sections containing fractures are used to search for the shortest paths between the designated shear entrances and exits. For anyone combination of entry point and exit point, the shortest paths corresponding to different cross-sections are different and cluttered. For the sake of safety and convenience, these shortest paths are simplified as a circular arc that is regarded as a potential slip surface. The fracture frequency is used to determine the probability of sliding along a prospective critical slip surface. The potential slip surface through the entrance point (0, 80) and exit point (120, 0) is identified as the final critical slip surface of the slope due to the maximum fracture frequency.

Numerical study of the mud loss in naturally fractured oil layers with two-phase flow model

Drilling in the naturally fractured formations always encounters severe mud loss, which endangers the drilling safety. In oil & gas layers, a larger quantity of mud loss has great adverse effects on later oil and gas production. Therefore, it makes the knowledge of mud loss very important, and after that, it is instructive to control the mud loss in the oil layer. In this paper, we first consider the mud loss process using a two-phase flow model. A discrete fracture network model is employed to describe the fluid transfer between fractures and matrix pores. And different relative permeability curves and capillary pressure functions are used in matrix pores and natural fractures. The finite element method and “upwind” scheme are employed to deduce the numerical discretized formulas. The correctness of our model is verified with the published literature. Then a sensitivity analysis is performed to study the laws of mud loss in naturally fractured oil-wet formation. Compared with single-phase flow models, the mud loss rate of two-phase flow model is lower. The two-phase flow model is used to identify the most influential factors on mud loss and predict the distribution of water saturation in oil formation during the drilling process. We find that: since the pressure increase zone is much larger than the water saturation increase zone, the size of the water invasion zone is unknown in the single-phase flow model. The most important influential factors of mud loss in the oil layer are fracture connection, drilling fluid density and matrix permeability. Once the natural fractures are connected to the wellbore, the mud loss rate is much higher than that when the natural fractures are not connected to the wellbore. The capillary pressure has little effect on the mud loss. The larger the capillary pressure, and the lower the mud loss rate. The model in this paper is instructive to learn the law of mud loss and contamination range of drilling fluids in oil layers.

Lithology-Controlled Hydrodynamic Behaviour of a Fractured Sandstone–Claystone Body in a Radioactive Waste Repository Site, SW Hungary

The fracture network modelling and hydrogeological assessment were performed in an 845 m deep borehole of the potential high-level waste repository formation and its caprock. The geometry of the fracture network was simulated using the discrete fracture network (DFN) modelling method, which is based on the geometric characteristics of the individual fractures. The hydrogeological evaluation was based on changes in porosity and permeability along the borehole using flow zone indicator (FZI) values that denote hydraulic flow units (HFU) within the rock body. Fracture network characteristics and hydrogeological features are mainly determined by the wellbore lithology, which can be divided into three zones. The sandstone body was intersected in the upper 300 m of the borehole, which forms a single HFU. The second zone was developed along with the transition zone between the sandstone and the underlying claystone bodies. Here the predominant rock type is claystone, but the characteristics of the fracture network are distinctly different from the deeper parts of this rock body. Below 400 m is the third zone, where distinct and extensive HFU-s could not form, probably due to different water–rock interaction processes that could have changed the porosity and permeability from point to point in the claystone.

Discrete fracture network (DFN) modelling of a high-level radioactive waste repository host rock and the effects on its hydrogeological behaviour

Fracture network modelling and a hydrological evaluation were performed in a well more than 900 m deep that penetrated the Boda Claystone Formation, a potential host rock for high-level nuclear waste disposal facilities in Hungary. The fracture network geometry was generated with a discrete fracture network algorithm, in which the permeability and porosity of the system can be calculated if the aperture of the fractures is known. The hydrological aperture of the fractures was estimated via an aperture calibration based on a comparison of the measured and modelled permeability values. Flow zone indices were calculated for numerous sections along the well, designating hydraulic units, in which the fluid flow-controlling properties are internally even.Based on the fracture network geometry and the hydraulic flow units, most parts of the well behave uniformly, while three narrow zones differ significantly. The first zone is located in the upper 100 m of the well probably formed due to weathering. The second zone is located at approximately 400 m, where a large-scale structural boundary is presumed. In the third zone at 700 m, a lithological change greatly affects the hydrological properties, but the influence of tectonic processes cannot be ruled out.

Enhancing acid fracture design in carbonate formation using a dynamic up-scaling procedure to convert discrete fracture network to dual continuum

For a low-permeability carbonate formation, the acid fracture process is simulated through coupling a commercial acid fracture simulator (GOHFER) to a finite volume reservoir simulator (IMEX). Unlike LGR (Local grid refinement) approach that suffers from severe convergence problems, a dynamic up-scaling procedure is employed to convert the discrete fracture network (DFN) model into a dual continuum model for our simulation. In this paper, multiple simulations are used to optimize the acid fracture schedule parameters, such as fluid volume, flow rate, perforation location, number of injection steps, and acid type, in order to maximize the effective fracture length. For four points perforation as the relevant operation case in this work, acid fracturing consumes 44 Mgal acid volume at a rate of 25 bpm at a 20 ft interval, resulting in a 3.5-fold increase in well productivity index (PI). Based on experimental relative permeability data, two models of fine and medium grid sizes are employed to be compared with the LGR method. The results have shown that these models are 22 and 54 times faster than LGR, however, the latter has a deviation of less than 5% which is quite acceptable.

Fracture Characteristics and Reservoir Inhomogeneity Prediction of the Gaoyuzhuang Formation in the Xiong’an New Area: Insights From a 3D Discrete Fracture Network Model

The Xiong’an new area in northern China is rich in geothermal resources. The Gaoyuzhuang Formation in the Proterozoic Jixian system is an important recently discovered geothermal reservoir. The Gaoyuzhuang reservoir has been affected by multiple stages of tectonic movement. The fractures are very well developed, resulting in strong heterogeneity in the reservoir’s porosity and permeability. Few studies have been conducted on the fractures and heterogeneity of the Gaoyuzhuang reservoir. In this study, interpretation of data from image logs was used to summarize the characteristics of the fractures, including the fractures’ strikes, dips, and lengths. The permeability distribution of the reservoir in the vicinity of the borehole was predicted. The hydro-thermal (TH) coupling model was used to reproduce the process of pumping tests, and the simulation results were found to be in good agreement with the field test data. In addition, the relationship between the fracture aperture and length was obtained, dmax=4.2×10−5⋅L0.5. The results of this study provide data and technical support for further reservoir research and evaluation of the geothermal resources of the Gaoyuzhuang reservoir.

A Rough Discrete Fracture Network Model for Geometrical Modeling of Jointed Rock Masses and the Anisotropic Behaviour

The geometry of the joint determines the mechanical properties of the rock mass and is one of the key factors affecting the failure mode of surrounding rock masses. In this paper, a new rough discrete fractures network (RDFN) characterization method based on the Fourier transform method was proposed. The unified characterization of the complex geometric fracture network was achieved by changing the different Fourier series values, which further improved the characterization method of the RDFN model. A discrete element numerical calculation model of the complex RDFN model was established by combining MATLAB with PFC code. Numerical simulation of the anisotropic mechanical properties was performed for the RDFN model with a complex joint network. Based on the results, the geometry of the joint network has a significant influence on the strength and failure patterns of jointed rock masses. The failure modes of the opening are highly affected by the orientation of the fracture sets. The existence of the rough fracture sets could influence the failure area of different excavation situations. The study findings provide a new characterization method for the RDFN model and a new characterization approach for stability analysis of complex jointed rock masses.

Dynamic fluid transport property of hydraulic fractures and its evaluation using acoustic logging

The existing acoustic logging methods for evaluating the hydraulic fracturing effectiveness usually use the fracture density to evaluate the fracture volume, and the results often cannot accurately reflect the actual productivity. This paper studies the dynamic fluid flow through hydraulic fractures and its effect on borehole acoustic waves. Firstly, based on the fractal characteristics of fractures observed in hydraulic fracturing experiments, a permeability model of complex fracture network is established. Combining the dynamic fluid flow response of the model with the Biot-Rosenbaum theory that describes the acoustic wave propagation in permeable formations, the influence of hydraulic fractures on the velocity dispersion of borehole Stoneley-wave is then calculated and analyzed, whereby a novel hydraulic fracture fluid transport property evaluation method is proposed. The results show that the Stoneley-wave velocity dispersion characteristics caused by complex fractures can be equivalent to those of the plane fracture model, provided that the average permeability of the complex fracture model is equal to the permeability of the plane fracture. In addition, for fractures under high-permeability (fracture width 10~100 μm, permeability ~100 μm2) and reduced permeability (1~10 μm, ~10 μm2 as in fracture closure) conditions, the Stoneley-wave velocity dispersion characteristics are significantly different. The field application shows that this fluid transport property evaluation method is practical to assess the permeability and the connectivity of hydraulic fractures.

Advancements towards DFKN modelling: Incorporating fracture enlargement resulting from karstic dissolution in discrete fracture networks

We present a multiscale approach to develop a deterministic Discrete Fracture Network (DFN) model including fracture enlargement due to karstic dissolution in low-porosity carbonate rocks. The result is a DFKN (Discrete Fracture and Karst Network) model from which it is possible to quantify, in a practical manner, the impact of fracture enlargement due to karstic dissolution on the fracture network transmissivity. We provide a case study using two large outcrops in the Jandaíra Formation, an intensely karstified carbonate platform of Turonian-Campanian age in the Potiguar Basin, Brazil. The study is based on scanlines collected both in field and from Unmanned Aerial Vehicle (UAV) imagery, allowing to cover three orders of magnitude in length and aperture of fractures. We found that if karstic dissolution is low, aperture-length data plotted on bi-logarithmic graphs can be fitted by straight lines, thus following the aperture-length power law typical of mechanical apertures. However, advanced stages of karstification modifies the aperture-length relationship due to the fact that dissolution causes similar fracture enlargement regardless of fracture length. As a result, the aperture-length curve gradually degenerates to a constant value as the degree of dissolution increases. This behavior is the key to quantitatively taking into account fracture enlargement, by simply adapting an aperture-length curve according to the known stage of karstic dissolution of the modeled carbonate rocks. We show that fracture enlargement can increase the transmissivity of a karstified fracture network by up to five orders of magnitude without changing its connectivity. Our approach can be useful to better understand the behavior of fluid flow in carbonate reservoirs with fracture networks affected by karstic dissolution.

Sensitivity analysis on the effect of natural fractures and injected fluid on hydraulic fracture propagation in a fractured reservoir

In this work, a series of numerical simulations based on an outcrop-data-based discrete fracture network (DFN) model of the shale reservoir in the Ordos Basin, China are conducted, trying to determine the dependence of hydraulic fracturing on the properties of natural fractures and injected fluid. The accuracy of the method is verified through the PKN solution for a homogeneous rock model and a laboratory test for a shale sample obtained from the study area. Parametric studies are performed for different field parameters, including natural fracture aperture, natural fracture frictional coefficient, fluid injection rate and fluid viscosity. The results suggest that the level of fracture complexity decreases as the aperture of natural fractures increase, and increasing the frictional coefficient of natural fractures can improve the hydraulic fracturing efficiency. Higher injection rates could translate into higher injection pressure and more energy available for rock failure near the wellbore, and therefore improve the fracture complexity. However, higher fluid viscosity leads to higher resistance to fracture propagation, and thus reduces the hydraulic fracturing efficiency.

Evolution of Fault-Zone Hydromechanical Properties in Response to Different Cementation Processes

Abstract Progressive cementation and sealing of fault-localized fractures impact crustal mass transport and the recovery of fault strength following slip events. We use discrete fracture network (DFN) models to examine how fracture sealing during end-member cementation mechanisms (i.e., reaction- versus transported-limited cementation) influences the partitioning of fluid flow through time. DfnWorks was used to generate randomized fracture networks parameterized with fracture orientation data compiled from field studies. Single-phase flow simulations were performed for each network over a series of timesteps, and network parameters were modified to reflect progressive fracture sealing consistent with either reaction- or transport-limited crystal growth. Results show that when fracture cementation is reaction-limited, fluid flow becomes progressively channelized into a smaller number of fractures with larger apertures. When fracture cementation is transport-limited, fluid flow experiences progressive dechannelization, becoming more homogeneously distributed throughout the fracture network. These behaviors are observed regardless of the DFN parameterization, suggesting that the effect is an intrinsic component of all fracture networks subjected to the end-member cementation mechanisms. These results have first-order implications for the spatial distribution of fluid flow in fractured rocks and recovery of permeability and strength during fault/fracture healing in the immediate aftermath of fault slip.

Numerical study of the effect of rock anisotropy on stresses around an opening located in the fractured rock mass

The fractured rock mass is a discontinuous, heterogeneous, anisotropic medium in nature, including single, multiple or a network of joints with regular or irregular patterns located between the intact rocks. Anisotropy, as an inherent feature of rock mass, is a factor that causes the mechanical properties of rock to be various in different directions, so it is very important and necessary to be considered in the design of rock engineering structures. This study aims to investigate the effect of rock anisotropy on induced stresses around an opening like a borehole located in the fractured rock mass. For this purpose, the combination of discrete fracture network (DFN) method to construct a geometric model of rock mass, with distinct element method (DEM) to the numerical analysis of stress, was used. All coupled DEM-DFN models were developed based on the field data of the rock and joint system of the Sellafield site in the UK. The discrete fracture network models with different degrees of anisotropy, including a circular opening with two different cross-section areas placed at the center of model, were used for stress analysis. In each type of model, by applying different boundary conditions, the magnitude of induced stresses at selected points around opening is calculated, and then the obtained results are compared with the Kirsch method as a popular closed-form solution. The numerical results of this study showed that the presence of irregular joints with different orientations in the model caused anisotropy in the magnitude of induced stresses obtained around the opening. This is an important issue that the Kirsch method is not able to consider in the calculations. • Fractured rock mass is a discontinuous, heterogeneous and anisotropic medium including network of joints with the stochastic pattern located between intact rocks. • Anisotropy is a factor that causes the mechanical properties of rock to be various in different directions. • Coupled DEM-DFN models are developed to numerical analysis of stress around a borehole. • Presence of irregular joints with different orientations caused anisotropy in the amount of induced stresses around the borehole.

Visualizing and Quantifying Generation, Propagation, and Sweep of Nanocellulose-Strengthened Carbon Dioxide Foam in a Complex 2D Heterogeneous Fracture Network Model

Summary There exist two main issues hampering the wide application and development of carbon dioxide (CO2) foam in conformance improvement and CO2 mobility reduction in fractured systems: (1) instability of foam film under reservoir conditions and (2) uncertainties of foam flow in complex fractures. To address these two issues, we previously developed a series of nanocellulose-strengthened CO2 foam (referred to as NCF-st-CO2 foam), while the primary goal of this work is to thoroughly elucidate generation, propagation, and sweep of NCF-st-CO2 foam in a visual 2D heterogeneous fracture network model. NCF-st-CO2 foam outperformed CO2 foam in reducing gas mobility during either coinjection (COI) or surfactant-alternating-gas (SAG) injection, and the threshold foam quality was approximately 0.67. Foam creation was increased with the total superficial velocity for CO2 foam and almost stayed constant for NCF-st-CO2 foam in fractures during COI. For SAG, large surfactant slug could prevent CO2 from early breakthrough and facilitate foaming in situ. The improved sweep efficiency induced by NCF-st-CO2 foam occurred near the producer for both COI and SAG. Film division and behind mainly led to foam generation in the fracture model. Gravity segregation and override was insignificant during COI but became noticeable during SAG, which caused the sweep efficiency decrease by 3 to 9%. Owing to the enhanced film, NCF-st-CO2 foam enabled mitigation of the gravitational effect, especially around the producer.

QDC-2D: A Semi-Automatic Tool for 2D Analysis of Discontinuities for Rock Mass Characterization

Quantitative characterization of discontinuities is fundamental to define the mechanical behavior of discontinuous rock masses. Several techniques for the semi-automatic and automatic extraction of discontinuities and their properties from raw or processed point clouds have been introduced in the literature to overcome the limits of conventional field surveys and improve data accuracy. However, most of these techniques do not allow characterizing flat or subvertical outcrops because planar surfaces are difficult to detect within point clouds in these circumstances, with the drawback of undersampling the data and providing inappropriate results. In this case, 2D analysis on the fracture traces are more appropriate. Nevertheless, to our knowledge, few methods to perform quantitative analyses on discontinuities from orthorectified photos are publicly available and do not provide a complete characterization. We implemented scanline and window sampling methods in a digital environment to characterize rock masses affected by discontinuities perpendicular to the bedding from trace maps, thus exploiting the potentiality of remote sensing techniques for subvertical and low-relief outcrops. The routine, named QDC-2D (Quantitative Discontinuity Characterization, 2D) was compiled in MATLAB by testing a synthetic dataset and a real case study, from which a high-resolution orthophoto was obtained by means of Structure from Motion technique. Starting from a trace map, the routine semi-automatically classifies the discontinuity sets and calculates their mean spacing, frequency, trace length, and persistence. The fracture network is characterized by means of trace length, intensity, and density estimators. The block volume and shape are also estimated by adding information on the third dimension. The results of the 2D analysis agree with the input used to produce the synthetic dataset and with the data collected in the field by means of conventional geostructural and geomechanical techniques, ensuring the procedure’s reliability. The outcomes of the analysis were implemented in a Discrete Fracture Network model to evaluate their applicability for geomechanical modeling.

Using effective medium theory to calculate permeability of rock with complex fractures

Dewatering in surface and underground mining is highly dependent on the permeability of fractured rocks, which is usually investigated by analytical and numerical methods. Analytical methods are used for simple, regular joint systems, while for complicated discrete fracture networks, numerical methods are more suitable. However, such numerical simulations demand considerable computation resources and are time-consuming for large-scale models. An attractive approach is to perform a rough mapping of the fractures into a lattice to apply effective medium theory. A computational code based on effective medium theory for calculating the overall permeability of fracture networks was developed. Several stochastic discrete fracture network models were generated with the fracture aperture distributed uniformly, correlated and uncorrelated with fracture trace length. The results showed that when the aperture was correlated or uncorrelated with trace length, the calculated mean permeability values by the effective medium theory method agreed well with numerical modelling. Also, approximated representative elementary volume was well fitted by numerical methods in the constant fracture aperture case. However, the calculated representative elementary volume by the effective medium theory method was smaller than the numerical estimation in the other two cases.

Numerical modeling of fluid flow in tight oil reservoirs considering complex fracturing networks and Pre-Darcy flow

Multistage hydraulic fracturing for horizontal wells is extensively applied in the exploitation of tight oil reservoirs. Complex fracture networks with different scales are created and hydrocarbon flowing in the matrix demonstrates a Pre-Darcy characteristic due to a strong interaction of rock and fluids in nanoscale pores. The modeling of complex fracture networks and a Pre-Darcy effect is required in numerical simulation of tight oil reservoirs. In this paper, stochastic fracture networks are employed to mimic actual fracturing networks according to the corresponding probability density functions (PDFs), and microseismic events are used to constrain the locations of fractures. The numerical simulation for tight oil reservoirs is developed by considering complex fracture networks and the Pre-Darcy effect based on an embedded discrete fracture model (EDFM) framework. Finally, the influences of fracture network properties including isolated fractures, connected fracture networks and the Pre-Darcy effect on well production are investigated via integrating the developed fracture network generation and tight oil numerical simulator. The results indicate that for isolated fractures with no direct connection with a wellbore, there is little influence on well performance. The influence becomes more obvious when isolated fractures connect to a wellbore, and an increase in fracture length and conductivity can enhance well performance especially when fractures lie off streamline lines. Similarly, for a connected fracture network, an increase in fracture network properties (number, length and conductivity) can promote well performance and fractures perpendicular to a horizontal wellbore can generate a maximal flow rate compared with other fracture orientations. The fracture network connectivity plays an important role in the fluid flow, and there are good relationships between well production and two proposed parameters, which can be treated as candidates to characterize the flow characteristics in fractured reservoirs.

Effect of fracture network on water injection huff-puff for volume stimulation horizontal wells in tight oil reservoir: Field test and numerical simulation study

Horizontal wells with volume fracturing stimulation technique is an effective technique to improve the single well oil production in tight oil reservoir, in which contact area of wellbore and formation is enlarged and large scale fracture network to improve the conductivity is formed. However, water injection is commonly avoided for fractured reservoirs because of the water channeling problem, which leads to the rapid decline of oil production. Water injection huff-puff is proved to be an effective method to improve the formation pressure and oil production. In the previous studies, researchers pay more attention on the huff-puff well, and the production of adjacent wells and its influencing factors are not taken seriously. Herein, the water injection huff-puff test is carried out in Huanjiang oilfield. Performance of both huff-puff well and adjacent wells are improved with oil production increased by 168.29 m3 and 876.82 m3, respectively. Two effect modes for adjacent wells are summarized considering the fracture network, well position and water breakthrough type. Moreover, reservoir numerical simulation is performed to investigate the influence of fracture network on the performance of water injection huff-puff technique. Two kinds of fracture network models of main fracture-connected model and main fracture-disconnected model are established. The results indicate that fracture network influences the performance of water injection huff-puff. For fracture connected model, shutting down adjacent wells during water injection period is an effective method to improve oil production in the whole field due to the fact that water channeling will occur leading to low formation pressure replenishing efficiency and sharp water cut rising of adjacent wells. For fracture disconnected model, formation pressure near the huff-puff well is significantly improved and so as the oil production. Water displacement effect is obvious for adjacent wells. The oil production of adjacent wells increases obviously without water cut rising sharply. In this model, shutting down adjacent wells during water injection is not recommended considering the long well occupation time.

Review of pseudo-three-dimensional modeling approaches in hydraulic fracturing

In the field of hydraulic fracture modeling, the pseudo-three-dimensional (P3D) approach is an efficient and practical computational tool serving as a compromise between two-dimensional and planar three-dimensional models. This review discusses the P3D modeling approach from its early developmental stage in the 1980s to the present. The evolution of P3D modeling is drawn over time based on the major differences in the governing formulation and assumptions considered by each model. The problems of equilibrium height growth and vertical viscous fluid resistance (i.e., non-equilibrium height growth) emphasize the primary differences among these models. Besides, the P3D-based complex fracture network models for shale oil and gas reservoirs accounting for the interaction between preexisting natural fractures and induced hydraulic fractures are discussed. Finally, in the application section, several simulations are reported to demonstrate the validation of the P3D numerical algorithm by comparing it with the Perkins–Kern–Nordgren (PKN) large and small asymptotic solutions, as well as the effect of time-dependent variable injection rates on the hydraulic fracture propagation. The results showed a good matching between P3D and PKN solutions and a significant effect of the wellbore variable injection rate on the evolution of the fracture length.

A particle-tracking formulation of advective–diffusive heat transport in deformable fracture networks

Modeling heat transfer in complex heterogeneous fractured system is key for geothermal energy applications. Discrete fracture network (DFN) modeling is the classical framework to reproduce the advective part of the transport, which is determined by the fracture connectivity and heterogeneity. This approach in general sacrifices the representation of the rock matrix, disregarding both its diffusive heat exchange with the fractures and the effects of its thermo-mechanical deformation on the fracture aperture. Here we propose a new semi-analytic formulation that can be implemented in a DFN simulator with particle tracking approach. The contribution of the rock matrix in terms of diffusive heat exchange and thermal contraction/expansion is analytically evaluated, which respectively impact the advective heat transfer and the fracture aperture variation. The method is proved to be accurate and robust. Results from simulations of cold fluid injection show that rock contraction affects the transmissivity, which accelerates the advective transport resulting in a faster recovery of cold fluid at the outlet. The methodology enables investigating the reservoir behavior and optimizing the geothermal performance while keeping the computational effort within reasonable values. This allows exploring the uncertainty in cases when the in situ characterization is poor, which is the spirit of the DFN modeling.