Simulation of Naturally Fractured Reservoirs. State of the Art

Résumé — Simulation des réservoirs naturellement fracturés. État de l’art : Partie 2 – Échanges matrice-fracture et spécificités des études numériques — Les réservoirs naturellement fracturés contiennent une partie significative des réserves en huile mondiales. La production de ce type de réservoirs constitue un défi pour les ingénieurs de réservoir. L’utilisation des simulateurs de réservoir peut aider 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 premières publications sur le concept double-milieu dans les années soixante. Cet article et le précédent (Partie 1) présentent les techniques actuelles de modélisation utilisées dans les simulateurs industriels. Il n’y a pas de réponse définitive pour simuler de manière optimale les échanges matrice-fracture, et, différentes méthodes sont mises en œuvre dans les simulateurs industriels de réservoir. Ce papier se concentre sur la modélisation de la physique des écoulements au sein des milieux matrice et fracture et des échanges entre ces deux milieux afin de mieux comprendre les différentes formulations proposées dans la littérature. Plusieurs problèmes particuliers liés à la simulation numérique Abstract — Simulation of Naturally Fractured Reservoirs. State of the Art: Part 2 – Matrix-Fracture Transfers and Typical Features of Numerical Studies — Naturally fractured reservoirs contain a significant amount of the world oil reserves. The production of this type of reservoirs constitutes a challenge for reservoir engineers. Use of 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 and the preceding one (Part 1) present the current techniques of modeling used in industrial simulators. The optimal way to predict matrix-fracture transfers at the simulator cell scale has no definite answer and various methods are implemented in industrial simulators. This paper focuses on the modeling of physical mechanisms driving flows and interactions/ exchanges within and between fracture and matrix media for a better understanding of proposed flow formula and simulation methods. Typical features of fractured reservoir numerical simulations are also described with an overview of the implementation of geomechanics effects, an application of uncertainty assessment methodology to a fractured gas reservoir and finally a presentation of a history matching methodology for fractured reservoirs.


INTRODUCTION
After describing in the preceding paper the mathematical formulation of the physics required for simulating flows in multi-scale fractured reservoirs, we now focus on the main specific issue of dual-porosity (dual-permeability) simulators, that is the modeling of matrix-fracture transfers at the simulator cell scale. The optimal way to predict these transfers has no definite answer and various methods are implemented in industrial simulators. However, rather than trying to detail the diversity of such methods, the physical mechanisms driving flows and interactions/exchanges within and between fracture and matrix media are focused on. The purpose is helping the reader in assessing the proposed flow formulas and simulation methods and selecting ones for his own field applications.
The first section presents an overview of matrix-fracture transfer formulations. The second section presents the implementation of geomechanics in fractured reservoir simulation. Quantification of uncertainties and history matching of fractured reservoirs are dealt with in the last two sections. The uncertainty assessment methodology is illustrated on a fractured gas reservoir case.

Matrix-Fracture Fluid Transfer Equations
Before detailing the expression of the matrix-fracture mass transfer flux, the equations (developed in Part 1) for the flow of each component in the fracture and matrix media are recalled to show the contribution of the matrix-fracture transfer flux, F mf kp , in the overall system of equations. Hereafter, the equations for the single permeability option of a dual-porosity model are given. The dual-permeability option would involve an additional transport term within the matrix medium equation (Eq. 2): (1) The expression of matrix-fracture transfer flux is based on the simplified representation of the dual-medium fractured reservoir, as introduced in the preceding paper (Part 1) and shown in Figure 1. Since fracture and matrix grids are identical and superposed, any matrix cell is associated to the In other words, a cell is a discretized element of volume including both fracture and matrix blocks. It is first assumed that, in a given cell, there are N identical matrix blocks of dimensions a, b and c, and that each one behaves as a block located at the center of the cell (Fig. 1). If the cell dimensions are Δx, Δy and Δz, then N = Δx. Δy. Δz/(a.b.c). Secondly, the transfer flux of component k in phase p between a matrix block and the surrounding fractures is expressed as the sum of the fluxes through each of the six faces, i, of the block. We finally obtain the following expression of the matrixfracture transfer flux per unit volume of matrix rock: (3) where a, b, c are the matrix block dimensions. One has to be aware that this block is an equivalent block defined at the model cell scale. That is, all blocks are assumed to be identical and parallelepiped in a given cell. Their dimensions, a, b, c, are the result of an upscaling procedure such as the one developed by Sarda et al. (1997). Discussion of that upscaling issue is found in Bourbiaux (2010).
The mass transfer flux of component k in phase p across a matrix block face i is split into two parts, the transport flux term due to convection-diffusion, f t , and the diffusion-dispersion flux term, f d : Whereas dispersion often predominates under the significantflow-velocity conditions encountered in single-medium reservoirs or in fractures, the low-velocity transport fluxes prevailing in matrix-fracture transfers tend to emphasize the role played by the molecular diffusion mechanism in such transfers.
The following two sub-sections detail the formulas of matrix-fracture transport and diffusion-dispersion fluxes.

Matrix-Fracture Transport Flux
The transport flux term can be expressed either by the direct application of Darcy law between matrix and fractures with the explicit input of matrix block dimensions or the input of a shape factor σ that takes into account the size and shape of the matrix blocks.  Dual porosity representation of a fractured reservoir.
If matrix block dimensions are input, then the Darcy law is applied between the center of the block s and the center of the six faces of the matrix block of dimensions a, b and c. As shown in Figure 2, the face centers in the X, Y and Z directions are respectively x -, x + , y -, y + , zand z + .
The flux, through a face i, of component k in phase p is expressed as in Sabathier et al. (1998): sp is the potential of phase p taken in the matrix medium at the center s of the matrix block. Φ f ip is the potential of phase p taken in the fracture medium at the center i of the face i limiting the block, with i = x -, x + , y -, y + , z -, z + . A i and l i = a, b or c, are respectively the cross-section area and the block length in each direction perpendicular to the face.
As the matrix-fracture flow rate is governed by the matrix permeability, the absolute permeability used for the calculation of this flow rate is K i m , the matrix permeability in the direction i. The formulation of Equation (5) allows taking into account permeability anisotropy of the matrix through different values K x , K y and K z in each direction X, Y and Z. Equation (5) is analogous to a conventional Darcy flow equation, with a specific matrix-fracture transmissibility that is formulated as a function of the dimensions of the equivalent (i.e. cell-representative) "matrix block". Using Equations (3, 5) the matrix-fracture coupling transmissibility between a cell of the matrix grid and the corresponding cell of the fracture grid is then written in the anisotropic case as: The splitting of matrix-fracture transfer flux into different terms referring to the direction of the flux (Eq. 3) gives the possibility to modify or annihilate the contributions of lateral and/or top/bottom faces (i.e. of horizontal and/or vertical transfers) to the overall matrix-fracture flux.
Considering a single-phase quasi-steady-state flow Warren and Root (1963) did not formulate the matrixfracture transfer flux as a function of equivalent block dimensions (Eq. 5) but directly expressed F mf as: (6) with K m the matrix permeability, p f and p m the fracture and matrix pressures, ρ the fluid density, μ the fluid viscosity and σ the well-known shape factor. σ was introduced as a constant matrix-fracture-exchange factor that only depends on the geometry and characteristic size of the matrix blocks assuming a steady-state matrix-fracture transfer. The shape factor dimension is the inverse of a squared length. Then, assuming a parallelepiped block of size a, Warren and Root (1963) determined an analytical solution of the pressure diffusivity in such a block and by identification with Equation (6), derived the following expression of the shape factor: where N is the number of flow dimensions (1 to 3).
As will be shown after for multiphase flow, the matrixfracture transfer is not pseudo-steady-state but involves transient periods of variable importance depending on the flow problem considered. That is, the shape factor is not a constant. This inadequacy of the pseudo-state transfer assumption explains why a certain controversy developed around the definition of the shape factor. Actually, matrix-fracture transfer involves a time convolution between fracture flow conditions (pressure, saturation, etc.) changing with time around the matrix blocks, and the actual flow taking place within the matrix medium. However, constant fracture boundary conditions are implicitly taken into account to formulate the matrix-fracture transfer flux of many dual-porosity simulators. However, this assumption is only made at the scale of a matrix block. Considering that a reservoir model cell contains many matrix blocks, the matrix blocks of a given cell are assumed to be completely surrounded by one of the phases present in the fracture medium, with a proportion of such blocks that depends on the fracture cell saturation of the phase under consideration. That proportion is formulated with the weighting function f p *f detailed in Section 1.4.2. Very soon, many authors became aware of the necessity to improve the very approximate matrix-fracture transfer formulas. They firstly tried to define the best approximate expression applicable in a steady-state transfer equation. All proposed expressions differ more or less from one another according to the approximation made in solving the exact transfer problem at the matrix block scale. Hereafter is presented a non-exhaustive list of such shape factor expressions. Kazemi et al. (1976) proposed a simple formula for the shape factor of parallepiped blocks with an isotropic permeability: The latter results from the finite difference expression of the exchanges, written between the six lateral faces and the block center. In Equation (5), the shape factor is internally calculated from a, b and c values and matrix anisotropy. We obtain for isotropic matrix the same value than Kazemi et al. (1976): As noticed by Thomas et al. (1983) and Ueda et al. (1989), the shape factor from Kazemi et al. (1976) needs to be adjusted by a factor of 2 or 3 for more reliable prediction. Coats (1989) proposed to use a constant equal to 8 instead of 4 in the above expression, probably considering that the average flow length involved in matrix-fracture transfer from block center to a given external face equals a(b, c)/4 and not a(b,c)/2. And he actually found that this approximation is closer from the exact expression derived from the analytical solution of the diffusivity equation than the previous definition given by Kazemi et al. (1976). Lim and Aziz (1995) solved the unsteady state pressure diffusivity equation for matrix blocks of given geometries. The analytical solution of the average matrix block pressure evolution consists in a series of exponential terms with timedependent exponents. As shown later on, the steady-state transfer rate formulation proposed by Warren and Root (1963) actually leads to an exponential evolution of the matrix average pressure with time. By identifying the latter expression to the first term of the analytical series, Lim and Aziz (1995) determined a shape factor constant equal to π 2 , close to 10. Quintard and Whitaker (1996) applied a mathematical technique consisting in averaging, at the overall scale of a matrix block, the single-phase Darcy flow equations describing local flows within the matrix block. The shape factor constant was then determined by minimizing the difference between the two pressure evolutions resulting respectively from the reference volume averaging method and from the steady-state transfer rate formulation. As Warren and Root (1963), but contrary to previous quoted authors, their shape factor expression involves a constant that is not exactly proportional to the number of flow dimensions, for instance 12 and 28.4 for 1D and 2D exchanges, instead respectively of 12 and 24.
Using a random walk approach, Noetinger and Estébenet (2000) determined the same shape factor constants as Quintard and Whitaker (1996).
Table 1 (Bourbiaux et al., 1999;Granet, 2000) summarizes the values of σa 2 , for a cubic block of lateral dimension a, exchanging fluids by one, two or three couples of opposite faces (1D, 2D or 3D flow transfer). Comparison between the volume averaging methods and other approaches can also be found in Landereau et al. (2001). It must be underlined that all shape factor expressions involved in the pseudo-steady formula of matrix-fracture transfer give an approximation of a physical process that is transient by nature. Bourbiaux et al. (1999) quantified numerically the transient nature of matrix-fracture transfer for a pressure diffusivity mechanism. They considered a cubic block of 3.048 m (10 ft) in lateral dimension surrounded by 4 vertical lateral fractures and impervious horizontal limits. The pressure relaxation of a matrix block set from initial time onwards at a pressure different from a fixed fracture pressure was simulated with a finely-gridded matrix block model. The latter gives an accurate prediction, i.e. a reference solution, of the evolution with time of both the matrix block pressure and the fluid transfer rate to the surrounding fractures. These pressure and rate data were introduced into the pseudo-steady-state formulation given above in order to determine the corresponding shape factor. Figure 3 gives the evolution versus dimensionless time: of the shape factor constant, that is σa 2 . As shown in this figure, the latter is not constant but decreases by a ratio close to 8 from the initial time to a late period of time over which it stabilizes to a constant asymptotic value. This asymptotic value, close to 20, actually corresponds to the definition proposed by Lim and Aziz (1995), that is the first-order (late-time) term of the analytical solution for this transient pressure diffusivity process.
In line with the previous study, Rangel-German and Kovscek (2003) derived a new time-dependent shape factor formula for water-oil matrix-fracture transfers. The timedependent shape factor is expressed as the area of the matrix Figure 3 Evolution in time of shape factor σa 2 for a square block (after Bourbiaux et al., 1999). block contacted by the wetting phase divided by the product of the bulk volume of the rock matrix times the distance between the fracture and the matrix loci where the saturation equals the average matrix saturation. The parameters in their formula are obtained from experimental results involving local saturation measurements or from analytical methods. The shape factor converges to the value of Lim and Aziz (1995) for the long times. Van Heel et al. (2008) analytically derived a timedependent (transient) thermal shape-factor that captures the heating of matrix blocks for all time-scales in the case of steam injection in fractured reservoirs. When this transient shape-factor is used in combination with an analyticallyderived viscosity correction (to capture the effect of the temperature profile inside the matrix blocks), the coarse grid dual porosity-dual permeability simulations accurately reproduce the fine grid single porosity simulations.
A number of authors have proposed analytical expressions as a function of time for the transfer functions to match the exponential form of the oil recovery observed in spontaneous imbibition experiments (Aronofsky et al., 1958;Zhang et al., 1996). Kazemi et al. (1992) incorporated these empirical transfer functions in dual-porosity simulators through fast numerical convolution or pseudo-capillary-pressure functions. Di Donato et al. (2007) incorporated transfer functions based on analytical expressions in a dual-porosity streamlinebased model for simulating flow in fractured reservoirs in which matrix-fracture transfers are accommodated as source/sink terms in the transport equations along streamlines. Lu et al. (2008) proposed a transfer function formulation based on 1-D analyses in the literature with a decoupling of the transfer rate into contributions from the different physical mechanisms, fluid expansion, diffusion, and displacement by capillary imbibition and gravity drainage. They used correction factors to capture accurately both the early-and late-time behavior.

Matrix-Fracture Diffusion-Dispersion Flux
Matrix-fracture compositional transfers essentially concern external gas drive processes. The diffusion-dispersion flux of the kth component of phase p across a matrix block face i is formulated as a simple discretized form of Fick's law written between block surface and block center in the selected direction: with D the molecular diffusion coefficient.
The diffusion transmissibility, computed by the simulator with respect to the grid structure (radial, Cartesian, CPG, etc.) as for the permeability transmissibility, is defined as: The exchange area, A i , and the matrix block dimension l i in the direction perpendicular to the face i, are used for this calculation. If water is present in the fracture limiting the matrix block, the area used for diffusion calculation through the matrix block face is corrected for the water fracture saturation value, S f w , as expressed in Equation (10). Implementing such a formula in a reservoir simulator entails some difficulties when one phase, oil or gas, is absent from one of the two media, either fracture or matrix. This situation is encountered during gas drive processes. The composition of the absent phase in one of the two media has then to be defined in order to be able to compute a diffusion flux between matrix and fracture. To that end, both phases are assumed to be in equilibrium at the matrix-fracture interface. Then, the composition of the absent phase is determined through an extended flash of the studied fluid system, calculated through an iterative phase mixing procedure. For instance, if a fully-gas-saturated fracture is in contact with a fully-oil-saturated matrix, then the oil at block boundary is assumed to be immediately in equilibrium with the fracture gas. This oil composition, O s , is determined through an extended flash calculation of the matrix oil composition under the cell thermodynamic conditions. And the matrixfracture diffusion flux in the oil phase is computed by considering the difference between the bulk matrix oil composition and the fracture-block interface composition, O s . Conversely, if oil-saturated fractures surround gas-saturated matrix blocks (for instance during recompression, either by waterdrive or water injection, of an already gas flooded field, leading to a rise of the gas-oil contact in the fractures), then the matrixfracture diffusion flux in the gas phase involves the concentration difference between the matrix-fracture interface equilibrium gas, G s , and the bulk matrix gas.
When both phases exist in the matrix, the diffusion flux is computed for each phase with an exchange area proportional to the matrix phase saturation and the difference of concentration in the matrix and in the fracture.
Another difficulty encountered in the simulation of nonequilibrium gas injection in fractured reservoirs is the assessment of diffusion coefficients. Effective diffusion coefficients can be determined by history matching laboratory diffusion experiments in a matrix block with a reference discretized model of that block (Le Gallo et al., 1997). Lenormand et al. (1998) proposed to calculate the transfer of a component by diffusion as a function of fracture geometry, fluid velocity in the fracture, and fluid compositions in fracture and matrix. The diffusion flux is calculated assuming a uniform fluid saturation at the interface between matrix and fracture and a laminar flow in the fracture. The resulting diffusion flux is given by a simple analytical formula. Calculations were then compared to experiments performed at reservoir conditions with pure methane injected in the fracture and binary or ternary mixtures in the matrix. The results were in very good agreement, without resorting to any adjustable parameters. Regarding hydrocarbon systems, the molecular diffusion process is the single mechanism taking place in the oil-saturated matrix blocks as long as this liquid phase keeps on dissolving some gas circulating in the fractures and expels excess liquid to these fractures as a result of swelling. If the matrix oil cannot or can no more dissolve injected gas components but on the contrary vaporizes light or intermediate components, then a two-phase situation establishes within the block. At that stage, gas saturation is modeled as a block-center average value and no more as a block-boundary value. Such saturation gradients are associated with capillary pressure gradients and may initiate flows within the block until capillary equilibrium is re-established across the blocks. Multiphase transfers then become driving mechanisms for matrix oil recovery. They are detailed hereafter.

Multiphase Transfers
Multiphase matrix-fracture transfers are associated with the movement of fluid contacts in the fracture network limiting the blocks. This movement may be initiated by the depletion process then sustained by injection of water or gas. The transfers associated with a water-oil contact rise and with a gas-oil contact descent are considered separately hereafter.

Water-Oil Transfers
Assuming a quasi-static pressure equilibrium between matrix blocks and fractures, the progression of the water-oil contact within the fractures may entail a transfer of water into the matrix blocks, thanks to capillary forces provided that the matrix medium is water-wet, gravity forces due to water-oil density difference, and sometimes also viscous forces in the case of low-conductivity fractures where the pressure gradient exerts a driving force across the matrix medium.
Among them, capillary forces require a careful assessment because they may or not contribute to the oil recovery from matrix blocks, depending on the wettability of the matrix pore surface. A water-wet matrix spontaneously imbibes water, oil being concomitantly expelled, whereas an oil-wet matrix is not spontaneously invaded by water and produces some oil under the sole effect of gravity forces and viscous forces if significant. Many fracture fields, mostly of carbonate nature, are characterized by intermediate-or neutral-wettability situations where the matrix medium has a limited and low affinity for both the aqueous phase and the hydrocarbon phase. In such situations for oil-wettability cases as well, gravity forces can considerably enhance the oil recovery from matrix blocks.
Capillarity-driven exchanges, denoted as spontaneous imbibition, involve countercurrent flows of water and oil between the heart of the matrix block and the lateral boundaries of that block that are set in contact with the aqueous phase. Additional gravity forces tend to drive water and oil flows in the same direction vertically. They impede the countercurrent production of oil on the bottom faces of blocks and enhance it on their upper faces. For weak capillary forces, gravity forces may even completely override the effects of capillary forces and lead to an upward cocurrent flow with water imbibing the lower boundaries of the blocks and oil produced on the upper boundaries.

Ultimate Oil Recovery
The ultimate saturation of a matrix block surrounded by water is the result of equilibrium between capillary forces and gravity forces at any position within the block. That is, the average saturation of a matrix block limited by water-saturated fractures can be formulated as follows: Although gravity forces are always driving forces for the oil production of a water-surrounded matrix block, these forces have no significant influence on the oil recovery at the matrix block scale in the case of a purely water-wet matrix, because always-positive capillary pressures (i.e. positive over the entire saturation range) are alone sufficient to bring the matrix block oil saturation to the microscopic residual oil saturation value. But for an intermediate-or oil-wet matrix medium, the oil desaturation of the matrix is not complete under the sole action of capillary forces, then the water-oil gravity head exerted locally within the blocks shifts the local water saturation to higher values, thus resulting in the average block saturation given by Equation (12) at equilibrium. As the block height increases, the difference between the final average saturation of the block and the "microscopic" residual oil saturation vanishes, because the capillary retention of oil on the upper limits of the matrix block then represents a decreasing proportion of the total block height.

Kinetics of Imbibition
Let us first consider a pure capillary imbibition process without any significant effect of gravity forces. For a dimensional analysis of the process, we assume that a front progresses from the boundary of the block toward its center, with water and oil flowing in reverse directions behind the front. Note that in countercurrent flow conditions contrary to cocurrent flow conditions, the oil phase necessarily remains mobile behind the front for being produced, that is, the oil saturation remains above the residual saturation, at least in the course of the displacement process. For simplification purposes, we also assume a one-dimension piston-type displacement inside a block of size l in contact with water on one of its face only. Applying mass conservation law and generalized Darcy law in both phases, the fraction of recoverable oil that is recovered at a given time, t, is formulated as: (13) with φ the matrix porosity, ΔS the matrix saturation variation, p cw the capillary pressure drop across the displacing front, and k rw and k ro the relative permeabilities of water and oil when flowing in countercurrent directions.
Equation (13) for , which is based on very simplified assumptions of the displacement process, is given for the sole purpose of physical dimensional analysis. It shows that the countercurrent production of oil varies as the square root of time. Experimental measurements of imbibition actually reveal such behavior during the initial stage of the imbibition process, as long as water has not reached the center of the block then behaving as a semi-infinite medium. The time required for recovering a given fraction of recoverable oil is proportional to the block size squared and inversely proportional to the product of matrix block permeability and capillary pressure, i.e. inversely proportional to the square root of permeability for porous media characterized by the same dimensionless capillary pressure curve (or Leverett function). This flow behavior can be experienced in water-wet fractured reservoirs with a high fracture density and conductivity, where both gravity and viscous flow effects can be neglected.
In other situations, especially neutral, intermediate-or oilwettability situations, the effect of gravity forces becomes preponderant on the kinetics of production in addition to the ultimate oil recovery as stated above. This kinetic impact is considered hereafter in the frame of gas-oil transfers.

Gas-Oil Transfers
Symmetrically to water-oil contact rise in the fracture network, the gas-oil contact moves down during depletion because of gas-cap expansion, segregation of liberated solution gas, apart from possible gas injection. The main difference lies in that the gas-oil capillary pressure, defined as p cg = p gp o , is always positive, since the liquid phase is always wetting in the presence of gas. In addition, the minimum threshold pressure, p c min , required for gas to enter the matrix, determines a minimum block height, H min , below which no oil can be drained from the matrix: (14) It is noteworthy to indicate that the capillary threshold height, H min , may increase by one to two, even three orders of magnitude as the reservoir pressure decreases from the saturation pressure down to low values. This observation stems from the fact that the capillary pressure is proportional to the interfacial tension, which is very sensitive to pressure, especially for volatile oil systems. This pressure dependency may impede the efficiency of a gas drive implemented after fractured reservoir depletion.

Ultimate Oil Recovery
As for the water-oil situation, the ultimate saturation of a matrix block initially saturated with oil and surrounded by gas is the result of equilibrium between capillary forces and gravity forces at any vertical position within the block. The average saturation of a fully-gas-surrounded matrix block is formulated as follows: Obviously, the impact of the threshold capillary pressure on the final average saturation of the block becomes negligible when the block height largely exceeds the capillary retention height, H min .

Impact of Capillary Continuity and of Flow Barriers on the Oil Recovery
The existence of contacts or porous "bridges" across fractures leads to the existence of a partial or total capillary continuity between successive blocks. Then, the effective block height to be taken into account for the estimation of the matrix oil recovery is equal to the height of the stack of blocks in capillary continuity and not the height of individual blocks. The ultimate oil recovery may then be much higher than what would be estimated considering the individual block height, whereas the kinetics of production may be restricted at the low-surface-area contact between blocks. Labastie (1990) performed both numerical and experimental investigations to determine under which conditions capillary continuity actually exists. He concludes that the capillary continuity between matrix blocks is established through porous medium bridges. When non-porous bridges bind the matrix blocks together, a capillary continuity may still be established via the fracture if its aperture is (very) small. Matrix capillary continuity modeling issue could then be dealt with via a fracture capillary pressure and via a matrix permeability/transmissibility reduction accounting for the existence of inter-block bridges. Sajadian et al. (1998) S gc S p dp geq m og g cg cg performed the same kind of experimental study and introduced a critical fracture aperture size, depending on the interfacial tension, the spreading coefficient, the rock wettability and the surface roughness.
Conversely, the existence of flow barriers in the matrix medium, for instance impermeable shaly laminations, may reduce the effective block height. All intermediate situations are frequent where heterogeneities reduce the flow ability while keeping the capillary continuity.

Kinetics of Drainage
Let us consider a pure gravity-driven flow in a matrix block of height l, with a negligible impact of capillary forces, i.e. no capillary retention at the base of the block. Again, the displacement is assumed frontal (or piston-type) but with both phases flowing in the same direction (cocurrent flow). The viscous pressure gradient in the gas phase is also neglected by comparison with the oil phase pressure gradient. Applying mass conservation law and generalized Darcy law in both phases, the recovered proportion of recoverable oil can then be formulated as a function of time as: (16) Equation (16) shows that the cocurrent recovery of oil normalized with respect to the amount of recoverable oil, , varies linearly with time. The normalized recovery rate, , is proportional to the matrix permeability, to the phase density contrast and to the oil mobility: keeping in mind that the vapor phase viscosity effects were assumed negligible.
The normalized recovery rate is also inversely proportional to the block height. This last observation simply stems from the fact that the gas-oil front moves down at the same speed whatever the block height. That is, the maximum kinetics of drainage of a given fractured rock volume is obtained for blocks with a minimum height. However, for recovery purposes, this height must remain at least one order of magnitude higher than the capillary threshold height.

Derivation of Transport Flux
The potential difference in matrix-fracture transfer Equation (5) can be written in a particular way in order to show the respective contributions of the involved physical mechanisms: This enables to split the potential difference into the four following terms (Sabathier et al., 1998): (18) where p cp f(m) is the capillary pressure of phase p, defined as the difference between the pressure of phase p and the oil pressure. Φ f ip is interpolated from the fracture grid values. ρ f* is the average fluid density in the fracture node expressed as: The four terms of the previous equation represent the respective contributions of the different physical mechanisms of transfer, which are: -expansion forces linked to a possible pressure difference between matrix and surrounding fractures; -capillary forces which may either enhance or prevent the oil transfer from matrix to fracture; -gravity forces which are involved in vertical transfers only; -viscous forces linked to fracture flows.
Fracture viscous forces are generally negligible, except when the fracture network is poorly differentiated from the matrix medium and has a low conductivity.
As for most dual-porosity simulators, the formulation described before cannot reproduce the detailed transfer process of fluids at matrix block scale because of the very simplified representation shown above in Figure 1. Therefore, the parameters of the formulation must be adjusted to the solution of a finely-gridded block model. However, this tuning can be based on a physical approach by adjusting independently the different mechanisms with scaling factors, which modify Equation (18) as: More specifically, setting a α factor to zero allows to test selectively the importance of the corresponding physical phenomenon and/or to avoid the computation of negligible physical contributions (the model is generally used with These α values can be tuned to reference experimental tests or single-porosity model results to ensure a higher predictability of the dual-porosity model. Moreover, as described later on, they can be defined analytically for simple matrix-fracture transfer mechanisms. Abushaikha and Gosselin (2008) built a set of reference cases to benchmark the dual-medium models, including the Sabathier et al. (1998) formulation, presented above. They concluded that the latter introduces a better representation of the gravity forces, by splitting the matrix-fracture transfer into a horizontal exchange term and a vertical one, and by using a vertical equilibrium approach. They noticed some weak points in the modeling of capillary imbibition transfer, but the improvements described in the next sections were not tested.

Available Options for Scaling Factors
The α coefficients, which are by default equal to 1, may be used by the user to modify the relative importance of the driving forces, either to test their impact on the matrix block production or to match the reference results given by laboratory tests or by a finely-gridded model. Four physical options can be selected.

Option 0
It is equivalent to α c = α g = α v = 0. This option is recommended for single-phase depletion as it takes into account the compressibility of pore volume and fluids alone.

Option 1 (PC)
It is equivalent to α g = α v = 0. Only the capillary forces are active. This option is recommended for water injection in small blocks.

Option 2 (PC-GR)
This option is equivalent to α v = 0. In addition to capillary transfer terms, exchanges due to density differences between fluids are computed in the vertical direction, as shown above in Figure 2. If several fluids are flowing in the fracture, it is assumed that, in a given cell, each individual block is surrounded by only one fluid and that the number of blocks immersed in gas, oil and water is proportional to the respective saturations (normalized to 1 -S wi ) in the fractures. The average gravity effect results in the ρ* definition. Many comparative tests have been made between options 2 and 3 (cf. below) and have demonstrated that this approximation is acceptable to account for gravity in practical cases, by comparison with other inherent approximations.
This PC-GR option is the most commonly used for field studies. One will notice that this formula involves three fracture-matrix connections: one horizontal connection that lumps the exchanges on all four lateral faces (x -, x + , y -, y + ) and two vertical connections for the bottom and top faces.

Option 3 (PC-GR-VI)
This option takes into account all the mechanisms of multiphase matrix-fracture transfer, that is, compressibility, capillarity, gravity and viscous drive. Six fracture-matrix fluxes are computed on the basis of the phase potential imposed by the fracture on each block face center, "i". As the phase potential is solved at the cell center "s", its value on a given block face center, "i", is extrapolated using its gradient in the considered direction, "s-i". This potential gradient, which results from the fluid flow in fractures, is estimated from a linear interpolation between the potential values in the neighboring fracture cells along direction "s-i" (see Quandalle and Sabathier (1989) for more details).
The reliability of these options is illustrated by the following example from Quandalle and Sabathier (1989). They simulated water injection in a fractured column (Fig. 4a). The column is made of eight 30-feet (9.144 m) cubic matrix blocks surrounded by fractures and is initially oil-saturated. Water is injected in the fractures at the bottom of the column while an equivalent fluid volume is produced at the top. The column was first simulated with a single-porosity model and a standard Cartesian grid where each matrix block is represented by one grid cell actually surrounded by fracture cells.
Results obtained with this single-porosity grid are considered as the reference because the grid reproduces exactly the Warren and Root (1963) ideal model with each matrix block represented by a single cell. The corresponding dual-porosity grid is made of 8 cells with 8 fracture unknowns and 8 matrix unknowns.
Two sets of runs were performed. The first set corresponds to a high fracture permeability of 100 Darcy (98.7 μm 2 ). The second one corresponds to a lower fracture permeability of 10 Darcy (9.87 μm 2 ). The cumulative oil recovery versus time obtained with the two grid systems for both cases are presented in Figures 4 and 5. Fractured column. Cumulative oil production, high-fracturepermeability case (after Quandalle and Sabathier, 1989). The three options described above, PC, PC-GR and PC-GR-VI, were used in the dual-porosity simulations for both cases. Gravity plays a significant role on production for both the high and low fracture permeability cases. Hence, the first option, PC, taking into account only the capillary forces, cannot reproduce the reference single-porosity model results, as shown in Figures 4 and 5. The PC-GR option predicts very well the high fracture permeability case as the viscous effects in the fractures are negligible. On the contrary, the PC-GR-VI option is necessary in the low-fracture-permeability case (Fig. 5).
In most matrix-fracture transfer cases, the dual-porosity model results have to be tuned to a reference solution, given by a fine-grid single-porosity model as above. This tuning is performed via the scaling factors defined before. If such a matching procedure is easier and more satisfactory (from a physical point of view) than using pseudo-k r /p c curves, it nonetheless remains time-consuming for large reservoir models with complex histories. For this reason, semi-analytic procedures were developed. They are described hereafter for two frequent types of transfer, water-oil capillary imbibition and gas-oil gravity drainage, and also evoked for a more complex transfer case involving molecular diffusion and a kinetics of reaction.

Specific Uses of Predictive Scaling Factors
Developments and results given in this section are also found in Sabathier et al. (1998).

Capillary Imbibition
The capillary imbibition mechanism is predominant for water-surrounded blocks characterized by water-to intermediate-wettability properties. If gravity contribution is negligible, the final matrix oil recovery is determined by the capillary equilibrium between fracture and matrix media, that is, by the matrix medium saturation corresponding to the same capillary pressure in the matrix as in the water-saturated fracture (generally zero). This value must be part of the model input p c curve. On the other hand, the kinetics cannot be accurately computed with a single node representation of the matrix block. To account for a gradual advance of water from block boundaries to block center, the capillary pressure gradient (Eq. 5,20) is no more computed with the constant block size value, l, but with a flow distance that is function of the matrix saturation (Sabathier et al., 1998).
The benefit of this improvement has been quantified for a single block of 10 feet (3.048 m) by comparing the dualporosity predictions of water-oil imbibition with the reference solution given by a finely-gridded single-porosity model. Rock-fluid data, taken at a pressure higher than the bubble point, are those used for the Sixth SPE Comparative Solution Project (Firoozabadi and Thomas, 1990) except for the water-oil capillary pressure curve which was modified to be representative of a strongly water-wet system. The fracture was assumed to remain completely saturated with water all through the imbibition process. Figure 6 shows that the kinetics of capillary imbibition predicted by the basic dualporosity formulation (with Kazemi shape factor) is delayed by comparison with the exact kinetics computed with a discretized block model. However, with the improved formulation incorporating a distance-versus-saturation function, the kinetics of transfer is much better reproduced (Fig. 6). Such an improvement has several advantages over the present methods used to tune dual-porosity models. Firstly, it is only based on the physics of imbibition. Secondly, it is internally applied by the model, avoiding any subgridding of the matrix or arbitrary and/or time-consuming tuning methods. The limits of this improved formula lie in the choice of the function tying the distance of exchange with matrix saturation. Cumulative oil production, low-fracture-permeability case (after Quandalle and Sabathier, 1989).

Figure 6
Cumulative oil recovery of a 10-feet cubic block by water-oil imbibition (after Sabathier et al., 1998).

Gravity Drainage
The exact average saturation of a matrix block at the end of a gravity drainage process can be determined by integration of the saturation values derived from the capillarity-gravity equilibrium at each block height, that is, for a gas-oil drainage and a zero fracture capillary pressure: The general dual-porosity formulation detailed in Equation (5) predicts a different equilibrium saturation value,S m* geq , which satisfies the following equation: However, a scaling factor of the capillary term, α c , or the gravitational term, α g , of the matrix-fracture potential difference can be included into the formulation in order to satisfy Equation (21) at equilibrium. In addition, a good approximation of the exact kinetics of drainage can be obtained by introducing a multiplying factor of the vertical matrix-fracture transfer flux, C tz , which is derived from the initial matrix-fracture potential difference responsible for gravity drainage. C tz can be expressed versus the α c scaling factor and the threshold capillary pressure, p m cg min , as: The use of such factors α c and C tz is equivalent to that of a pseudo-function of capillary pressure. The proper values of these factors are dependent on block dimensions and fluid properties which may be subject to significant changes versus pressure during field exploitation. Therefore, a dynamic procedure was incorporated at cell scale into the dual-porosity simulator to compute at each time step the values of the scaling factors detailed before, whatever the block geometry and the density contrast and interfacial tension between phases.
The single-block gravity drainage case used for the Sixth SPE Comparative Solution Project was taken as a demonstration case (Sabathier et al., 1998). As previously for the capillary imbibition case, the reference solution for matrix oil recovery was computed on a finely-gridded single-porosity model with 1000 matrix cells. Oil recovery results versus time are shown in Figure 7. The basic dual-porosity formulation (with Kazemi shape factor) predicts an oil recovery of 52% of initial oil in place, instead of 40.2% which is the actual average matrix saturation inferred from the equilibrium between capillary and gravitational forces or computed with the single porosity model. The improved dual-porosity formulation which incorporates the previous equilibrium calculation gives an exact prediction of matrix oil recovery and quite an acceptable prediction of production kinetics. Cumulative oil recovery of a 10-feet cubic block by gas-oil gravity drainage (after Sabathier et al., 1998).
The limits of that improved formula stems from an adjustment that concerns the final average saturation of the matrix block and the initial matrix-fracture transfer rate alone, without considering the exact transient rates. A more accurate prediction could be obtained by introducing a matrix-saturation dependent scaling factor, but with the remaining question of its general validity for any matrix-fluids system. The necessity of such improvement should also be evaluated with respect to the practical accuracy required for the transient states of matrix blocks.

Complex Transfers
The difficulties are even greater if we deal with multiphase transfers of fluids in thermodynamic non-equilibrium because mass transfer of components occurs at the interface between phases in addition to convective and diffusive transfers occurring within each phase. Bourbiaux and Lacroix (2008) proposed a reliable formulation for simulating multiphase compositional transfers between a matrix block filled with oil and a fracture filled with non-equilibrium gas, without having to discretize the matrix block. The modeling is based on the underlying principles: -a component transfer takes place through the interface between oil and gas so as to establish local equilibrium between both phases: this local equilibrium is obtained instantaneously compared to the time required to establish a global equilibrium between fracture and matrix fluids, by means of diffusion phenomena at block scale; -a molecular diffusion occurs within each phase as a result of the composition difference between the phase interface and the phase taken as a whole; -due to volumetric changes within each phase, the liquid phase expands (swelling) or shrinks (vaporization) within the matrix block. ) ) dp cg m Actually, this formulation assumes that a vapor-liquid front with quasi-static equilibrium at the interface is established from initial time onwards and moves at a speed controlled by the diffusion transfers within each phase. Such a formulation, coupled with the modifications described hereunder, allows simulation of the matrix-fracture exchange flows when the fracture is entirely saturated with gas and the matrix entirely saturated with oil.
For a liquid-saturated matrix block in contact with a gas-saturated fracture, the diffusion flux for either phase (vapor or liquid) is expressed by Equation (10). If matrix and fracture are saturated with different phases, Equation (10) written for vapor phase predicts no flux because concentration C m kp is nil. In the same way, no diffusion flux in the liquid phase is predicted because C f kp equals zero. In addition, the distance of exchange is fixed (l/2) and does not take into account the fluid dynamics described above.
Then, the following improvements were made: -to compute a gas diffusion flux between fracture and matrix, a fictive matrix gas concentration, C m kg , is adopted as long as a vapor phase is absent from the matrix medium. This fictive concentration is the one in equilibrium with the liquid phase composition of the matrix. Reciprocally, the same procedure can be used to compute a diffusion flux in the liquid phase, the fictive concentration being this time a fictive fracture liquid concentration; -as for the capillary imbibition case, the distance of exchange l/2 in Equation (11) is no more constant but variable in order to take into account the progression of the vapor phase into the matrix block. Consistently with these improvements regarding diffusion transfers, reaction kinetics and thermal conduction formulas were also developed to correctly upscale the oxidizinginduced matrix-fracture transfers.
All the modeling features of this complex transfer were validated against single-porosity simulations in 1-D, 2-D and 3-D with a binary system, a ternary system and a complex system largely inspired from the Ekofisk field (Lacroix et al., 2004). These simulations were carried out to identify and model the physical mechanisms controlling matrix-fracture transfers during air injection in light-oil fractured reservoirs, first at the matrix block scale, then at the field scale. They involve reactions in addition to diffusion transfers. Results in 2-D XZ and 3-D can be found in Lacroix et al. (2004).

Relative Permeability Calculation at Matrix-Fracture Interface
The fractures surrounding a matrix block are considered as flow boundary conditions for the block. The matrix block is not subgridded and specific transmissibilities are used to compute the fluxes between the block and the face centers. The basic expression of a phase (volumetric) transmissibility is , as it appears in Equation (5). The absolute permeability, K m , is the possibly anisotropic matrix cell permeability. The area, A, of each matrix block face and the distance of exchange, l, are computed from the block dimensions, which constitute specific input properties assigned to each cell or group of cells. Phase viscosity, μ, is a function of the pressure and composition, the latter variables being taken in the upstream cell. A specific point, developed hereafter, concerns the choice of the relative permeabilities, k r . This choice depends on the phase flow direction, from matrix to fracture or inversely.

Matrix to Fracture Flow
When a phase is flowing out of the matrix to the fractures, its relative permeability corresponds to the saturation of the matrix cell, i.e. that of the upstream cell as usually; however, a physical reason for this choice is that the flow takes place in the matrix medium whereas the fracture only plays as a boundary condition. One has to be aware that standard dualporosity models generally use conventional relative permeability data derived from forced displacement experiments, although spontaneous displacements, such as capillary imbibition, are involved in matrix-fracture exchanges. However, some laboratory studies (Bourbiaux and Kalaydjian, 1990) showed that the coupling between two immiscible phases would be increased in countercurrent flows by comparison with cocurrent flows, hence that distinct cocurrent and countercurrent relative permeabilities, or a full matrix of relative permeabilities, would ensure a higher predictability of matrix-fracture transfers in heterogeneous reservoirs where the magnitude of cocurrent and countercurrent exchanges may change from one reservoir location to another. As most simulators use a single set of k r , the best choice would then be to adopt the k r curves that enable to match a representative laboratory experiment of matrix-fracture transfer. Flow hysteresis may be used as this phenomenon may occur for complex field histories, involving for instance successive depletion and re-pressurization phases leading to various trapped gas saturation values in matrix blocks according to their past saturation history.

Fracture to Matrix Flow
When a phase is entering the matrix from the fractures, the relative permeability would depend on the upstream (fracture) saturation according to the conventional upstream scheme used for k r , whereas the whole flow takes place within the matrix block. For the previous reason, the fractureto matrix relative permeability of a given phase p is: with f p *f a weighting function of the fracture phase-p saturation. The definitions of the matrix relative permeabilities  (24) depend on the nature of the phase as described hereafter.

Gas and Water Phases k *m
rp is the matrix phase relative permeability that corresponds to a complete matrix flood (conditions which are supposed to prevail at the block limits). Consequently: When hysteresis is not activated, the maximum saturations for the gas and water phases are, S gmax = 1 -S wi -S org and S wmax = 1 -S orw -S gc . When hysteresis is taken into account in the matrix, specific values are defined. f *f p is a function of the normalized phase fracture saturation, S *f p , defined as: Taking f *f p = S *f p for gas and water phases, corresponds to the hypothesis that, within a cell, the number of elementary blocks that are surrounded by gas and water is proportional to the mobile phase saturation, which is consistent with the assumption of phase segregation within the fracture network.
Other types of function f *f p may however be specified in this formulation.

Oil Phase
In the case of an imbibition or re-imbibition of matrix blocks by oil, k *m ro is defined similarly to the other phases as: Equation (20) for the matrix-fracture transfer shows that the fluxes may result, on the one hand, from the pressure difference between the fracture and matrix nodes (first line of the equation) and, on the other hand, from the capillarity, gravity and viscosity forces. Oil pressure difference between the fracture and the matrix blocks necessarily leads to the (forced) invasion of the oil-surrounded matrix blocks by this phase; hence the following oil weighting function is assigned to that pressure contribution to the matrix-fracture transfer: This allows a fracture-to-matrix oil flow during reservoir re-pressurization scenarios involving a re-saturation of the fracture network with oil.
For the other fractions of the potential, the oil relative permeability is computed with a weighting function that is function of the following normalized fracture oil saturation: (31) where S o mf is a minimum fracture oil saturation, specified as input, below which the oil re-imbibition is no longer allowed.
As above for water and gas phases, other types of function f¨f o may be given as model input. Moreover, they can be specified separately for horizontal and vertical flows.

Simulation of Oil Re-Imbibition
As already stated, gravity drainage is the main mechanism of oil production in gas-drive fractured reservoirs. However, another physical mechanism has to be taken into account because the oil flowing out of a matrix block generally reimbibes, at least partially, the lower neighboring blocks if already saturated by a vapor phase.
Although experiments can give evidence of that capillarity-driven phenomenon (Lefebvre du Prey, 1976), the necessity to simulate it for a given field case study cannot be assessed a priori but has to be evaluated from the history match of production by gas-oil gravity drainage. Actually, reimbibition delays the oil recovery from matrix blocks. A few methods to reproduce re-imbibition effects on production with dual-porosity simulators are given hereafter. Por et al. (1989) incorporated the block-to-block re-imbibition phenomenon in Shell reservoir simulator by adding connections between the matrix and fracture nodes of superposed cells, and by using connection-dependent relative permeabilities. This way, the matrix oil produced into the fractures at a given cell location can re-imbibe the matrix blocks of the cell located below. The gravity driving force applied on matrix blocks is obtained via a half-block-height shift between fracture and matrix nodes. Fung (1991) used pseudo-capillary potentials and the dual-permeability option to compute matrix-to-matrix flow terms in the vertical direction and thus takes into account the re-imbibition phenomenon. The pseudos are calculated a priori if a vertical equilibrium assumption can be made for the fluid distribution in the matrix blocks. When this assumption is not valid, the pseudos can be determined, as usual, from fine-grid simulations.
Regarding IFP approach, a re-imbibition function may be used to model block-to-block flow interactions within dualmedium cells. Fracture-to-matrix oil re-imbibition is not simulated between neighboring cells; however, the adoption of the dual-permeability option gives a possibility to reproduce the same overall effect, i.e. a "dynamic" oil retention by the matrix medium. Regarding intra-cell re-imbibition, an oil inflow (re-imbibition) rate into matrix blocks is computed in addition to the outflow rate resulting from gas-oil gravity drainage. This inflow rate is estimated as the following fraction α of the outflow rate: where α max is the maximum value of the reimbibition fraction, and f(S * o ) is a vertical reimbibition function that is defined versus the normalized matrix oil saturation S * o .

IMPLEMENTATION OF GEOMECHANICS IN FRACTURED RESERVOIR SIMULATION
Geomechanical effects can be particularly pronounced in some reservoirs, such as poorly compacted reservoirs and fractured or faulted reservoirs, and need in this case to be taken into account. The set of equations underlying a thermohydro-mechanical problem may differ according to the level of coupling between thermal-hydraulic and mechanical phenomena. After space-time discretization, the coupled problem can be written in the following matrix form (Settari and Walters, 1999): (33) where κ is the stiffness matrix, u the displacement vector, L the coupling matrix between mechanical and flow unknowns (respectively deformations and pore pressures), P is the vector of reservoir unknowns (i.e., pressures, saturations and temperature), F is the force vector resulting from the applied boundary stress tensor, L T is the transposed matrix of L, E is the flow matrix and R is the right hand side of the flow equations including the source/sink terms. The symbol Δ t represents the change over time step, i.e.: with n the index of time discretization. Note that in the conventional reservoir simulation notation (Aziz and Settari, 1979) matrix E is decomposed as E = T -D where T is the symmetric transmissibility matrix and D is the accumulation (block diagonal) matrix.
For fractured reservoirs, this set of equations is written for the matrix (denoted as "rock" in the present section) medium (grid) and for the fracture medium (grid). Equation (33) accounts for the geomechanical equilibrium of the considered medium, whereas Equation (34) represents the fluid mass balance of that medium. Note that, to simplify the presentation, the stiffness matrix κ and the flow matrix E are considered here to be linear operators. In general however, E and κ are non-linear operators accounting for the non-linear nature of the hydro-mechanical problem to be solved.
Equations (33) and (34) are coupled through the coupling matrix L. On the one hand, pore pressure, saturation and temperature changes in the rock and fracture systems, Δ t P, affect the stress equilibrium equation. On the other hand, geomechanical strain changes, Δ t u, have an impact on flow through the modification of the flow transmissibility matrix T, resulting from: -rock permeability variations; -fracture and fault permeability variations due to normal and shear displacements. Three (mechanical) coupling levels are usually defined (none, partial and full). Settari and Walters (1999) describe the different coupling levels that can be used to solve the whole set of Equations (33-34): • No mechanical coupling: resolution of thermal and mass fluid balance in the reservoir without any calculation of stress equilibrium (conventional reservoir simulation); • The partially-coupled approach is based on an iterative explicit coupling between a conventional reservoir flow simulation (mass fluid balance) and a geomechanical simulation (mechanical equilibrium). In this methodology, the pore pressure and temperature variations in the rock and fracture systems computed by the reservoir flow simulator are used as input load data of the geomechanical simulator for solving stress and strain variations in the reservoir. Stress-dependent reservoir properties (porosity, matrix permeability, fracture conductivity or equivalent fracture network permeability) are then updated in the reservoir model. From this point, two configurations can then be encountered. The first one -the explicit method -consists in using the updated stress-dependent reservoir properties to simulate the next time step of the reservoir flow simulation. It is well suited for low levels of coupling. The second one -the iterative method -consists in solving again the present time step of the reservoir flow simulation using the updated stress-dependent reservoir properties to improve the consistency between flow simulation and geomechanical simulation at the end of the step. This partially-coupled approach has the advantage of being flexible and benefits from the high developments in physics and numerical techniques of both the reservoir and the geomechanical simulators. Depending on the studied physical mechanisms, this coupling methodology can be iterative or not but in all cases it must be designed in order to ensure a consistent and stable process; • The fully-coupled approach simultaneously solves the whole set of equations (mass fluid balance and mechanical equilibrium) in one simulator. It leads to consistent descriptions but published works show that the hydraulic or geomechanical mechanisms are often simplified by comparison with conventional uncoupled geomechanical and reservoir approaches.
Research activities on the integration of geomechanics in reservoir simulation have been mainly focused on singleporosity reservoirs with the aim of predicting compaction and subsidence of weakly-compacted reservoirs (Longuemare et al., 2002). Only few references (Koutsabeloulis et al., 1994;Guttierez et al., 1997;Lewis and Ghafouri, 1997;Lewis and Pao, 2002;Marmier et al., 2006;Bagheri and Settari, 2008) are related to the integration of geomechanics in doubleporosity reservoir simulation. The assessment of such effects is usually carried out on very simple reservoir configurations and the applied methodology is not yet commonly implemented for the study of a real field case. However, if they actually occur, fracture-deformation phenomena have a great impact on the performance of naturally fractured reservoirs because the fracture permeability is highly sensitive on the fracture aperture, as shown by Bagheri and Settari (2008). Simulating a field depletion case, these authors showed that the fracture permeability was reduced as the result of the closure of fractures due to an increase in the effective stress with pore pressure depletion. That permeability reduction was up to four times of the initial value, especially at the reservoir top. As a result, the oil rate was predicted to decline more rapidly than in the case of a simulation ignoring geomechanical effects. On the contrary an increase of permeability was observed during the production of the huge gas fractured field of Lacq Profond in France (Rolando et al., 1997;Paux and Zhou, 1997). Surprisingly, together with the depletion of approximately 600 bars due to the 40-year gas production, a significant increase in productivity was recorded at the wells during the same period. This phenomenon was explained by the poro-elastic deformation undergone by the Lacq reservoir since the beginning of production, as evidenced by measurements of microseismic activity within the field. An average variation of permeability with depletion (20 times greater than the initial value) was input in the dual-porosity model and resulted in a satisfactory match of the 40-year period of production.

UNCERTAINTY QUANTIFICATION
Both fracture and matrix flow parameters have a great impact on the hydrocarbon recovery from fractured reservoirs. The large number of parameters and the high uncertainties on some of them entail potential risks in the assessment of fractured reservoirs recovery performance. The statistical theory, and especially the experimental design approach (Feraille et al., 2003), is clearly well-suited to identify the main uncertain parameters influencing production, to quantify their respective impact on production forecasts, and thus to provide a sound technical support to reservoir management decisions.
Let us consider a production response of interest (typically the cumulative oil production, the recovery factor, the gas-oil ratio), and some uncertain deterministic reservoir parameters that can be either "physical parameters" or "production constraint parameters". The principle of the experimental design consists in defining a minimum number of scenarios to be simulated (i.e. a minimum set of "experiments" to be performed) in order to estimate the production response of interest as a function of the uncertain parameters. The uncertainty assessment of this production response then involves the following two steps: • Identification and sorting of the parameters that mostly influence the production response. This step is crucial since it allows to eliminate parameters that have a negligible impact on the response and to focus on the influential ones. This uncertainty analysis is synthesized in a graph, denoted as the Pareto plot. The Pareto plot gives the relative contribution of each uncertain parameter on the response. All the uncertain parameters that have an influence magnitude greater than a given threshold are considered as influential, while the remaining terms are interpreted as negligible. The use of Pareto plot will be illustrated in the example detailed hereafter; • Construction of a model (generally a polynomial function), denoted as the "proxy" model, that expresses the response of interest as a function of the most influential uncertain parameters. Depending on the intrinsic complexity of the response, the latter may be represented by a simple linear function of uncertain parameters, or may on the contrary require to be fitted to more complex polynomial functions. The quality of the proxy model is quantified in terms of accuracy, i.e. to which extent the proxy model fits the simulation results of selected scenarios (the results of selected "experiments") and of predictability, i.e. to which extent the proxy model predicts a field response simulated with parameters chosen at random within their range of uncertainty. The proxy model can then replace the reservoir flow simulator (the "experiments") to quantify a probabilistic field response in the pre-defined uncertainty environment. This probabilistic response is obtained using a sampling method, such as Monte-Carlo, with the fast model approximation of the field value. Therefore the risk on field value due to the influential parameters is quantified as a probability density function. The associated percentiles (P10, P50 and P90) can then be easily found.
An application of this uncertainty assessment methodology to a fractured gas reservoir is presented hereafter and demonstrates the interest of such an approach.

Background
The presence of natural fractures in tight gas carbonate fields often leads to early water breakthroughs leading to premature well shut-ins and to a low field gas recovery. The recovery assessment from such fields requires to take into account the role played by water imbibition of the matrix which, depending on the fracture density and rock properties, can significantly delay water breakthrough. Assessing the impact of such spontaneous imbibition phenomena on the field recovery requires first, experimental measurements (Egermann et al., 2004) to quantify the involved matrix petrophysical properties, then field-scale simulations of gas production on a representative model of the reservoir.
The simulated reservoir example, considered as representative, is a fractured carbonate reservoir with a tight vuggy and dolomitized matrix, of less than 1 md, and a strong underlying aquifer in a tilted structure. The effective fracture permeability is in the order of 20 mD. The reservoir thickness is 50 m. Gas is produced by a single up-dip well perforated over the entire reservoir thickness. Depletion is performed at a constant rate, such that the Original Gas In Place (OGIP) would be recovered after 10 years of water-free production. Production is however stopped as soon as some water, corresponding to a water-gas ratio of only 5 × 10 -6 vol/vol at surface conditions, is produced, that is before 10 years. Fractures provide the main contribution to well productivities but are also responsible for matrix gas bypassing and early water breakthroughs leading to premature well shut-ins.
The matrix petrophysical properties measured on a set of representative samples were also input in the simulation model. Correlations were formulated between these measured properties and were then introduced in the simulation model in order to minimize the number of uncertain parameters in the uncertainty analysis. For the case under consideration, satisfactory linear correlations could be derived between the porosity and the logarithm of permeability with a φ-K law: φ = a × ln(K) + b, and between the gas relative permeability end-points (residual gas saturation and relative permeability values) and the initial water saturation. A reference matrix capillary pressure curve was also defined for given porosity and permeability values and was modified accordingly (i.e. using the same Leverett function) when the latter varied from one simulated scenario to another.

Uncertainty Study
The question raised was to assess the impact of uncertainties of both the fracture parameters and the matrix petrophysical properties on the gas recovery at water breakthrough. The COUGAR TM (2008) uncertainty assessment software was used for this study.

Defining the Uncertain Parameters Space
The first step consists in an inventory of the involved uncertain parameters and in an estimation of the range of values that may be covered by each of these parameters. In the present case, 15 possibly-influent uncertain parameters were identified. They are listed in Table 2, with their respective uncertainty ranges. These ranges were defined according to the laboratory experimental results or from available literature data otherwise.

Setting Up an Experimental Design
As explained before, this step aims at best sampling the uncertain parameters space. A first-order experimental design was selected. That is, it was assumed that the contribution of a given uncertain parameter to the simulated response is linear (i.e. quadratic effects are ignored). The possible impact of interactions between the uncertain parameters was also 279 Reference maximum capillary pressure 10.0 15.0 bars investigated but the number of selected simulations was insufficient to interpret each interaction. Therefore the interactions were gathered in different interaction groups, named "alias groups", and the impact of each "alias group" on the production behavior was quantified. Fifteen "alias groups" were defined. For instance, the first "alias group", named PERMF:POROSF11, includes several interactions between PERMF and POROSF11, BLOCK12 and KRWM11, etc. The content of each "alias group" is automatically and optimally defined by the uncertainty software. With such conditions, the experimental design required the simulation of only 32 scenarios to assess the field gas recovery in the uncertainty space defined by 15 uncertain parameters. Each of these 32 simulated scenarios represents a given combination of minimum, mean and maximum values of the uncertain parameters normalized between -1 and +1 values.

Sensitivity Study
The results of these properly-selected simulations ("experiments") are summarized in Figure 8. The recovery factor, defined as the cumulated gas production at well shut-in divided by the Original Gas In Place (OGIP), is shown as a function of the water breakthrough time. Note that the proportionality between the recovery factor and the water breakthrough time is only the result of a gas production rate taken in proportion to the OGIP. The main observation is that recovery values range from 11 up to 85%, that is, the recovery factor is very sensitive to the uncertain parameters under consideration. The mean value, equal to 43%, looks realistic for such fractured gas reservoirs where early water breakthrough leads to rather low recovery efficiencies. Recovery factor versus water breakthrough time.
Next step consists in explaining the origin of this high sensitivity. Figure 9 represents the Pareto plot for the recovery factor. Each uncertain parameter is represented by a bar with a length that is proportional to its influence on the studied response, and with a pale blue (resp. dark blue) colour showing the positive (resp. negative) impact on the response. For instance, the PERMF uncertain parameter is represented by a dark blue bar showing that the gas recovery factor decreases when PERMF increases. It is ranked fourth according to the influence on the gas recovery factor, with a contribution to the recovery factor variability equal to 7.6%.
Among the 15 parameters under consideration, the fracture porosity (POROSF11) is by far the most influent one with a positive effect. It explains 23% of the variability of the recovery factor. Then, each of the following three parameters contributes for 5 to 10% to the variability: -the matrix maximum relative permeability (KRWM11), with a positive effect; -the fracture permeability (PERMF), with a negative effect; -the block size (BLOCK12), with a negative effect; as well as three alias groups involving the fracture permeability and other uncertain parameters.
This sensitivity diagnosis can be easily interpreted. Indeed, the recovery factor is higher when: -the fracture porosity is increased since most of the fracture gas is recoverable; -the fracture permeability is decreased as such a decrease slows down water progression to the production well, leaving more time for its imbibition into the matrix medium; -water imbibition into the matrix medium and gas expulsion to the fractures are accelerated thanks to a high matrix relative permeability of the water phase, and to a large contact area between small blocks and fractures. From the previous considerations, it follows that very low recovery factors, i.e. very early water breakthrough times (less than 700 days), are obtained for cases combining low fracture porosity and high fracture permeability.
Note that the coefficients related to the matrix petrophysical correlations are not influential for the recovery efficiency, meaning that the laboratory measurement uncertainties were within a range having no significant impact on the gas recovery factor.

The Proxy Model of the Gas Recovery Factor
Taking profit from the sensitivity analysis results, the experimental design technique was applied again to determine a proxy model. The latter was used to obtain a probabilistic estimation of the gas recovery factor in the uncertainty space under consideration. The proxy model building involved two steps: -a second experimental design focused on the main parameters, including the parameterization of a polynomial function modeling the gas recovery factor as a function of those parameters (the proxy model);  Table 2; -it involved both the linear impact and the quadratic impact of each of the above eight parameters, as well their mutual interaction, that is, a full second-order experimental design was used. This more accurate experimental design included 81 reservoir simulations. The corresponding Pareto plot for the recovery factor is shown in Figure 10.
The main results of these 81 runs, i.e. the recovery efficiencies, were matched to a "Response Surface Model", i.e. a polynomial function of the eight uncertain parameters, including interaction and quadratic terms. The most influential parameters remain the same as those identified during the sensitivity study.
The response surface fits the 81 scenarios very satisfactorily, as indicated by the accuracy coefficient value of more than 95%, with respect to the actual responses of these sce-narios. In addition, it can reliably predict another scenario among the studied uncertainty space, with predictability also close to 95%, with respect to the actual response of such a scenario. Hence, it was suitable for a probabilistic assessment of the gas recovery factor.
The Monte-Carlo sampling technique was then used to determine the probability curve of the recovery factor. In the absence of any specific information, a uniform law was selected for each uncertain parameter, that is, any value between the minimum and the maximum value was assumed equiprobable. Ten thousand random scenarios were taken to sample this uncertainty space. Note that the use of uniform laws maximizes the resulting uncertainty on the recovery factor since the value of each uncertain parameter is taken at random over its entire uncertainty range of values.
The resulting cumulated probability curve for gas recovery factor is shown in Figure 11. The percentiles P10 (respectively P90) obtained from the cumulative probability function are the recovery factor values, equal to 53% and 20% respectively, that have a probability of 10% (respectively 90%) to be exceeded.
The corresponding probabilistic distribution with the percentile values P10, P50 and P90 added is shown in Figure 12. The gas recovery factor is poorly-defined, as it ranges from 20 to 53% if we discard the 10% lowest and the 10% highest recovery scenarios. Sensitivity diagnosis for recovery factor.
In conclusion, for the fractured field under consideration, all the parameters determining the amount of fracture gas available for production were found to control the value of the gas recovery factor: a high fracture porosity, a low fracture equivalent permeability as well as matrix medium parameters increasing the kinetics of spontaneous imbibition phenomena, such as a high water relative permeability and a small block size.

Conclusions
It is necessary to underline that uncertainty diagnosis has to be carried out with great care regarding the consistency between the analysis tools and assumptions and the intrinsic nature of studied uncertainties. The questions to answer to that respect are the following: -is a linear Response Surface Model (RSM) accurate enough? Is a quadratic design needed? For each parameter? -is a non-linear design needed? For all responses?
Actually each field case study involves some specific features and the response surface accuracy can hardly be assessed a priori. Examples show that the uncertainties on many reservoir case studies can be predicted with a sufficient accuracy using classical experimental designs (linear or quadratic). Very heterogeneous reservoirs or non-linear responses with threshold effects (such as the watercut or the   Probabilistic recovery factor with uniform laws (density of probability as a function of recovery efficiency).
gas-oil ratio) may even sometime be accurately modeled with a classical, i.e. linear or quadratic, experimental design. On the contrary, a non-linear sensitivity behavior may result from the modification of parameters expected to have a linear impact on the studied response. Combined impacts of multiple parameters are most probably the reason for such an unpredictable behavior of the reservoir.
The way to assess the response surface accuracy is to carry out confirmations runs with additional simulations and to draw a cross plot between the Response Surface Model predicted values and the actually-simulated ones. If a discrepancy is observed, a refined experimental design is implemented, which may be quadratic instead of linear, non-linear instead of linear (Scheidt and Zabalza-Mezghani, 2004).
On a general standpoint, the presence of fractures in a reservoir introduces additional sources of uncertainty on the field production. This justifies the adoption of a probabilistic approach of the field production based on a careful sensitivity diagnosis of uncertain parameters with dedicated computation techniques. Moreover, quantitative uncertainty analysis should constitute the preliminary step of the field history matching in order to determine the most relevant parameters impacting fractured field production.

HISTORY MATCHING
As mentioned in the previous section, a great amount of uncertainties is associated with the numerous data entered in fractured reservoir models. These uncertainties may have a great impact on the reliability of the production forecast. The known production history of fractured fields can be used to improve the quality of the prediction by means of a history matching process. The values of some selected data of the reservoir model are then iteratively modified to match the simulation results to the field-measured production data. Mathematically, the history matching procedure can be seen as an optimization problem that consists in minimizing an objective function quantifying the mismatch between simulated and observed production, and also 4D seismic data if recorded. This formulation allows using automatic history matching processes instead of a manual one.
Adjusted geological variables are mainly permeabilities, porosities, transmissibilities, anisotropy and also matrix block size in the case of fractured reservoirs. Traditional history matching is done on the field flow model, resulting from an upscaling procedure, and not on the geological model. However, this geological model constitutes the most reliable information support of matched variables since it results from the integration of geological, geophysical, and well data through conditional geostatistical simulation (Deutsch and Journel, 1998). For this reason, an innovative integrated methodology for constraining both the 3-D geological model and the field flow model to well data, production history and 4D seismic data was developed (Fenwick and Roggero, 2003;Mezghani et al., 2004). The proposed approach allows to history match complex fractured and non-fractured reservoir models in a consistent way by updating the originallyuncertain field pieces of information rather than up-scaled parameters that intrinsically include modeling approximations and can hardly be downscaled to the original field input. Advanced parameterization techniques are used to modify either the geostatistical model directly or the fluid flow simulation parameters in the same inversion loop.
As a preliminary step, the relevant inversion parameters have to be selected and their respective impacts on production quantified according to a sensitivity study based on the experimental design technique (Feraille and Roggero, 2004). Then, the history matching itself is performed with the most significant parameters using an automated inversion procedure. In this step, the gradual deformation method (Hu, 2000) is used to constrain geostatistical models while preserving their spatial variability (reproducing their covariance). This technique may be combined with other inversion methods in order to adjust other deterministic parameters simultaneously. A gradient-based optimization is performed to minimize the objective function. Gradients are computed by numerical approximations and are updated during the optimization process. This makes possible the history matching of any kind of parameter, for single-porosity as well as dualporosity reservoirs. Finally, the history matching procedure is repeated with several model realizations to quantify uncertainty on production forecasts.
Using a similar approach, Jenni et al. (2004,2007), applied the gradual deformation concept to history match stochastic object-based models of large-scale (field-scale) fractures, namely (sub-)seismic faults and fracture swarms. Such stochastic models are constrained by fracture/fault-related seismic attributes, fault-related strain field, structural information (curvature, etc.), etc. They used an upscaling procedure for performing fluid flow simulations in the presence of networks of large-scale fractures. To deal with this situation, a specific gradual deformation method was designed to gradually move and/or deform the stochastic fractures in the reservoir field while preserving their consistency with the geostatistical constraints, i.e. the distribution of fracture density and orientation, and the location of seismic fractures. The location of sub-seismic faults can then be constrained to field production data using an inverse simulation procedure, as described above. Such a methodology was applied to a water-drive North African field, where the presence of subseismic faults was found to determine the water breakthrough time and the watercut evolution of wells. The initial realization of the stochastic fracture network and the final realization obtained after a global optimization are presented in Figure 13. Some of the wells were well matched, as shown in Figure 14 for well P1. To further improve the match of each well, a local calibration was performed. Accordingly, Figure 13 Initial realization of the stochastic fracture network (left) and final realization of the stochastic fracture network after global optimization (right) (after Jenni et al., 2004).

Figure 14
Comparison of the water-cut curves (function of time in years) of the reference network (--), the initial realization (--) and the realization after global optimization ( ) (after Jenni et al., 2004). the reservoir field was subdivided into four zones, each zone containing an injector and a producer. The local calibration was then performed sequentially zone by zone, giving a satisfactory match for all wells.
This application shows that an object-based stochastic model of large scale fractures, constrained only to orientation and density maps of fracture sets, leads to flow predictions that remain far from the actual field behavior. The history matching procedure based on a gradual deformation of fracture objects gives satisfactory results while preserving the underlying geological features.
The present methodology is still limited to large-scale vertical fractures that cross the entire reservoir. Extension to completely three dimensional fractures is currently underway. Another important issue to be addressed in future works concerns the simultaneous match of main fracture properties.

CONCLUSIONS
Naturally fractured reservoirs are among the most complex hydrocarbon reservoirs encountered, in terms of geological characterization and flow behavior. These reservoirs hold a significant share of the remaining oil in place worldwide, and their contribution is expected to increase in the future development projects. As a major step of the integrated reservoir assessment workflow, reservoir simulation has to provide reliable answers to the engineers in charge of the production of fractured reservoirs. This paper and the preceding one gave an overview of the currently-used numerical simulators and know-how regarding fractured reservoir simulation. This paper is more particularly focused on the modeling of matrixfracture transfers. That topic is extensively reported in the literature, especially the controversial question of shape factors. Actually, despite the availability of reliable upscaling solutions for certain matrix-fracture transfer mechanisms, improvements are still needed for complex transfers, as involved in EOR processes.
Geomechanical effects are not often taken into account in fractured reservoir studies but they may have an impact, favorable or not, on the oil recovery if fractures are stresssensitive. Nowadays, the geomechanical unknowns are most often solved separately from the fluids flow ones, for CPUtime limitation. However, the trend is toward the simultaneous resolution of intrinsically-coupled geomechanical and flow systems, provided that efficient numerical procedures are developed to permit full-field studies.
History matching is one crucial step in reservoir studies and much effort has been made to automate the process. Newly-developed inversion techniques allow matching not only the reservoir model but also the underlying geological one. These techniques can be applied to dual-porosity simulation models to calibrate specific fractured reservoir parameters, such as the fracture permeability or the block size. In addition, gradual deformation methodologies give the possibility to match the geological fracture network model to field production data.
Complementarily to history matching, uncertainty assessment is a necessary step to strengthen the reliability of fractured reservoir studies. Recently-developed analysis Realization after global optimization methods and software can provide at low cost quantitative estimates of the impact of fracture and matrix parameters on oil recovery, along with probabilistic recovery curves. To conclude, a great contribution from these diagnosis and simulation methodologies is expected in the near future for the development of the challenging oil reserves that fractured reservoirs are still holding.