Simulation of Naturally Fractured Reservoirs. State of the Art

Résumé — Simulation des réservoirs naturellement fracturés. État de l’art : Partie 1 – Mécanismes physiques et formulation des simulateurs — L’utilisation des simulateurs de réservoir fracturé aide l’ingénieur de réservoir à mieux comprendre les principaux mécanismes physiques, à choisir le procédé de récupération le mieux adapté et à l’optimiser. Des progrès sensibles ont été réalisés depuis les pre-mières publications sur le concept double-milieu dans les années soixante. Cet article présente les techniques actuelles de modélisation utilisées dans les simulateurs industriels. Après une description des pro-cédés de récupération et des principaux mécanismes physiques associés, suit un historique de la simulation des réservoirs. La formulation mathématique générale des simulateurs double-milieu est ensuite décrite. L’article se termine sur une présentation de la simulation numérique de l’écoulement dans les fractures et les failles et de la modélisation des puits. La formulation des échanges matrice-fracture ainsi que les techniques d’analyse des incertitudes et de calage d’historique contraint par la géologie et l’écoulement sont traitées dans un second article, associé à celui-ci, Partie 2 : Échanges matrice-fracture et spécificités des études numériques. Abstract — Simulation of Naturally Fractured Reservoirs. State of the Art: Part 1 – Physical Mechanisms and Simulator Formulation — Use of fractured reservoir simulators can help reservoir engineers in the understanding of the main physical mechanisms and in the choice of the best recovery process and its optimization. Significant progress has been made since the first publications on the dual porosity concept in the sixties. This paper presents the current techniques of modeling used in industrial simulators. Following a description of the recovery processes and of the main physical mechanisms involved, a history of the fractured reservoir simulation is presented. Then the general mathematical formulation of dual porosity simulators is described. The paper ends with a presentation of the numerical simulation of flow in fractures and faults and of well modeling. The matrix-fracture transfer formulations, as well as techniques for uncertainty analysis and geology&flow-constrained history matching, will be addressed in the companion paper, Part 2: Matrix-Fracture Transfers and Typical Features of Numerical Studies.


INTRODUCTION
Petroleum field development has a very high cost and a strong degree of uncertainty.It follows from that the need of numerical simulations to estimate reservoir performance under a variety of production schemes.Continuous technology improvements made reservoir simulations easier and more realistic.Today and since about twenty years, results from reservoir simulations are used for major reservoir development decisions.Reservoir simulation can answer several crucial issues such as the choice of the "best" recovery process (technical and economical) for a given reservoir, the calculation of production profiles and of expenses (operating costs and capital expenditure), and the risk evaluation.Numerical simulations have become a reservoir management tool at all stages of the reservoir life.This is particularly true in the case of fractured reservoirs.
The dual-continuum, or dual-porosity, approach was first proposed by Barenblatt and Zheltov (1960) in order to simulate the flow behavior in fractured porous media.Since this date, this concept has been applied by Warren and Root (1963) to model the transient well test responses of fractured reservoirs.Then, it was adopted as the underlying model of all industrial fractured reservoir simulators, including threedimensional multiphase compositional ones.Continuous progresses have been made, especially in the modeling of matrix-fracture transfers.Dual porosity simulators are widely used in the petroleum industry for field-scale simulations of recovery processes in fractured reservoirs.
Recently, major advancements were achieved in the simulation of the behavior of fractured reservoirs with the development of efficient methodologies for integrated reservoir studies (Bourbiaux et al., 1997(Bourbiaux et al., , 1998(Bourbiaux et al., , 2002(Bourbiaux et al., , 2005;;Sabathier et al., 1998;Cacas et al., 2001).The main steps of an integrated workflow are the following: -constrained modeling of the fracture network, based on the analysis, the interpolation and extrapolation of fracture information acquired in wells and seismic surveys, sometime completed by outcrop analogue data; -characterization of the hydrodynamic properties of this natural network from flow-related data; -choice of a flow simulation model suited to the role played by fractures and faults at various scales and involving upscaled parameters derived from the flow-calibrated geological fracture model; -simulation of reservoir flow behavior on the basis of a physical assessment of multiphase flow mechanisms prevailing in transfers within and between media.
The first three steps are addressed in a companion paper (Bourbiaux, 2010), which discusses the specific framework of fractured reservoir simulation, as an introduction to this paper.The present paper will focus on the fourth step and present the state of the art in the domain of dynamical simulation of fractured reservoirs with dual-porosity and dualporosity dual-permeability models.The first section describes the recovery processes applied to fractured reservoirs and the main drive mechanisms involved.The second section presents a history of the development of fractured reservoir simulators.The third section describes the features of dual-medium simulators -mathematical formulation of the physics, alternative approaches in fractured reservoir simulation, simulation of specific phenomena in fractures, simulation of conductive faults and the different issues of well modeling.Other topics like modeling of matrix-fracture transfers, modeling of geomechanics effects, quantification of uncertainties and history matching are addressed in the companion paper, Part 2: Matrix-Fracture Transfers and Typical Features of Numerical Studies.

DESCRIPTION OF RECOVERY PROCESSES IN FRACTURED RESERVOIRS
There are major differences between recovery performance of fractured and non-fractured reservoirs.The high contrast of capillarity between the matrix and the fractures is the main cause of these differences.A review of the recovery mechanisms in fractured reservoirs can be found in (Firoozabadi, 2000).As pointed out by this author, one of the main characteristics of fractured reservoirs is high rate wells in the early life of the field, due to the high effective single-phase permeability of the combined matrix-fracture porous media.The flow behavior contrast between fracture and matrix media is emphasized under two-or three-phase flow conditions.Considering gas-oil gravity drainage for instance, capillary forces hinder positive gravity displacement effects on matrix oil recovery; moreover the oil drained from the matrix to the fracture can partially or totally re-imbibe neighboring blocks under the effect of capillary forces.Finally gas-oil displacement turns out to be a complex recovery mechanism as it is governed by capillary, gravity and viscous forces and affected by compositional effects including transfer between phases and diffusion.Depending on their matrix block sizes and matrix/fissure permeabilities, fractured reservoirs can be produced using several recovery processes, primary recovery, gas drive (gas cap expansion and solution gas drive combined with possible convection phenomena), waterflood, miscible or immiscible gasflood and enhanced oil recovery methods.The main drive mechanisms involved in these processes are schematized in Figure 1.Depending of the nature of the recovery process, one or several mechanisms, among imbibition, water drive, gravity drainage, reimbibition and diffusion, can contribute to production.Although it motivated research studies in the laboratory, the detailed physics of these mechanisms is not reviewed and only their expression and consequences for modeling are reported herein.

Primary Recovery
When the production is initiated, the total compressibility of the fluids and the fractured rock is one of the key factors that affect recovery performance.Knowledge of formation compressibility can be very important when the oil is highly undersaturated (Firoozabadi, 2000).High compressibility allows economical depletion of fractured reservoirs with no matrix porosity.The duration of this initial phase depends on the magnitude of the pressure gap between the initial reservoir pressure and the saturation (or bubble point) pressure.
For gas fractured reservoirs, expansion is a predominant mechanism of production due to the high fluid compressibility.In the case of active aquifer the duration of this pressure depletion phase in the matrix is a function of the height of the gas column above the original gas-water contact.

Gas Cap Expansion
Gas cap expansion is dominated by gravity drainage.For thick and/or high-relief fractured reservoirs consisting of matrix blocks with a large effective height and/or high permeabilities, gas cap expansion can be efficient thanks to a high gravity head and a sufficient kinetics of drainage.Microscopic efficiency of gas displacement is generally higher than that of water because the oil phase can spread over water in the presence of gas.This process was observed in the laboratory and was found to lead to very low residual oil saturation (Lange, 1998).

Solution Gas Drive
As depletion goes on, the pressure falls below the bubblepoint pressure in the lowest-pressure regions, that is, in the

Water
Oil Gas

WOC
upper regions of the reservoir and also close to wellbores.Gas bubbles nucleate within the oil phase.When the bubble point pressure is reached within the matrix blocks, gas bubbles appear within the pore network of the matrix.As long as these bubbles grow while remaining immobile, an oil phase somewhat impoverished in gas is expelled from the matrix blocks and conveyed to the production wells.At that precise stage, a decrease in the production Gas-Oil Ratio (GOR) may then be observed if significant.Actually very soon, gas bubbles coalesce and form a mobile phase, for a minimum value of the gas saturation, called the critical gas saturation S gc .As pointed out by Firoozabadi et al. (1992) knowledge of critical gas saturation is important for estimating recovery in a solution gas-drive reservoir.Measured values of the critical gas saturation in the literature range from 2% to 27% PV.For solution gas-drive reservoirs, and particularly for fractured petroleum reservoirs, higher values of critical gas saturation mean higher oil recoveries.When the gas is mobilized, gas drive becomes a prevailing recovery mechanism of the matrix oil in the oil zone.However, the oil-to-gas mobility ratio, that is already unfavorable because of the low gas viscosity, rapidly gets worse because of a rapidly-increasing gas-to-oil relative permeability ratio, due to gas saturation increase.In addition, for reservoirs subjected to convection phenomena (see Sect. 3.3.1),the diffusion of matrix solution gas to the fractures causes a shrinkage of the oil phase within the matrix blocks which further decreases the matrix oil phase saturation and mobility.Eventually, the oil recovery from solution gas drive alone generally remains low.At that time, the evolution of well effluent differs according to the reservoir structure and fracture intensity: -in the case of high-dip, thick and intensely-fractured reservoirs such as Iranian fields, the gas expelled by the matrix blocks can segregate within the fracture network and form a secondary gas-cap.The latter will expand via the fracture network as production goes on, and initiate gas-oil gravity drainage of the gas-surrounded matrix blocks.In that situation, the production GOR stays at a lower value than in conventional reservoirs; -the previous situation typical of well-fractured reservoirs may not occur in the presence of a low-relief reservoir with a poorly-connected and/or low-permeability fracture network.Such reservoirs then behave like conventional ones and rapidly deliver a high-GOR effluent.Such a production behavior calls for a strategy of pressure maintenance, denoted as "secondary recovery", that involves the (re)injection of a fluid phase, either water or gas.Solution gas drive is generally regarded as an ineffective recovery process for fractured reservoirs, except for hard-toproduce fields where other recovery mechanisms, driven by capillarity and/or gravity forces, are ineffective.This is the case for tight, viscous and oil-wet reservoirs where the matrix can neither be imbibed by water nor effectively drained by gas.

Waterflood
Injection of water represents the majority of the projects of secondary recovery compared to gas injection.As the injected water preferentially flows through the fracture network, a high water-saturation boundary condition is established on the matrix blocks.The displacement of oil by water in the matrix medium is then due to three mechanisms: -spontaneous capillary imbibition if the matrix rock is water-wet; -viscous displacement under the pressure gradient generated by fracture flows; -and gravity effects linked to water-oil density difference.
For preferentially water-wet fractured reservoirs with small blocks, oil production by water imbibition is efficient.Indeed, due to capillary forces, water spontaneously invades the matrix blocks thus producing oil by countercurrent and cocurrent flows (Bourbiaux et al., 1990).However, intermediate to oil-wet field situations are often encountered in practice.In such fields, the waterflood efficiency will depend on the possible role played by other recovery mechanisms, including gravity forces and also viscous drive due to fracture flow.For instance, Firoozabadi (2000) quotes intermediate-wettability field examples, such as the Ekofisk chalk field in the North Sea, where water injection turned out to be very efficient although laboratory spontaneous imbibition measurements would predict very little oil recovery.Tang and Firoozabadi (2000) conducted laboratory experiments on Kansas outcrop intermediate-wet chalk samples which highlighted the contribution of viscous forces to water injection efficiency.Beliveau et al. (1993) presented an evaluation study of more than 20 years of waterflood history in the naturally fractured Midale field in Canada, a heavily fractured vuggy limestone overlain by a less fractured chalky dolomite.Waterflood performance is dominated by these oriented vertical fractures, which typically are spaced 1 to 4 ft apart.A dual-porosity simulator was used for modeling the mechanisms of fracture-matrix fluid transfer, i.e. capillary imbibition, gravity drainage, pressure equilibration and, in the case of Midale, viscous displacement due to the low fracture conductivity.Sensitivity runs showed that imbibition effects are directly proportional to capillary pressure and matrix permeability, and inversely proportional to the square of fracture spacing.The simulations also suggested that gravity drainage plays only a minor role in fracture-matrix fluid transfer in this thin reservoir.The different well behaviors were matched by varying only fracture spacing and fracture network conductivity.This suggests that a good characterization of the fracture system made by a multidisciplinary team was critical to understanding the Midale reservoir behavior.Ultimate waterflood recovery was predicted to be 24% OOIP.The study also showed that opportunities exist to increase sweep efficiency through infill drilling, especially with horizontal wells to take advantage of the natural fractures.
After waterflooding a reservoir, a large amount of oil remains trapped in the matrix, both at the pore level of the matrix medium by microscopic capillary trapping effects and, at a large scale, within poorly-swept zones.For instance, the picture of the Midale reservoir after waterflooding (described above) is that of a high-permeability conduit in the lower part of the reservoir overlain by a rich-oil zone.This led to study a tertiary recovery process, and a miscible CO 2 -flood pilot demonstrated that a substantial recovery is obtainable.Field-scale CO 2 -flood recoveries are predicted to be an incremental 20% OOIP owing to the specific Midale geology and favorable gravity effects.
Gas fractured reservoirs are also concerned by waterdrive.Indeed, significant quantities of natural gas are produced from low-permeability carbonate and shale fractured reservoirs (Holditch, 2006).Tight gas carbonate fields are often faced with early water breakthrough in the presence of fractures connected with an active aquifer.Examples of tight naturally-fractured gas reservoirs with active water can be found in the foothills of Alberta and Northeastern British Columbia in Canada.Initial gas recovery stems from the mechanism of expansion of both fracture and matrix gas.In the presence of an active aquifer, subsequent gas recovery involves an imbibition process, which may be very slow.Few papers are the subject of imbibition recovery of gas.Pow et al. (1997) studied imbibition rate and ultimate recovery on representative core obtained from a tight carbonate reservoir in the foothills of Alberta.They concluded that in some cases, the best operating strategy may be to produce this kind of reservoir at the highest rate as possible until the wells water out, and then to shut-in wells to allow gas to re-accumulate.Egermann et al. (2007) performed spontaneous water-gas imbibition experiments on samples from a vuggy carbonate field and showed that water imbibition into vuggy carbonates involves an important gas trapping and the slow imbibition kinetics depends on matrix block size and permeability.

Gas Injection
By comparison with water injection, gas injection is most often a compositional recovery process as the injection gas is rarely in full equilibrium with the in-situ oil.Depending on the oil and injected gas compositions and on the thermodynamic conditions in the reservoir, miscibility of the injected gas with the in-situ oil can be achieved in theory, either through a direct miscibility or through a dynamic process that gradually enriches the in-situ oil or the injected gas in intermediate components until miscibility conditions are reached.However, achieving such favorable displacement conditions in a fractured reservoir requires that the injected gas does not bypass the matrix oil but is driven into the blocks by gravity forces and molecular diffusion combined with thermodynamic transfers between gas and oil phases.The same requirements have to be satisfied in an air injection process to ensure that oxygen consumption by matrix oil is concomitant to gas progression within the fracture network.
In many practical situations, such favorable thermodynamical conditions cannot be satisfied and the gas injection process remains essentially a multiphase process with gravity drainage being the main mechanism driving oil from the matrix into the fracture network.However, diffusion and gasoil thermodynamical exchanges leading to oil swelling or vaporization contribute to an additional recovery of matrix oil.Oil is first recovered by a swelling effect if the oil phase is undersaturated.Gas dissolves into the oil phase, therefore increasing its volume and expelling the oil out the blocks.Later on, when vapor phase saturation is established in the blocks, ongoing flow of fresh injected gas in the surrounding fractures entails the onset of a vaporization process that leads to the extraction of the light components from the oil phase.
Because the liquid phase is always wetting a rock surface, gas-oil capillary forces always contribute to oil retention within matrix blocks and hinder the positive role played by gravity forces.However, their magnitude varies strongly with reservoir pressure, vanishing close to the gas-oil miscibility pressure.Fracture viscous flow effects are less important than for water injection because of the low gas viscosity, and can generally be neglected except in poorly-fractured media.
Gas injection can be carried out using an external gas source like CO 2 , nitrogen, air or re-injecting the produced gas plus an additional make up gas when needed.Daltaban et al. (2002) compared from numerical simulations the efficiency of miscible and immiscible gas injection into the Chuc reservoir in the Gulf of Mexico.This reservoir contains two sets of fractures, namely a network of microfractures due to diagenesis and tectonically-induced large-scale fractures (up to hundreds of meters).Miscible separator-gas injection and immiscible nitrogen injection (N 2 miscibility was as usual not achieved at the reservoir conditions) were studied.Numerical results show that nitrogen injection yields a lower recovery than miscible gas injection, with a recovery of only 16% of the remaining oil compared to 27% with separator gas.This low recovery factor is due to the thermodynamic nature of the displacement process.Nitrogen injection is a vaporizing gas drive but because N 2 -enriched gas flows rapidly in the reservoir and through the fractures, an early gas breakthrough occurs resulting in a premature production decline and low oil recovery.
Miscibility in single-porosity reservoirs can be assessed from well-known experimental and numerical procedures.The 1-D miscible flow process has been shown to be dictated by phase behavior alone and not by rock or other fluid properties.This is not the case for fractured reservoirs where the physics behind miscible injection is more complex and not well known.Uleberg et al. (2002) proposed a method based on compositional reservoir modeling for determination of Minimum Miscibility Pressure (MMP) in fractured reservoirs for multi-dimensional systems.Instead of recording the final recovery, they record all grid block physical properties such as saturation, equilibrium phase compositions, densities and gas-oil interfacial tensions at reported times throughout the simulation run.The gas-oil interfacial tension is used as the miscibility indicator.The method allows studying developed miscibility locally in specific regions of the reservoir and thus is well suited for fractured reservoirs to estimate miscibility in given matrix blocks.The authors concluded that minimum miscibility pressure/enrichment levels for fractured reservoirs are significantly higher than for conventional single-porosity reservoirs, due to the impact of fractures that play as major heterogeneities, and to necessary transfers into the matrix by molecular diffusion.They recommended detailed fine-grid compositional modeling to capture all the complex recovery mechanisms of non-equilibrium gas injection in fractured reservoirs.
Immiscible CO 2 injection is also an efficient EOR process for fractured reservoirs as illustrated by the case of Bati Raman heavy-oil field operated by the Turkish Petroleum Corporation (TPAO) since year 1986.The quantitative understanding of the oil recovery mechanisms and sweep efficiency by CO 2 were attempted through an integrated field study (Combe et al., 1997;Faure et al., 1997).This approach involved laboratory measurements, reservoir characterization and numerical simulation studies at various scales from the core to the full-field scale.Two oil recovery mechanisms are associated with the high amount of CO 2 which can be dissolved into the oil: oil swelling and oil-viscosity reduction.The penetration of CO 2 into the oil-saturated matrix blocks occurs through phase transfer and liquid diffusion processes, the efficiency of which necessitates a sufficiently large contact area between the injected gas and the oil-saturated rock.This contact area was expected to be large in Bati Raman, because of the numerous fissures present in the reservoir, rapidly filled by the injected gas.

Conclusion
The existence of a capillary continuity between blocks (block-to-block transfer function) and the description of the exchanges between the fractures and the matrix blocks (fracture-block transfer function) are critical issues for the efficiency of the three processes described above, gas cap expansion, waterflood and gasflood.Capillary continuity has a great impact on the ultimate oil recovery.Assuming capillary continuity between the matrix blocks, the stack of blocks is drained as a single block having the same height as the stack.Without continuity, each matrix block drains separately until the capillary hold-up height of each individual block is reached at the equilibrium between capillary and gravity forces.Sensitivity runs are often necessary to highlight the main recovery mechanisms acting within a fractured reservoir.Conway et al. (1996) performed such a reservoir study on the A-block reservoir of Gorm Field, a tight chalk oil reservoir subject to both gas and water injection.They showed the effect on the production performance of capillary continuity, oil-water imbibition, gas relative permeability hysteresis and viscous gradients within the fracture network.In this fractured reservoir with low to moderate effective fracture permeabilities, flow-induced pressure gradients within the fracture network can indeed become significant, contrary to typical fractured reservoir situations.

Enhanced Oil Recovery
The efficiency of matrix-fracture transfer mechanisms may vary considerably from one fractured reservoir to another depending on the flow properties of the matrix porous medium and the characteristics of the fracture network.Hence, the oil recovery from fractured fields covers a very large range, thus leaving enhanced oil recovery opportunities for many of them.
Enhanced Oil Recovery (EOR) methods, such as gas injection combined with water injection (Water Alternate Gas), chemical or thermal methods are also feasible methods for producing complex fractured reservoirs.As previously mentioned, waterflood with capillary imbibition is an efficient recovery process for water-wet fractured reservoirs with small blocks.Nevertheless, after waterflooding a reservoir, a large amount of oil remains trapped in the matrix, both at the pore level of the matrix medium by microscopic capillary trapping effects and, at a large scale, within unswept zones.In the case of oil-wet or poorly water-wet fractured reservoirs, one must also mention a "mesoscopic" trapping due to the vertical equilibrium between capillary and gravity forces at the matrix block scale.Moreover, unfavorable conditions, such as a high oil viscosity or an oil-wet matrix can reduce or annihilate the effect of the capillary imbibition mechanism.A way to enhance capillary imbibition is to change the water and oil properties (water or oil viscosity, water-oil interfacial tension) by injecting steam, hot water, polymer or surfactant.Babadagli (2001) studied by laboratory-scale experimentation the performance of those three methods, surfactant, polymer and hot injection in the case of heavy oil.The three methods yielded a higher and faster capillary imbibition of matrix blocks, compared to waterflooding.Two issues should be taken into consideration for field application, the cost efficiency of the process and the upscaling of the phenomenon from laboratory to field scale.Manrique et al. (2007) presented an overview of EOR experiences in carbonate reservoirs in the United States.They concluded that CO 2 flooding was the dominant recovery method and that chemical methods have made a relatively small contribution in terms of total oil recovered.EOR projects by chemical methods were developed essentially during the period between 1980 and 1990, with polymer flooding as the most important method.Energy demand and high energy prices will sooner or later lead to an increase of EOR projects in the future, in comparison with the past two decades.
Polymer or gel injection is more related to wellbore treatment than to EOR.The main objective is to counteract the detrimental impact of fractures on the macroscopic sweep efficiency and to control the injected fluid mobility within the fracture network.Wassmuth et al. (2004) conducted experiments for application to the Weyburn reservoir, showing that, after gel placement in fractured cores, the miscible CO 2 flood could achieve oil recoveries similar to those achieved while flooding non-fractured cores.Field-scale polymer injection cases reported in the literature are related to the period between 1960 and 1990 (Manrique et al., 2007).Hochanadel et al. (1990) mentioned the injection of an anionic polymer in some wells of the Townsend Newcastle Sand Unit fractured reservoir, located in Wyoming, to reduce the channeling effect of fractures.Regarding matrix oil, a low-concentration surfactant injection can yield a significant increase of the matrix oil recovery according to lab-scale experimental results (Chen H.L. et al., 2000;Spinler et al., 2000) but is seldom applied at field scale, due to economic risks as well as the difficulty in finding surfactants still efficient in high salinity and temperature conditions.However, Chen H.L. et al. (2000) mention single and multi-well pilot tests at Yates that showed an oil recovery improvement and a water-oilratio reduction in response to well treatment with a dilute surfactant solution.
Several field applications of steam injection in fractured carbonate reservoirs were reported: Lacq Supérieur, France (Sahuquet et al., 1982), Qarn Alam, Oman (Macaulay et al., 1995;Al-Shizawi et al., 1997;Zellou et al., 2003), Teapot Dome, Wyoming (Doll et al., 1995), South Belridge, California (Johnston et al., 1995), Yates, West Texas (Snell et al., 1999(Snell et al., , 2000)).Most projects are related to medium and light oil reservoirs, 32°API for Teapot Dome, 30°API for Yates and to some heavy oil reservoirs, 21°API for Lacq, 16°API for Qarn Alam.The main conclusions of these steam injection projects are the following: -a reservoir with a strongly fissured character can be efficiently produced using a steamdrive process; -as for any improved oil recovery method, but especially for a steam process, the successful design and implementation of the project depends on an accurate characterization of the fracture system; -a significant additional recovery is obtained by steam injection in medium-oil fractured reservoirs compared to primary and waterflood recoveries.
As pointed out by Hoffman et al. (2003), displacement mechanisms for steam injection into heavy oils are quite different from those into light oils.The heavy oil main recovery mechanisms are thermal expansion and accelerated gravity drainage by viscosity reduction, whereas for light oils, the incremental oil recovery is mainly linked to the vaporization of light hydrocarbon components, and to a lesser extent to thermal expansion and oil viscosity reduction effects on gravity drainage.At late production times and as the distillate bank approaches the producing wells, vaporization becomes the main mechanism.Vaporization results in a very low residual oil saturation.Thermal expansion is an important recovery mechanism in fractured reservoirs due to the thermal conduction which conveys heat in reservoir areas not contacted by steam flow.Dehghani et al. (2001) studied the feasibility of steam injection for three light-oil reservoirs including a vuggy carbonate reservoir.As Hoffman, they concluded that light-oil steam projects considerably differ in their implementation from conventional heavy-oil steamflood because distillation effects allow the use of wider well spacing and greater injection rates.
Water injection may be very efficient in some weakly water-wet chalky fractured reservoirs, as mentioned by Firoozabadi (2000).Many oil-wet naturally fractured reservoirs are nevertheless not prone to water injection because water will not imbibe into the matrix but flow preferentially through the fractures, resulting in very low oil recovery values.Al-Hadhrami and Blunt (2001) proposed to inject steam or hot-water in oil-wet fractured carbonate reservoirs.The steam or water heats the rock, which then undergoes a thermally-induced wettability reversal.Hot water can then spontaneously imbibe into the water-wet rock matrix, resulting in favorable oil recoveries.The process was analyzed, using an analytical solution in 1-D, with application to the Ghaba North field conditions in Oman.The authors showed that 30% oil recovery could be expected with steam or hot waterflooding, compared to only 2% with natural aquifer drive, but no field application was reported.As pointed out by Roosta et al. (2009), wettability alteration of fractured reservoirs towards more water wetness leads to more oil recovery due to capillary imbibition of water into matrix blocks, and prevents the oil re-imbibition into adjacent blocks.Thermally induced wettability alteration has been the subject of several studies, but there remain some uncertainties regarding the effects of temperature on wettability and relative permeabilities and the mechanism of wettability alteration (Ayatollahi et al., 2005).
Among the alternative EOR processes suited for fractured reservoirs, air injection appears as a promising emerging technique, especially for light-oil reservoirs.Turta and Singhal (2001) reviewed the field projects of Air-Injection Processes (AIP) and addressed the reservoir aspects of air injection as an enhanced oil recovery technique for light-oil reservoirs.They noticed that the Low Temperature Oxidation (LTO) process prevails for light-oil fractured reservoirs.They referenced several field projects in dolomite and limestone reservoirs in the Williston basin of North and South Dakota, U.S.A.For instance, air injection in a fractured reservoir was tested in the extensively-fractured CAPA Madison reservoir, North Dakota, where air was injected at the end of a waterflood.After 1.5 years of pattern flooding in this watered-out reservoir, the Air-Oil Ratio (AOR), that is the most important economic parameter, was 20 000 scf/bbl (3562 vol./vol.),twice that of other Williston basin AIP projects, and the process was discontinued.One of the possible reasons for this failure is the very low porosity (11%) of this reservoir.Turta and Singhal (2001) presented preliminary screening criteria for air-injection-based processes.Using these criteria, fractured reservoirs are not considered as good candidates for in-situ Combustion (ISC) or High Temperature Oxidation (HTO).Since the injected air flows preferentially in the high permeability fractures, the air contacts only oil present in these fractures or in their immediate vicinity and the combustion process is not sufficiently sustained.However, Schulte and de Vries (1985) conducted experiments showing that the burning process is governed by diffusion of oxygen from the fractures into the matrix.The main oil production mechanisms were found to be thermal expansion and evaporation with subsequent re-condensation of vaporized oil fractions from the matrix.A cone-shape combustion front developed in the matrix core sample and swept it completely.The experiments were qualitatively matched by two-dimensional numerical simulations.The process was also explained by a simple analytical model.From this analytical formulation, the authors concluded that the ISC process is feasible in naturally fractured reservoir with maximum fracture spacing in the order of 1 m (3.28 ft), which means a very dense fracture network.
To date, air injection is not a widespread process, probably due to safety aspects, the cost of compression facilities and the complexity of involved mechanisms.Nevertheless, this process is worth being considered for improving the oil recovery of some fractured reservoirs where air injection can significantly enhance matrix-fracture transfers.Among them, the fractured chalky reservoirs of the North Sea, as found in the Ekofisk field, could present favorable perspectives of application.Jensen et al. (2000) showed that, among various EOR methods applicable to the waterflooded Ekofisk field, only the hydrocarbon Water Alternate Gas (WAG) and Air Injection methods show favorable perspectives of application and are worth being assessed further.Indeed, the incremental oil recovery values predicted from reservoir simulations were 6.5% OOIP for Air Injection and 5.6% OOIP for CO 2 WAG.As mentioned above, major safety problems may occur in the case of early breakthrough of unconsumed oxygen at production wells via the fracture network.Assessing such a risk requires an efficient and reliable simulation methodology.Okamoto and Bourbiaux (2005) discussed the main physical mechanisms of air injection in a light oil fractured reservoir and presented a detailed workflow for an efficient numerical modeling of this process.

Conclusion
Over the last decade, gas injection has been the dominant EOR method for crude-oil fractured reservoirs.When CO 2 source is available, CO 2 flooding remains the major EOR process for fractured reservoirs.Application of EOR processes, other than CO 2 and polymer flooding, has been very limited in fractured reservoirs.This can change in the future.Alkali-Surfactant-Polymer (ASP) is a proven technology that may become a choice method for extending the life of mature fractured reservoirs.

STATE OF THE ART OF DUAL-POROSITY SIMULATORS
The dual medium concept was introduced by Barenblatt and Zheltov (1960).Applied to a fractured reservoir, this concept considers that the fracture medium on the one hand, the matrix medium on the other hand, behave like two flow continua interacting together.This concept was applied by Warren and Root (1963) to interpret well tests in the simplified framework of a fractured reservoir made up of a single flow continuum, the fracture medium, that locally interacts with porous parallelepiped matrix blocks only acting as a fluid source.The first attempt to apply the dual-porosity approach to simulate three-dimensional multiphase flows in a fractured reservoir was presented by Reiss et al. (1973).The matrix-fracture exchanges were represented by time-dependent source or sink functions.These transfer functions were derived from laboratory experiments or/and from the numerical simulation of the oil recovery mechanism in a finely-gridded single-block model.Rossen (1977) used the same approach with a semi-implicit handling of the source terms.Both models assumed that the oil recovery from a matrix block is only due to capillary and gravity effects and is only function of the elapsed time after water (or gas) comes in contact with the matrix block via the surrounding fractures.
The effects of compressibility and of fracture viscous flows on matrix-fracture transfers were neglected.Actually, it is very difficult to derive transfer functions that are only timedependent and valid at each stage of the field exploitation.Kazemi et al. (1976) presented a two-dimensional twophase water-oil model in which the matrix-fracture flow rate is related to the potential difference between matrix and fracture.The model accounts for imbibition but no gravitational term is included.He introduced a shape factor definition that is consistent with the finite-difference formulation of matrixfracture exchange in a single-cell matrix block.
Various numerical models were developed afterwards to simulate three-dimension flows and especially to take into account the impact of gravity forces between fracture and matrix.Additional physical mechanisms were also studied, such as block-to block interactions.Several features were developed in order to extend the use of dual-porosity/dualpermeability simulators to more complex recovery processes, involving compositional, chemical, thermal, and geomechanical numerical simulations.All these aspects will be detailed in the following paragraphs.

Use of an Explicit Gravity Term
As explained in Section 3, conventional equations of matrixfracture exchanges assume 1-D, horizontal flow between gridblock centers of the matrix and fracture.Several methods were considered to include the effect of gravity in the flow terms, via an explicit gravity term, or through the use of pseudo relative permeability and capillary curves (see Sect. 2.2), or also by subgridding the matrix blocks (see Sect. 2.3).Gilman and Kazemi (1983) extended the two-phase Kazemi's model to three dimensions and improved the formulation of gravity forces between fracture and matrix by using two different depths for the fracture and matrix nodes.As pointed out by Sonier et al. (1986), this approach implies that the gravity forces applied on matrix blocks would remain constant although the gravity head on the matrix blocks changes with the fluid saturation in both media.Litvak (1985) and Sonier et al. (1986) improved this modeling approach.Their matrix-fracture exchange formula introduced a transient (time-dependent) gravity term calculated from the height of the matrix blocks and from the saturation within the matrix and surrounding fractures.Assuming a gravity-driven Vertical Equilibrium (VE) of fluids in both media, fracture and matrix saturations can indeed be converted into hypothetical fluid contact levels that are used to estimate the gravity head applied on the matrix blocks of the considered cell.Bossie-Codreanu et al. (1985) included the impact of gravitational forces on matrix-fracture transfers by using a specific cell arrangement.The basic feature is that an elementary volume of the fractured reservoir is simulated by several cells: the matrix is concentrated into one "matrix cell" surrounded in each direction by "fracture cells".The "fracture cells" offer a continuous path for fluid flow while "matrix cells" are discontinuous.This formulation, based on the dual porosity concept, allows the direct calculation of matrix-fracture flows, taking into account capillary and gravity effects.Quandalle and Sabathier (1989) also formulated the gravity contribution to matrix-fracture transfers through the estimation of a gravity head.The latter involves a given matrix block height per cell and time-dependent fluid densities in the fracture medium.In addition, the matrix-fracture transfer is expressed for each face of the parallelepiped blocks as a function of different terms representing the respective contributions of all involved physical mechanisms of transfer, including viscous transfers linked to fracture flows.

Use of Pseudo-Parameters
Thomas et al. (1983) modeled the effect of gravity on matrixfracture transfers by using pseudo-relative permeability and pseudo-capillary pressure curves based on the vertical equilibrium assumption.An additional term was also added in the matrix-fracture flow terms to account for the effect of the pressure gradient across the matrix block and for diffusion of gas into the undersaturated oil of matrix blocks during nonequilibrium gas injection.Dean and Lo (1988) demonstrated that the effect of gravity drainage could be included in pseudo-capillary-pressure terms for both the matrix and fracture without the need to include explicitly a gravity term.They determined these pseudo-capillary-pressure terms by history matching a reference fine-grid model of a single matrix block.Rossen and Shen (1989) used pseudo-capillary-pressure curves for both the matrix and fracture.The fracture curve can be determined directly from rock properties and matrixblock dimensions, while the matrix curve can be obtained from the results of a single simulation of a fine-grid model of a single matrix block.Coats (1989) used pseudo-functions as Rossen and Shen (1989) and discussed different formulations to simulate gravity effects.
The drawback of these approaches is to generate dynamic pseudocurves in an automatic and general way.As discussed in Bourbiaux (2010) the determination of effective multiphase permeabilities, or of pseudo-relative permeabilities, has no general solution.Except for situations of quasi-instantaneous matrix-fracture transfer times, the solution for pseudo-relative permeabilities is far from unique as it closely depends on the flow mechanism under consideration, the fracture flow conditions in terms of rate, or saturation, composition, temperature, and on the field flow history.Different conventional methods exist to compute those parameters, such as for instance the numerical resolution of steady-state flows on fine-grid models with various fractional flows of the fluids in presence.However, none of them has a general scope of application.

Subgridding Matrix Blocks
The necessity of a detailed simulation of the matrix-fracture transfers using gridded matrix blocks was underlined by Saidi (1983), who implemented a model where the matrixfracture transfers in each reservoir zone characterized by given rock properties are simulated using a gridded matrix block model surrounded by gridded fractures.The first twenty-two years of the Haft-Kel Iranian field were successfully matched using this approach.
A subgridded two-phase and heat flow model was proposed by Pruess et al. (1985) with the Multiple Interacting Continua (MINC) approach.Chen et al. (1987) developed a three-dimensional compositional thermal simulator allowing the rock matrix block to be subdivided into a two-dimensional (r, z) grid in order to study more precisely the effects of gravity, capillary pressure, and mass and energy transfer between fractures and matrix blocks.They concluded that oil recovery predictions in steamflooding are sensitive to capillary pressure values, the number of cells used for subgridding, and the matrix block size.Gilman and Kazemi (1988) proposed a matrix-block grid refinement with a subdomain approach, which gives accurate results and assumes a vertical segregation in the fracture only.
The drawback of the method is the computational cost which limits its range of application to only small reservoir models.Vossoughi et al. (1997) presented an improved solution method of the MINC procedure which results in a reduction of computation time and memory requirement.Famy et al. (2005) proposed an optimal sub-gridded matrix block model for the capillary imbibition case, by taking into account the physical specificities of this mechanism while minimizing the computational cost.They developed an optimized one-dimensional sub-gridded matrix block model on the basis of the local saturation evolution within the block.The methodology has been validated by comparison with reference fine-grid simulations, for various block shapes and rock-fluid properties, and for anisotropic flow conditions.The reference solutions are reproduced very accurately.Naimi-Tajdar et al. (2006) implemented the MINC approach in a parallel, 3-D, fully implicit, equation-of-state compositional model.The matrix blocks are discretized into subgrids in both horizontal and vertical directions.This subgrid scheme was used to model three-dimensional flow in the matrix rock with two-dimensional subgrids with good accuracy in the cases of water and gas injections.

Taking into Account Matrix Blocks Interaction with Dual-Porosity/Dual-Permeability Formulation
The preceding models deal with the representation of systems in which the matrix blocks are smaller than the gridblocks and for which fluid transfer in the matrix between gridblocks can be neglected.If it is not the case, such transfer must be considered in the flow equation with a dual-porosity/dual-permeability formulation.Because the matrix is connected between gridblocks, partial elimination, as described by Thomas et al. (1983), cannot be used to reduce the number of unknowns of the linear system to be solved.Blaskovich et al. (1983), Hill and Thomas (1985), Quandalle and Sabathier (1989) presented the physical and numerical aspects of the dual-porosity/dual-permeability approach.This formulation is useful to take into account "block-to-block" matrix flows due to partial matrix continuity in some parts of the reservoir.Nowadays, this formulation is available in most fractured reservoir simulators.Gilman and Kazemi (1988) used a dual-permeability approach in the vertical direction alone to accurately simulate the matrix flows involved in gravity-driven matrix-fracture transfers, thanks to a refined vertical discretization.Por et al. (1989) presented a dualporosity/dual-permeability with additional connections between the matrix node of a given cell and the fracture nodes of neighboring cells for modeling the block-to-block interactions and the role of capillary contacts during gravity drainage processes.

Other Dual-Porosity Simulation Features
Gilman and Kazemi (1983) added polymer and tracer equations in fracture and matrix media.Compositional dualporosity simulators started being developed in the late 80's with Quandalle and Sabathier (1989), Coats (1989) and Peng et al. (1990).Firoozabadi and Thomas (1990) presented the results of the Sixth SPE Comparative Solution Project on dual-porosity simulators.Two problems, gas-oil gravity drainage in a single matrix block, and depletion, gas or water injection in a cross-section model, were submitted to ten participating organizations owning most of the industrial dual-porosity model versions available at that time.Results differed considerably from one simulation model to the other, because of different formulations used to simulate matrix-fracture exchanges.The authors concluded by the need of further development of the physics and numerical modeling of fractured reservoirs.
More recently, Winterfeld (1996) presented a dual-permeability thermal simulator to simulate steam injection in the Yates field.Sabathier et al. (1998) presented two major improvements of their simulator to predict the matrix-fracture transfers by capillary imbibition and gravity drainage.
Recent developments in fractured reservoir simulation concern the possible impact of geomechanical parameters on the fracture flow properties, porosity and permeability (Chen and Teufel, 2000;Bagheri and Settari, 2008), the need of high-performance simulators with the use of parallelization techniques for giant fractured field simulation (Al-Shaalan et al., 2003) and the use of streamlines for waterflood simulations with three forms of the transfer function (Di Donato et al., 2003), and also the implementation of PEBI grids in order to keep a better track of the actual/geological fault/fracture network than conventional dual-porosity models (von Pattay and Ganzer, 2001).
The results of experimental studies of matrix-fracture transfer mechanisms -such as spontaneous imbibition (Bourbiaux and Kalaydjian, 1990) and diffusion (Le Gallo et al., 1997) -also enabled to formulate and validate physically-sound formula of these transfers in reservoir simulators.

Mathematical Formulation
The simulation of naturally fractured reservoirs is based on the dual-medium concept initially introduced by Barenblatt and Zheltov (1960).In this approach, the fractured reservoir is assumed to be consisted of two continua, the fracture and the matrix, that exchange fluids together, which is a dualporosity dual-permeability model by definition.Warren and Root (1963) proposed a simplified geometrical representation of the reservoir in order to facilitate the formulation of matrix-fracture transfers.The latter consists in an array of identical matrix blocks delimited by an orthogonal set of equidistant fractures oriented along the main directions of flow within the reservoir.Furthermore, Warren and Root did not consider any communication between matrix blocks, that is, a single-permeability version of the dual-medium modeling approach.As shown in Figure 2, the numerical model of a dualmedium reservoir consists in two identical superposed grids, a fracture grid and a matrix grid.Fluid flows are simulated between neighboring cells of the fracture grid as well as a transfer flux between the two superposed fracture and matrix cells at any grid location.A convenient way to express that matrix-fracture flux for a given couple of matrix and fracture cells is to consider the fractured medium as identical matrix parallelepipeds delimited by orthogonal sets of regularlyspaced fractures.Actually, matrix-fracture flow transmissivities can then easily be formulated from the dimensions of this "equivalent" matrix block.Figure 2 is obviously a coarse representation of the actual fractured medium that is constituted of blocks of various sizes and shapes.This representation introduces further inaccuracy into the matrix-fracture exchange prediction made by dual-porosity simulators, in addition to the inaccuracy resulting from the upscaled expression of transfers for non-discretized matrix blocks.Although various strategies have long been studied to improve matrix-fracture exchange predictions from dual-porosity modelsincluding transient formulas, use of two equivalent blocks, or sub-gridding -most industrial simulators remain based on the initial representation of fractured reservoirs as proposed by Warren and Root.
The equations governing three-phase, multi-component, three-dimensional flow in naturally fractured reservoirs, are written for the fracture and the matrix systems as follows.

Flow Equations for Fracture System
For a unit volume of reservoir, the mass balance equation for each component k, including water, is expressed as: (1) With the exception of the matrix-fracture mass transfer term F mf kp , the fracture system equations are the same as those used in single-porosity simulators.

Flow Equations for Matrix Blocks
In the case of a dual-porosity single-permeability model, the mass balance equation for each component k, including water, is expressed as: (2) Dual porosity representation of a fractured reservoir.
In the case of a dual-porosity dual-permeability model, the mass balance equation for each component k, including water, is expressed as: (3) In the previous Equations (1-3), superscript f refers to the fracture, superscript m to the matrix and subscript p to the phase.φ (φ f ) are the matrix (fracture) porosity defined as the pore volume of medium m(f) divided by unit volume of both media.

C M
kp , S M p are (for medium M = f, m) the mass fraction of component k in phase p and the saturation of phase p.

F mf
kp is the matrix-fracture mass flow rate of component k in phase p per unit bulk volume of reservoir.The formulation of this term will be described in the companion paper (Part 2).

Phase velocity
The phase velocity is expressed in both media as usually with the Darcy equation: (4) where Z is depth (positive, increasing downwards), g is the algebraic value of gravitational acceleration projection on Z axis.K is the absolute permeability tensor of the medium, p p the pressure of phase p, μ p the viscosity of phase p, k rp the relative permeability of phase p. Full tensor permeabilities are often required to describe complex reservoirs, especially fractured reservoirs.Use of positive definite symmetric full tensors requires more complex discretization schemes.That is why most simulators assume the permeability can be represented by a diagonal tensor whose principal axes coincide with the local coordinate axes.
Note that in Equation (4), K is the "equivalent" permeability of the flow continuum under consideration, that is the permeability of the medium that is scaled up to the simulated flow scale.The permeability of the fracture large-scale continuum is obtained in the Warren and Root model by homogenizing the transport equations on the fracture network, assuming impervious boundaries.Using the homogenization theory (Bourgeat et al., 1999) or large scale volume averaging techniques (Quintard and Whitaker, 1996), the upscaled flow properties involved in the set of Equations (1-4) were derived for single-phase flow conditions.

Additional equations
Twelve additional equations are added to the flow equations, six for the fracture and six for the matrix, expressing the sum of saturations: the sum of compositions in each phase, water (w), oil (o) and gas (g): (6) and capillary pressure relationships:

Equilibrium equations
If N c is the number of components k, 4N c equilibrium equations complete the set of equations, 2N c for the fracture and 2N c for the matrix: where K kgo and K kgw are respectively the gas-oil and gaswater equilibrium constants for the kth component.

Diffusion-dispersion flux
The diffusion-dispersion flux of a component k in a phase p in Equations (1-3) is expressed for both media, fracture and matrix, as follows: (11) for M = f and m.
with the diffusion tensor given by: β is the dispersivity of the medium.τ x , τ y and τ z are the pore tortuosities in x, y and z directions.The diffusion tensor is assumed to be diagonal with principal axes coinciding with the local coordinate axes.The same remark as for the permeability tensor in Equation ( 4) is valid.The diffusion-dispersion term is written in the same way for both media, but the molecular diffusion is essentially prevailing in the matrix medium whereas the dispersion is essentially prevailing in the fracture medium.
In a fractured medium, diffusion phenomena essentially concern the matrix medium where phase velocities are very slow, whereas dispersion concerns the fracture medium.The simulation of dispersion can hardly be predictive because the dispersion flux of a given component in a given phase depends on both the dispersion coefficient and the relative permeability function via the phase velocity term.Scaling of these parameters to cell size remains an open question for natural fracture networks.The phase velocity is increased if linear relative permeability curves are used and then the dispersion flux may be overestimated if these linear relative permeability curves are non valid.For instance, Hamon et al. (1991) used non-linear fracture relative permeability curves for matching the Meillon fractured gas field, as linear fracture relative permeabilities frequently resulted in too early water breakthrough in the simulations.
For each component k, the set of Equations ( 1) to ( 12) represents (6N c + 12) equations involving (6N c + 12) unknowns, that comprise, for each medium, the three phase pressures p w , p o , p g , the three phase saturations S w , S o , S g , and the 3N c mass fractions, C kw , C ko , C kg , of the N c components denoted k.
The equations and the representation of Warren and Root previously described have been presented for an oil fractured reservoir but are also valid for a gas fractured reservoir.The main principles are the same.Some particular mechanisms may be taken into account, such as non-Darcy flow effects in the fractures near the wells.The production of unconventional tight gas reservoirs involves specific flow and transfer mechanisms, that are the Klinkenberg effect and gas desorption from pore walls on the shale matrix, which were taken into account by Kucuk and Sawyer (1980) to simulate well tests in Devonian shale gas reservoirs.

Alternative Approaches in Fractured Reservoir Simulation
In some particular cases, a fractured reservoir study may be carried out with a single porosity model or a discrete fracture network model instead of a dual-porosity numerical model.These alternatives approaches are presented in this section.
Single-porosity and dual-porosity models are based on homogenized representations of the fractured reservoir either as a single equivalent medium or as two equivalent fracture and matrix media.
If the matrix medium does not contain or produce any oil at all or if a very dense and well-connected fracture network ensures fast matrix-fracture transfers compared to the largescale transport through the fracture network, then the singleporosity model can be reliably used, after a homogenization procedure.That situation of quasi-instantaneous matrix-fracture transfers is often encountered during the single-phase depletion phase of many fractured reservoir production histories.Another application example of the single-porosity approach can be found in Van Lingen et al. (2001).The single-porosity formulation was successfully applied to a Middle East carbonate reservoir, produced by peripheral water injection and with significant primary matrix productivity, and a large spacing between large-scale disconnected fractures (in the order of one or several hundred meters).The key element in the homogenization procedure for obtaining the rapid water advance due to the presence of conductive fractures was the use of pseudo-relative permeability curves for grid blocks containing fractures.These curves were determined through an analytical procedure, based on the local fracture and matrix properties and based on the assumption that the fracture volume of a grid block is filled with water prior to the imbibition of water into the matrix.
On the contrary, for typical dual-medium reservoirs with a high flow-property contrast between a connected fracture network and a porous matrix medium, the contributions of fractures and matrix to the reservoir production cannot generally be simulated as the contribution of a single equivalent medium because the low-permeability matrix medium flow response is delayed to a variable extent with respect to that of fractures.A dual-porosity model is then required to simulate separately fracture flow and matrix flow.For such fractured reservoirs, the notion of "matrix block" is introduced, which at any reservoir location represents a fluid "source" (or "sink") term that feeds the fracture medium.Dual-porosity models are extensively used in the petroleum industry to simulate the various production schemes of this last kind of fractured reservoir.Dual-porosity dual-permeability models are adopted for the most general situations where the matrix medium also contributes to a certain extent to large-scale flow.
Some authors recently implemented different approaches to simulate flow directly on the discrete fractured medium generated by fracture characterization software without any prior homogenization.They used discrete element models of the fracture network alone or of both the fracture network and the matrix blocks.Sarda et al. (2002) used dual-medium flow equations to simulate slightly-compressible flow on a Discrete Fracture Network (DFN) model with an explicit treatment of matrix flows and matrix-fracture transfers.Basquet et al. (2004) developed a version of that simulator dedicated to highly-compressible gas flow.That version takes into account the high fluid compressibility via a pseudo-pressure, and the non-Darcy flow near the wellbore via a rate-dependent skin.Granet et al. (1998Granet et al. ( , 2001) ) developed a 2-D discrete model of fractures and matrix to simulate single-phase and two-phase flow on the realistic geological images of fractured reservoirs generated by a stochastic fracture modeling software.They used an unstructured grid made of two kinds of elements, linear elements for the fissure network and triangular elements for the matrix with an appropriate finite volume scheme.Bourbiaux et al. (1999) applied this model to the simulation at local scale of the water imbibition of a fractured medium, involving two fracture sets generated by a stochastic fracture modeling software as shown in Figure 3.The fractured medium horizontal dimensions are equal to 200 m (656 ft).The size of the zoom is 20 × 16 m 2 (65.6 × 52.5 ft 2 ).
As illustrated in Figure 4 for water/oil imbibition, this approach can provide reference solutions to validate or calibrate the matrix-fracture transfer formulations used in dualporosity simulators.
Figure 4 Matrix oil recovery for the geological image (water/oil imbibition) (after Bourbiaux et al., 1999).Kim and Deo (2000) used a Galerkin finite element formulation for the same kind of simulations on 2-D discretefracture models.
These approaches are very interesting but remained limited to 2-D networks for a long time.Actually, mesh generation for 3-D networks was difficult and numerical simulation very CPU-time consuming.However

Gravity-Induced Phenomena
The internal complexity of fracture space lets one predict a widely-spread distribution of fracture conductivity values, henceforth, a heterogeneous distribution of flows among the fractures of a natural network.At field scale, single-phase flows in the fracture network are analyzed and modeled using Darcy law, considering the fracture network behaves like a continuous "fracture" medium.The permeability of this medium is the equivalent permeability of the fracture network considered as a whole at the scale of the drainage area of a given well.However, attempts are made in practice to have access to the "conductive fracture network" in order to predict flow anisotropy and preferential flow paths through the reservoir.Interference tests and tracer tests are the two main types of multiple-well tests enabling such a characterization.The main transfers taking place in a vertically-connected fracture network are segregation and natural convection.The geological image: overall view and zoom (with the matrix discretization with triangles) (after Bourbiaux et al., 1999).

P Lemonnier and B Bourbiaux / Simulation of Naturally Fractured Reservoirs. State of the Art Part 1 -Physical Mechanisms and Simulator Formulation
Provided that the structure of the reservoir is favorable, i.e. characterized by a large thickness and/or a high dip, and assuming a well-connected, continuous and conductive fracture network, the solution gas liberated from the fracture and matrix oil during the reservoir depletion can segregate within the fracture network and accumulate in the original gas-cap or form a secondary gas-cap.Thus, the field energy, associated with the gas expansion potential, is saved contrary to more conventional situations where gas and oil are driven to the wellbore according to their respective mobilities.
Natural convection is a complex phenomenon that may take place in thick or high-relief and highly-fractured reservoirs, like Iranian fields (Saidi, 1996) or Mexican Cantarell field (Manceau et al., 2001).As explained by Peaceman (1976a) and Saidi (1987), the driving force for the convection process is an inversed oil density gradient with depth, with heavier fluid on top of light fluid.This involves a gravity segregation process with gas liberation at the gas-oil contact, provided that a dense well-connected fracture network makes possible the vertical circulation of fluids.
More precisely, the existence of an inversed density gradient within the oil column is expressed as: a condition that may be expressed versus thermal and pressure gradients as: with c o the isothermal compressibility of oil defined as: and β o the isobaric thermal expansion of oil defined as: In practice, in addition to the fluid PVT behavior, one needs to dispose of reliable estimations of pressure and temperature gradients within the oil column in order to estimate the above convection criterion.Reiss (1980) quotes an Iranian field example where β o = 1.15 × 10 -3 vol/vol/°C -1 , c o = 1.02 × 10 -4 vol/vol/bar -1 , dT/dz = 3.64 × 10 -2 °C/m and dp/dz = 0.069 bar/m, leading to the following negative normalized density gradient versus depth: (13) One has to be aware that this figure represents an oil density increase of only 0.35% of its value over one hundred meters.Considering usual uncertainties on the terms involved in such a calculation, the question of convection occurrence remains often unsolved in practice.
Other factors favorable to the initiation of convection are a large reservoir thickness, high vertical fracture permeability, and also a low oil viscosity, to allow the onset of significant density-driven flows.Actually, convection flows in the fracture network generate a compositional gradient between the matrix medium and the fractures at a given depth.This gradient is responsible for matrix-to-fracture molecular diffusion phenomena that tend to annihilate convection processes if too weak or slow.
A convection process can remain active during field depletion because the heavier oil obtained at the gas-oil contact by solution gas liberation can migrate downwards through the oil column, while the gas-saturated oil of deeper matrix blocks diffuses to the fractures and segregates upwards.These convection-segregation flows that take place in a highly-permeable fracture network tend to homogenize the oil properties along the vertical depth, and lead to further liberation and segregation of gas to the gas-cap, thus explaining the observed decrease in production GOR from undersaturated reservoir zones (as shown in Fig. 5).Therefore, convection diagnosis is based on the following production observations, among others: -homogeneous or low gradients of oil-property and temperature vertical profiles; -the evolution of the produced oil composition with time for wells completed in the oil zone: decline of the bubble point pressure and of the production GOR with time.
There are very few publications about modeling of convection in fractured reservoirs (Peaceman, 1976a, b;Saidi, 1987;Ghorayeb and Firoozabadi, 2000;Manceau et al., 2001).All these publications are related to numerical simulation in two-dimensional (x, z) phenomenological models, Figure 5 Typical trend of production GOR of a fractured reservoir subjected to convection and segregation phenomena (after Manceau et al., 2001).
except for the last one of Manceau et al. (2001), that deals with a field production history.Peaceman (1976a) modeled the convection phenomenon caused by a fluid density inversion using a perturbation analysis approach, and identified some of the significant variables that affect the process of convection within the fissures with and without matrixfissure transfer.The mathematical theory was corroborated by numerical calculations (Peaceman, 1976b).Peaceman established some formulae in a two-dimensional vertical gasoil system for the rate of finger growth and the size of the convection cell, which depend on the magnitude of the density inversion, the oil viscosity, the matrix-fissure transfer rate and the ratio of matrix to fissure pore volumes.He concluded that convection cells, if formed, are as high as possible.Saidi (1987) made some comments on the derivation of the equations used by Peaceman and highlighted some problems.Saidi presented analytical solutions of convectiondiffusion in two-dimensional radial-circular and cross-section domains.He gives the velocity distribution in convective cells as well as the width of the convective cells for different boundary conditions.The effects of fracture parametersincluding fracture aperture or conductivity, fracture network density and connectivity -on the vertical distribution of fluid composition were studied by Ghorayeb and Firoozabadi (2000) in a rectangular two-dimensional (x, z) fractured porous medium model containing a two-component singlephase fluid under a prescribed linear temperature gradient with depth.Several simple configurations of fractured porous media were considered.Numerical results confirm that a pronounced convective motion takes place in the most open, i.e. most conductive, fractures.Ghorayeb and Firoozabadi observed from the simulation of 2-D fractured medium configurations that the main flow loops within the fractured medium follow major surrounding fracture, whereas other nested loops showed slow fluid movement.All these publications present the effect of convection on the behavior of fractured porous media but do not give any ways of including this process in 3-D numerical simulators.Whereas previous studies were 2-D and phenomenological, Manceau et al. (2001) simulated a fractured field production history taking into account the assumed occurrence of convection within the oil column.The convection phenomenon was mimicked by a strongly-amplified vertical diffusion process within the oil phase saturating the fractures.That is, a dual-porosity simulator was used with a very large vertical diffusion coefficient in the fracture grid.Increasing vertical diffusion results in a rapid change of oil properties over the whole thickness of the reservoir.As gas is liberated at the gas-oil contact, the changes of the oil composition are transferred to the oil directly below, which reproduces exactly the effect of convection.Therefore, convection was taken into account by implementing in the reservoir simulator a multiplying factor for the vertical oil diffusion coefficient in the fractures.
Moreover, because convection involves liberation of the light oil components at the gas-oil contact, the Pressure/ Volume/Temperature (PVT) representation of the fluids has to be as precise as possible and the use of a compositional model is strongly advised.
The previous approach with a multiplying factor of 10 6 was successfully used by Manceau et al. (2001) to simulate the production history of the major Akal reservoir located in the giant Cantarell oil complex in Mexico.The convection cells or "rolls" were not explicitly modeled as described in the papers above mentioned of Peaceman, Saidi, Ghorayeb and Firoozabadi, due to a lack of precise description of the fracture network properties.Only the results of the natural convection process were simulated and history matched.For Cantarell field case, convection was simulated as an accelerated diffusion process between superposed fracture cells, which is equivalent to a convection process ending with fluids vertical equilibrium.The simulated process has to be tuned to match the field compositional production history.That is, the overall kinetics of the process was simulated without trying to reproduce the physical loops.
Sensitivity tests led to the following conclusions: -convection does not much influence the evolution of the water-oil contact but contributes to the development of a larger gas cap in conjunction with segregation; -convection has a direct and strong impact on the fluid composition distribution within the reservoir, in conjunction with segregation, matrix-fracture diffusion and thermodynamical transfers; -in consequence, convection has a strong influence on the evolution of produced fluid composition.It was clearly observed, both on real field data and on the model, that the wells located below the gas cap had a GOR that significantly declined with time, while, on the reservoir edges, the GOR remained constant.For field production purposes, one has to keep in mind that convection tends to homogenize the fluid composition of the oil zone and leads to further liberation and segregation of gas to the gas-cap, thus explaining the observed decrease in production GOR from undersaturated reservoir zones.

Relative Permeabilities for Multiphase Flow in Natural Fractures
Experimental quantification of multiphase flows in fractures is difficult and hardly any data are reported in the literature.
Whereas major faults may have to be modeled individually as no-flow boundaries or as preferential flow paths depending on their conductivity, all other fractures are most often represented as an equivalent medium at the resolution scale adopted for field assessment.Darcy generalized law then applies to this equivalent medium, involving capillary pressure and relative permeability functions.Capillary pressures are often neglected, or taken as a constant value.This constant value may be considered as the capillary pressure of the fracture represented as two walls with a constant spacing, a very simplistic assumption indeed.
As pointed out by Deghmoum et al. (2000), reliable fracture relative permeability data are required for predictive dual-porosity reservoir simulations.Chen et al. (2004) showed that two-phase flow experiments in a single fracture do not yield linear relative permeability curves.However, such measurements of relative permeability hardly reflect the reservoir flow conditions because of laboratory limitation and of non-representative core samples.Due to these limitations, fracture relative permeabilities are generally assumed to be straight lines.This assumption may lead to erroneous results as it corresponds to the absence of any mutual interaction between fluids when flowing, which is only consistent with negligible capillary forces.Hamon et al. (1991) observed that upscaling of fine grid waterflood simulations results in non-linear fracture relative permeability for the network even if the individual fracture relative permeability is assumed to be linear.They used non-linear fracture relative permeability curves for matching the Meillon fractured gas field, as linear fracture relative permeabilities frequently resulted in earlier water breakthrough and in smaller pressure drop than observed in the field.Rossen and Kumar (1992) developed a percolation model for two-phase flow in a single fracture that is based on a gravity segregation model defined as the limit of the gravitycapillary equilibrium model.The range of application of this gravity segregation model is quantified by the following dimensionless parameter: (14) where Δρ is the density difference between phases, g is the gravitational acceleration, H is the fracture height, γ is the interfacial tension, and ε 0 is the mean value of the fracture half-aperture.For large values of H D , relative permeability curves approach the straight-line form, assumed in many reservoir simulations.
Then, Rossen and Kumar (1994) conducted a simulation study for a two-phase flow model in a single fracture.They studied the effects of fracture relative permeabilities on the performance of the waterflood case proposed in the Sixth SPE comparative Simulation Project on naturally fractured reservoirs (Firoozabadi and Thomas, 1990).The simulation results showed that low fracture relative permeabilities can significantly reduce the oil production rate as soon as the flow rate target of production wells can no more be satisfied.For the SPE 6 waterflood example, this occurs for a reduction in total fracture fluid mobilities by a factor greater than about five.However, altering fracture relative permeabilities was shown to be roughly equivalent to modifying the fracture absolute permeability.Therefore, if the field-estimated value of the fracture network permeability is an effective permeability inferred from multiphase flow, then using incorrect relative permeabilities have a minor impact: an underestimation of the absolute permeability would roughly balance an overestimation of the relative permeabilities.However, if the fracture absolute permeability is determined from a singlephase flow or from measurements of hydraulic apertures in core samples, then the effects of assuming wrong relative permeabilities can be significant.The main conclusion is that it is therefore preferable to adopt effective permeabilities rather than trying to measure absolute and relative permeabilities separately.Rossen and Kumar (1994) also showed that fracture relative permeabilities also affect reservoir sweep efficiency by altering the ratio of viscous to gravity forces in the fracture network.
Fourar and Lenormand (1998) interpreted two-phase flow measurements in a single fracture using a simple model based on viscous coupling between fluid phases.The fracture is modeled by two parallel planes with a small aperture.Using Stokes' equation, they obtained the following analytical equations for the two-phase liquid-gas relative permeabilities: (15) (16) These equations show that the relative permeability of the non-wetting fluid depends on the viscosity ratio R μ .
For the case of high flow rates in fractures, Fourar and Lenormand (2001) developed a k r model for two-phase flow with inertial effects, based on a generalization of the singlephase Forchheimer's equation by introducing a function F as a multiplying factor of the velocity.The function F can be derived from Equations (15-16), established in a previous paper (1998).The formulation was validated using experimental results with air and water.

Simulation of Conductive Faults
There exist several methods for modeling conductive faults in a reservoir.An example of such a reservoir, Fateh Mishrif in the Arabian Gulf, is shown in Figure 6 which presents fluid-convective fault sections indicated by dashed faultplane lines.The most straightforward method is to use very narrow gridblocks to model both the fault and the vicinity of the fault.This method leads to a good physical representation of the fault zone but unfortunately to unacceptable computation costs (Hearn et al., 1997).
A more flexible approach consists in modeling the faults as part of the fissure grid system of a dual-medium model (Gilman et al., 1995).However, this method cannot dissociate the roles of fissures and faults in multiscale, fractured and faulted reservoirs.
Conductive faults could also be modeled by normal-size gridblocks that include both the fault and the neighboring reservoir medium (Trocchio, 1990).The fault gridblocks are assigned equivalent permeabilities to simulate the fluid flow in the faults.This last method entails three drawbacks: -building such a model is time-consuming for a complex network; -it overestimates breakthrough times; -coarse grids do not allow simulating development scenarios with a close well spacing in densely-faulted zones.Lee et al. (1999) modeled conductive faults similarly to wells in a conventional reservoir simulator.They introduced a transport index, similar to the productivity index in well modeling, to formulate the fluid exchanges between the fault and the gridblock enclosing the fault.As the wellbore flow equations imply no accumulation terms, this approach does not model the volume and the distribution of fluids within faults, and cannot predict the gravity exchange between the fault and the neighboring medium.Moreover, like the previous method, this approach is not adapted for the simulation of infill drilling scenarios in densely-faulted regions.
To overcome the drawbacks of the methods described above, Henn et al. (2004) developed a specific model of conductive faults that ensures a precise and flexible description of fault geometry and fluid displacements, while keeping a high numerical performance.This model is based on an explicit representation of conductive faults and on the assumption of fluid gravity segregation within the faults, with negligible effects of capillarity.The vertical distribution of fluids within faults then simply obeys the Vertical Equilibrium concept (Coats et al., 1967) reduced to the sole effect of gravity forces.Under these conditions, phases in the fault are vertically segregated according to their densities.Thus, conductive faults can be modeled with a unique gridblock in the vertical direction.The horizontal flow of segregated phases along fault planes is formulated with generalized Darcy's law using relative permeability and capillary pressure pseudo-functions.Calculation of these pseudo functions is detailed in Henn et al. (2004).
Although such a model solves pressure and saturation unknowns at the scale of entire fault columns, the location of contacts between segregated fluids within each fault cell is re-calculated from saturation values in order to correctly formulate the exchanges of phases with the matrix background: this method is denoted as "a local approach" of fault exchanges with the surrounding medium (Fig. 7).
Coupling a conductive fault model with a dual-medium model leads to a triple-medium model (Fig. 8).The question is then to assess whether the fault-matrix exchanges have a significant contribution, compared to the matrix-fissure and fault-fissure exchanges.
A phenomenological study was carried out to evaluate the influence of the connections between the fault and matrix media.Fluid flows in the fissured medium, and fault-matrix and fault-fissure exchanges, depend on matrix block sizes and on fissure and matrix permeabilities.The fissure-to-matrix Figure 7 Local approach (after Henn et al., 2004).Triple-medium model (after Henn et al., 2004).
equivalent permeability ratio K * fm , as defined by Quintard and Whitaker (1996), is used with the assumption that e f << e m : (17) K i is the local permeability of medium i, i = m, f where m and f correspond to the matrix and fissure media respectively, e m is the block dimension and e f is the fissure aperture.
The numerical simulations performed by comparison with a reference solution showed a small influence on the oil recovery from all the matrix blocks.
Figure 9 shows that, whatever the value of K * fm , the oil recovery is not delayed very much when fluid exchange is annihilated between the fault and the matrix media.Moreover, that delay decreases when K * fm increases.That is, for most gas-oil flows, fault-matrix exchanges are negligible and the triple medium schematized in Figure 8 is reliable.
The triple-medium simulation approach was also validated for a gas-oil displacement in a three-dimensional sector of a field crossed by a dense conductive fault network (see Fig. 6 above).The injector and the producer are located within the reservoir fault panel far from the fault plane.Wells are defined in the fissure grid.At the initial time, the matrix and fissure media are saturated with oil.As the main fault-fissure exchange process is gravity drainage, production results are independent of the horizontal size of gridblocks edging the fault plane.Therefore, the dual-porosity model is coarsely gridded in the horizontal direction X.Two models, where the fault is explicitly gridded in the vertical direction, were built and tested as reference solution models: -model 1: fault is gridded with 16 gridblocks in the vertical direction; -model 2: fault is gridded with 4 gridblocks in the vertical direction.
As the fault is finely gridded, model 1 gives an accurate prediction of gravity drainage and the correct two-phase flow response at the production well.The evolution of the production gas-oil ratio at producer versus time is plotted in Figure 10 for the three models: reference models 1 and 2 and conductive fault model.
The comparison of models 1 and 2 shows the influence of the vertical gridding of the fault.A coarsely-gridded model overestimates the production gas-oil ratio.The conductive fault model and model 1 give identical results.Concerning gas-oil fault-to-fissure exchanges, the main process is gravity drainage.To model it properly, an accurate calculation of the pressure and saturation values within the fault gridblocks is required.Thanks to the local approach, the conductive fault model allows estimating these values accurately and predicts Oil recovery from the matrix medium.Water-oil simulation (after Henn et al., 2004).
an accurate physical response in agreement with the finelygridded model (model 1).In addition, the numerical performances of that new model are good compared to those of the time-expensive model 1.By using the Vertical Equilibrium concept in conductive faults, the total number of gridblocks in a full-field model crossed by a dense conductive fault network is significantly reduced without any detrimental effects on model predictability.In particular, the multiphase flow behavior of producers is quite well predicted.The numerical stability provided by a fully-implicit numerical scheme also contributes to model performance in view of full-field long-term simulations.This fault modeling approach is designed for both single and dual-porosity simulations.The single medium could represent either a matrix medium or a homogenized fissured medium.Applications concern any reservoir study where conductive faults are suspected to influence field response.

Wells Modeling
A well is defined in a dual-porosity model as in a single-porosity one.The only particularity is the definition of the well perforation location.A well may be perforated either in the matrix medium, or the fracture medium, or also in both media.The well is represented by a source (sink) term in each cell crossed by the perforated well section.This source term is formulated through the use of a well Productivity Index (PI) that relates the well section flow rate to the difference between the average well cell pressure and the wellbore pressure, assuming a steady-state radial flow.The PI is generally calculated by using Peaceman formula (Peaceman, 1983) for vertical or horizontal wells, or by using an extension of this formula for inclined wells (Mochizuki, 1995).In the case of a fractured reservoir, the geometry of flows around wellbore is strongly conditioned by the geometry and properties of the fracture network in the near-wellbore region.Hence, the PI values based on Peaceman formula with the cell upscaled fracture permeabilities are often unrepresentative since the upscaling procedure is generally based on flow simulation in a parallelepiped model with lineartype pressure boundary conditions, which do not capture the near-well flow behavior.As a result, the well productivity calculated by a dual-porosity flow simulator can be very different from that calculated on a near-wellbore Discrete Fracture Network (DFN) model.More representative PI values can be derived from the simulation of near-wellbore flows on a discrete element model that includes the actual representation of the geological near-wellbore fracture network and adopts the proper flow boundary conditions.This PI determination procedure was shown to very much improve the simulation of wells in fractured reservoir simulation models (Ding et al., 2006).This improvement was demonstrated through the following synthetic example.
Two sets of fractures at different scales were generated in a reservoir of 7600 × 6400 m 2 and 10 m thickness.The first set consists in major sub-seismic faults at the reservoir scale (Fig. 11), whereas the second one is a set of small-scale systematic joints (Fig. 12).The reservoir is represented as a dualporosity model with a horizontal gridblock size of 100 m.
A horizontal well, 1000 m in length, is intersected by two faults as shown in Figure 11.Fault conductivities are 2500 mD.m (2.47 × 10 -12 m 3 ) and the average conductivity in diffuse fractures is 250 mD.m (2.47 × 10 -13 m 3 ).
A single-phase flow is simulated during 5 days with a flow rate of 500 m 3 /day.Simulation results are compared in Figure 13.The dual-porosity model with near-well upscaling gives results much closer to those of the reference DFN model than the dual-porosity model with conventional upscaling.At the end of simulation, the dual-porosity model error on the wellbore pressure is 0.9% for the PI based on near-well upscaling and 49.3% for the conventional PI.Wellbore pressure (after Ding et al., 2006).

CONCLUSIONS
A significant amount of the world hydrocarbons reserves is located in naturally fractured reservoirs.These reservoirs exhibit low matrix permeability and high fracture permeability and often oil-wet to mixed-wet rock properties, which results in low oil recovery.A good understanding of the flow behavior is required for identifying and optimizing the best recovery process.As highlighted in this paper, the multiple physical mechanisms involved in oil and gas recovery processes of fractured reservoirs make their simulation a complex and challenging issue.As illustrated by the state-ofthe-art section, literature on the subject is very rich and significant advances have been made since the sixties for improving the reliability and efficiency of fractured reservoir simulators, especially in the modeling of matrix-fracture transfers.Accordingly, dual-porosity simulators are widely used in the petroleum industry for field-scale simulations of recovery processes in fractured reservoirs.Waterflood and gasflood are the main processes of secondary recovery of fractured oil reservoirs, but there are opportunities for EOR methods as thermal and chemical methods.
Significant quantities of natural gas are produced from low-permeability carbonate and shale fractured reservoirs.The depletion of tight gas naturally-fractured reservoirs often gives rise to a water influx and consequently simulation requires a good understanding of the water-gas imbibition mechanism taking place in such fields.
Most naturally fractured reservoirs contain multiscale fractures/faults.In this context, modeling of conductive faults encounters physical and numerical difficulties.Several methods were presented in this paper with a focus on one method based on the vertical equilibrium concept.
Well modeling is in most of the cases based on Peaceman formula for Productivity Index (PI) calculation.We discussed in the last section the use of a more representative method derived from Discrete Fracture Network (DFN) simulation models of the near-wellbore region.
The matrix-fracture transfer formulations, as techniques for uncertainty analysis and geology&flow-constrained history matching, will be addressed in the companion paper, Part 2: Matrix-Fracture Transfers and Typical Features of Numerical Studies.

Figure 1
Figure 1Main drive mechanisms in fractured reservoirs.

P
Lemonnier and B Bourbiaux / Simulation of Naturally Fractured Reservoirs.State of the Art Part 1 -Physical Mechanisms and Simulator Formulation

P
Lemonnier and B Bourbiaux / Simulation of Naturally Fractured Reservoirs.State of the Art Part 1 -Physical Mechanisms and Simulator Formulation Figure 2 ρ Mp is the density of phase p in medium M = f, m, i.e. a function of pressure and composition.isthe velocity of phase p in medium M = f, m.is the molecular diffusion and dispersion flux of component k in phase p in medium M = f, m.Q Mp is the volumetric injection/production rate of phase p in medium M = f, m per unit volume of reservoir.The rate is positive in production and negative in injection.

P
Lemonnier and B Bourbiaux / Simulation of Naturally Fractured Reservoirs.State of the Art Part 1 -Physical Mechanisms and Simulator Formulation D kp is the molecular diffusivity of component k in phase p.
, the situation seems to change since Geiger et al. (2009a) recently developed a 3-D three-phase three-component Discrete Fracture and Matrix (DFM) model.Their DFM model is constituted by unstructured hybrid element meshes associated to a finite element -finite volume transport scheme.Geiger et al. (2009b) implemented these techniques in a massivelyparallel simulator.They presented high-resolution DFM simulations with thousands of fractures and up to 5 million elements.

Time
Figure 3

Figure 12 Sub
Figure 12Sub-seismic faults and systematic joints (the latter are shown over the bottom part of the figure only) (afterDing et al., 2006).

P
Lemonnier and B Bourbiaux / Simulation of Naturally Fractured Reservoirs.State of the Art Part 1 -Physical Mechanisms and Simulator Formulation

Figure 13
Figure13 NOMENCLATURE C kp Concentration of component k in phase p (kg/kg) D Molecular diffusivity (m 2 /s) F kp Matrix flow rate per volume of component k in phase p P Lemonnier and B Bourbiaux / Simulation of Naturally Fractured Reservoirs.State of the Art Part 1 -Physical Mechanisms and Simulator Formulation