Numerical Simulation of Three-Phase Flow in Deforming Fractured Reservoirs

The mathematical formulation of a three-phase, three-dimensional fluid flow and rock deformation in fractured reservoirs is presented in this paper. The present formulation accounts for the significant influence of coupling between the fluid flow and solid deformation, an aspect usually ignored in the reservoir simulation literature. A Galerkin-based finite element method is applied to discretise the governing equations in space and a finite difference scheme is used to march the solution in time. The final set of equations, which contain the additional cross coupling terms as compared to similar existing models, are highly nonlinear and the elements of the coefficient matrices are updated implicitly during each iteration in terms of the independent variables. A field scale example is employed as an alpha case to test the validity and robustness of the currently formulation and numerical scheme. The results illustrate a significantly different behavior for the case of a reservoir where the impact of coupling is also considered.


INTRODUCTION
Productivity of a reservoir, to some extend, depends on the microseismic fissures that exist in the oil formation.The more intense the extend of fissuring, the more economic the reservoir is as fissures promote flow in the oil bearing strata.Due to the random statistical distribution of the fissures, direct numerical modeling of each individual fissures is a formidable task due to the large number of data sets involved in one hand, and the exceedingly high computational cost on the other.Hence, a continuum approach based on the double porosity model introduced by Barenblatt et al. (1960), and later extended by Warren and Root (1963) is more suitable in this scenario.
Over the past few decades, double porosity models with varying degree of sophistication have been developed to simulate flow through fractured porous media.Taking into account the deformability of the oil formation, Aifantis and co-workers (Aifantis, 1985;Wilson and Aifantis, 1982;Khaled et al., 1984) extended Barenblatt's model to what is now known as deformable double porosity formulation.Alternative of Aifantis' formulation were also given latter by others in this area, e.g.Valliappan and Khalili (1990), Cho et al. (1991), Bai et al. (1993), Berryman and Wang (1995), Ghafouri and Lewis (1996), Chen and Teufel (1998), in chronological order.The majority of these works, however, concentrated on single-phase flow phenomena in a deformable geo-material, which has its special utility in civil engineering application like footing analysis etc. Direct application of such models to a deformable oil reservoir is an idealization far from realistic as far as the petroleum industry is concerned.More recently, Lewis and Ghafouri (1997) and Bai et al. (1998) have extended the double porosity concept to include multiphase flow.To the authors' knowledge, the models developed by these two groups represent the only existing simulators able to analyze multiphase flow in deformable fractured geo-material, with particular emphasis on petroleum reservoir mechanics.However, a few inadequacies exist in these models.The formulation given in (Lewis and Ghafouri, 1997) neglected the fracture deformation while those in (Bai et al., 1998) only considered a two-phase oilwater problem.
The objectives of the present paper are threefold.First, to eliminate the deficiencies in the present multiphase double porosity model.Second, to recast the governing equations based on a multiphase flow formulation and finally, to present a numerical analysis of the proposed model.

PRELIMINARY CALCULATIONS
The model presented here conceptually consists of two separate models: -the deformable skeleton; -and the multiphase fluid flowing through the porous media.
The skeleton description is based on elastic theory, while the flow model is based on the double porosity concept.Two overlapping flow regions are considered, one representing the fissured network and the other representing the porous blocks.It is postulated that an elementary scale exits, which is sufficiently large for the formulation to be valid.In the following derivation, the subscripts α = 1 and 2, refer to the porous block and fissured block, respectively; and subscripts s, w, o and g stand for the skeleton, water, oil and gas phases.Other symbols and notations will be defined locally as they appear in the script, and the assumptions stated explicitly as we progress.
The porosities for both flow regions is defined as: (1) where Ω α is the pore volume, and Ω b is the bulk volume.We use the notation: (2) to denote the total derivative of a quantity and v s is the barycentric velocity of that quantity.Small perturbation is tacitly hypothesized assumed throughout the paper.In a multi-phasic domain containing oil, water and gas, conservation of mass required that: (3) The interacting motion of each phase (not including the skeleton constituents) can be linked to the evolution of the partial pressures of each phase to their saturation values.To this end, a capillary pressure-saturation curve is used, which is an approximation of both the drainage (drying) and the imbibitions (wetting) curves.In a saturated oil reservoir, the fluid pressure values at any point are related by their capillary pressure relationships.In general, capillary pressure is defined as the difference between the nonwetting and the wetting phase pressures, i.e.: (4) Therefore, for water-wet oil reservoirs, the following expressions are used: For oil-water system: (5) For oil-gas system: (6) From Equations (3), ( 5) and ( 6), the following saturationpressure relationships can be derived: where: (10) and: ( 11) are the slope of the S ςα vs. p c plot for an oil-water and an oilgas system, respectively.
To relate the skeleton motion due to the fluid, the effective average pore pressure is defined, where it is an average pore pressure weighted by the saturations for both the matrix and fractures systems viz: (12) On incorporating Equations ( 7)-( 9) into ( 12), the following equation is obtained for the average pore pressure in a three-phase flow regime: (13) where:

Skeleton Deformation
The linear momentum balance equation for the skeleton can be written as: (17) where σ σ is the total stress, and f is any arbitrary external and/or body force per unit volume.The total stress σ σ can be expressed in terms of the effective stress σ σ' and the average pore pressure pα , according to Biot's effective stress law as: (18) In Equation ( 18), the phenomenological constants α αα are termed the modified Biot's coefficients (or modified Biot-Willis parameter).As pointed out by Berryman and Wang (1995), the form given by Equation ( 18) is not seriously in doubt.An important characteristic of this equation is that each point in space now has two weighted fluid pressures associated with it, and therefore these pressures are not the true microscopic pressures in the fluid, but the average over some representative volume.
To obtain an explicit expression for the phenomenological constants in terms of physically measurable parameters, one can use Betti's reciprocity theorem.From the analysis carried out in (Khalili and Valliappan, 1996) for an isotropic elastic medium, it can be shown that: (19) in which K, K p and K m are the drained bulk moduli of the skeleton, fissured block and the solid constituent, respectively.It is not difficult to show from Equation ( 19) that: (20) where is the well-known Biot coefficient for a nonfissured material.This shows that Equation ( 18) lies within the framework of Biot poroelasticity (Biot, 1941).It is worthwhile mentioning that if the volume of fissures is reduced to zero, i.e.K p = K, then Substituting Equation ( 18) into (17), one obtains: Here, the issue of whether Equation ( 21) should be cast into an incremental or finite form will not be discussed.For a more detail treatment of this, we refer the reader to (Schrefler and Gawin, 1996).

Flow Model
For mltiphase flow in porous continua, the following generalized mass balance equation at reservoir conditions holds (Pao, 1998;Pao et al., 2001): (22) in which the L(ς) is the logical operator, defined as: The subscript STC denotes the quantity evaluated at stock tank conditions, and v ςα is the absolute velocity of the fluid phase ς in a medium α.The density of each phase may be expressed at the STC via the relation: (24) where B ς is the formation volume factor of the fluid phase ς.In Equation ( 22), R soα is the solution oil-gas ratio in the medium α, and Γ ςα is the so-called leakage function.Its expression emanates from the fact that during some infinitesimal time dt, a fluid of volume: (25) per unit of bulk volume is transferred, at reservoir conditions, between the porous block and the fissured network.Substituting Equation ( 24) into (22) yields: (26) The relative velocities of fluid phase ς in the medium α can be described by: (27) On the other hand, the relative velocity of the fluid can be linked to Darcy's law via: (28) in which K ςα is termed the effective permeability of phase ς.In a multiphase system, the simultaneous flow of the fluid will cause each fluid to interfere with the flow of the other fluid components.This effective permeability value must be less than the single-fluid permeability of the medium.Then, one can define the relative permeability as: (29) Substituting Equations (2), ( 27), ( 28) and ( 29) into (26), one obtains: From Equation (30), various alternative formulations exist and one can recast the mass conservation equations in terms of pressure, capillary pressure or saturation (Aziz and Settari, 1979).From a literature survey, it was found that Schrefler and Simoni (1991) made an extensive search on this issue and concluded that, for a nonfissured two phase flow problem, the best convergence was found using the combination of displacements and pressures as primary unknowns.Based on their conclusion, Equation (30) will now be recast in terms of fluid pressures as primary unknowns.Note that the time derivative of the formation volume factor and the solution oilgas ratio can be rewritten as: (31) and: (32) which pose no serious difficulty.Therefore, attention is now focused on the porosity term.Making use of the definition in Equation (1), we have: Equation ( 33) states that the increase in porosity is due to the increase in the pore volume minus the increase in the bulk volume.Note that the second expression on the RHS of Equation ( 33) is actually the volumetric strain of the skeleton, i.e.: (34) Upon substitution of Equation ( 34) into (30), and after some manipulation, we have the following eight unknowns (p w1 , p o1 , p g1 , p w2 , p o2 , p g2 , Ω 1 , Ω 2 ) but only six equations.One way to obtain the two additional equations is by relating the unknowns Ω 1 and Ω 2 with the primary field variables p ςα , and the skeleton displacement, u.In order to do so, we should assume that the eight unknowns could be lumped as (p -1 , p -2 , Ω 1 , Ω 2 ) according to the assumption made earlier in Equation ( 12).The reason for this is because it is impossible to analyze, e.g. the effect of p o on-alone without involving p w and p g , and vice versa.By making this assumption, we can now consider a representative volume of the fissured porous medium subjected to the stress conditions as shown in Figures 1 to 3.
Figure 1 Stress state Case I corresponds to an external hydrostatic pressure of d σM , an internal matrix pressure of dp -1 and an internal fissure pressure of dp -2 (adapted from Khalili and Valliappan, 1996).
Figure 2 Stress state Case II corresponds to an equal external and internal matrix pressure of dp -1 (adapted from Khalili and Valliappan, 1996).

Figure 3
Stress state Case III an equal external and internal fissure pressure of dp -1 and a pore pressure of 0 (adapted from Khalili and Valliappan, 1996).Applying Betty's reciprocity theorem to Figures 1 and 2, Figures 1 and 3, and Figures 2 and 3, it can be shown that (Khalili and Valliappan, 1996): According to Equations ( 35) and ( 36), the pore volume change per unit of bulk volume in a fissured porous medium is due to three components on the RHS.The first term is the volumetric change due to a change in the overall bulk volume of the solid skeleton.It depends on the overall bulk compressibility of the material, which is dictated by the multiplying coefficient (modified Biot's coefficients) in front of the volumetric strain.In the case where the material is incompressible, and are equal to unity, otherwise, they are less than one.The second term is the volumetric change due to a change in the fluid pressure occupying the pores or fissures.It depends on the matrix/grain compressibility of the solid constituent.The third term is the volumetric change due to the pressure difference between the pores and fissures.It is responsible for the internal deformation of the material.It should be noted that even under the situation where the matrix is incompressible, there is still a pore volume change due to the first and third term.
The coupled formulation is obtained by substituting Equations ( 35) and ( 36) into (30).Assuming that in most practical problems v s → 0, then: (37) the following compact form of the mass balance equations for each fluid phase is finally obtained: For the water phase in the pore matrix: (38) From Equations ( 38) and onwards, the four index subscripts which appear in the parameter indicate the relationship of each phase in either the matrix or the fractured medium to the other phases, i.e. λ w1o2 should read lambda (λ) relating water in matrix (w 1 ) to oil in fracture (o 2 ).In Equation ( 38), m is a hydrostatic vector defined as: (45) By the same argument, the water equation in the fractures can be obtained as follows: (46) where: In order not to bombard the reader with equations, we give only the expression for water equations here while the complete coupled equations can be found in the Appendix.Before proceeding further, it is important to make a few remarks at this stage.

Remark 1
The governing equations derived here are significantly different from those obtained either by Lewis and Ghafouri (1997) or Bai et al. (1998).In their derivation, the cross coupling terms λ w1w2 , λ w1o2 , λ w1g2 , λ w1w2 , etc. are zero.This implies that the multiphase flow fields in their equations are independent of each other in the matrix and fractures, while in the present formulation these are coupled.These cross coupling terms are present because the internal deformation of the skeleton is a function of the pressure differential between the matrix and the fissures.

Remark 2
The meaning of these cross coupling terms and their influence in Equations ( 38) and ( 46) can be explained qualitatively by looking at the sign in front of these coefficients.For example, the quantity (volume of water per unit of bulk volume per unit of time): (52) in the matrix equation (see λ w1w1 ) implies that the liquid water is attempting to drain from the matrix due to its saturation derivatives S'' w1 .This drainage is possible due to the matrix compressibility , and the internal deformation χ/K, due to the pressure increment ∂p w1 .On the other hand, the quantity: (53) in the matrix equation (see λ w1w2 ) is trying to prevent this drainage from happening in the matrix due to the saturation derivatives S'' w 2 in the fractures.This interference arises because the fractures are trying to prevent the matrix from deforming independently.

Remark 3
The definition of the relative flow vector of each phase, defined by Equation ( 27) is essential in dictating the final form of the coupled multiphase formulation.For example, in Bai et al. (1998), they define the relative flow vector as: (54) Note the difference in the term v sα in Equation ( 54) as compared to Equation ( 27).  , , , , } 1 1 1 0 0 0 What v sα actually implies is that: (55) i.e. two relative flow vector are defined, one with respect to the porous block and the other one with respect to the fissured block.By imposing the stress equilibrium condition: (56) and assuming that the strains are additively decomposable, i.e.: (57) one will find, after following the outlined derivation procedure, the following multiplying coefficient in front of the volumetric strain term: (58) in the matrix Equation ( 38) in which C 1 and C 2 are the elastic tangent moduli of the matrix and fractured block, respectively.Notice that the quantity: defines another compressibility coefficients.It is difficult to explain the physical meaning of Equation ( 59) since already measures the compressibility of the porous block.
The same argument also holds for the fissured block.
Because of this confusion, we should use the relative flow vector of each phase as defined by Equation ( 27), which defines the relative velocity with respect to the solid skeleton, in a macroscopic sense.

Remark 4
The fissured block is modeled as blocks of orthogonal fractures, e.g.parallel plates and as a sugar cube model.The geometric parameter, γ can be written for quasi steady-state condition as (Warren and Root, 1963): (60) and: (61) n = 1,2,3 is the number of normal sets of fractures and d 1 , d 2 , d 3 are the fracture intervals in each direction.

INITIAL AND BOUNDARY CONDITIONS
Normally, the initial conditions of a reservoir system can be defined by specifying the initial distribution of fluid pressure within the reservoir and/or its saturation, depending on the unknowns used in the formulation and the updating procedure.Therefore, the initial conditions for a threedimensional reservoir system can be defined as follows: (62) where p 0 ςα is the fluid pressure at position vector x at time zero and S 0 ςα is the corresponding saturation.For the equilibrium equations, the initial conditions can either be a prescribed initial displacement or a prescribed initial stress, i.e.: (64) For the sake of simplicity in this paper, we should assume that u 0 = σ σ 0 = 0.In the case where the initial stress is nonzero, the above procedure can also be applied, but now the calculated Cauchy stress has to be added to the initial stress value to obtain the in situ stress.This additivity of the stresses is reasonable in this case because we have assumed an elastic constitutive equation for the skeleton and thus the theory of superposition can be applied.
Generally, the flow boundary condition may be written based on Darcy's equation where the flow rate q ςα of the ς-phase, which can be either flow in or out of the system, must satisfy the following condition: (65) In practice, two different types of flow boundary conditions are usually applied, firstly by prescribing the flow rate: (66) or secondly, by specifying the value of pressure at the boundary as: (67) The boundary segments ∂Ω qα and ∂Ω pα have to satisfy the conditions such that: (68) For the equilibrium equations, one requires the displacement boundary conditions, i.e.: (69) and the traction boundary conditions.For the sake of simplicity, the traction boundary condition is set equal to zero, i.e.: (70) Again, the prescribed boundary segments ∂Ω u and ∂Ω t have to satisfy the conditions such that: (71)

NUMERICAL DISCRETISATION
A Galerkin-based Finite Element Method was applied to the governing equations obtained in the Section 2, where the displacements, u, and the fluid pressures, p ςα , in the two overlapping continua, namely the matrix and the fractured network, are the primary unknowns.By first recasting the governing equations in a weak form and approximating the unknown vector as: where , the following equation is obtained: (73) Applying the finite difference implicit θ-time stepping scheme to Equation ( 73) and letting: (74) and: (75) the temporally discretised form of Equation ( 73) will ultimately be given as: (76) where n is the iteration counter.The details of the Finite Element Method will not be elaborated here but can be found elsewhere, e.g.Lewis and Schrefler (1998), Pao et al. (1999), Pao (2000), Hughes (2000).All the coefficients matrices appearing in A and B are nonlinear and are dependent on the values of the sought unknowns.Therefore, iterative procedures are performed within each time step to obtain the final solution.In fact, Equation ( 76) has been implemented in the code CORES (COupled REservoir Simulator) developed at the University of Wales Swansea (Pao, 1998;Pao et al., 1999;Pao, 2000).In order to test the validity of the formulation, the alpha test suite for the implementation now follows.

NUMERICAL EXAMPLE
The selected numerical example is a pseudo threedimensional, field scale case, which was a part of the sixth SPE comparative solution project for double porosity simulators (Firoozabadi and Thomas, 1990).This example was chosen not only because the information required for the simulator is complete, but also a comparison of the solution technique with other simulators is also possible.A linear section of reservoir was modeled (Fig. 4).Vertically, five layers each having a height of 50 ft [15.2 m] was used.Horizontally, the reservoir was divided into ten 200 ft [61 m] grid blocks.A uniform thickness of 1000 ft [305 m] was used in the y-direction.The layer description for the cross-section is given in Table 1.Whereas both the shape factors and fracture permeabilities may be determined in terms of the matrix block size, the corresponding values given in Thomas et al. (1983), and presented in Table 1, were used directly.The oil production rate was calculated as follows: (77) and 506 where ∆p o2 is the well drawdown in psi and I p is the productivity index, defined as: (78) in which k ro2 is the oil relative permeability and µ o2 (cp) and B o2 (ResB/STB) are the viscosity and formation volume factor of the produced oil, respectively, at well bottom hole pressure.The coefficient J has dimension of ResB-cp/daypsi, or Darcy-ft, and the corresponding values for each layer are given in Table 1.
The water and gas rates were obtained in terms of the mobility ratio values as follows:  Saturation versus relative permeabilities in the oil-gas system for the matrix continuum.
Figure 6 Saturation versus relative permeabilities in the oil-water system for the matrix continuum.
Capillary pressure versus saturation in the oil-gas system for the matrix continuum.
Figure 8 Capillary pressure versus saturation in the oil-water system for the matrix continuum.
where M wo2 and M go2 are the mobility ratios of water to oil and mobility ratio of gas to oil, respectively, and are calculated as follows: (81) (82) The production rates are calculated using data evaluated at the beginning of each time step, i.e. an explicit loading is assumed.Depletion runs were carried out to a maximum of 10 years, or whenever production declined to less than 1 STB/D [0.16 STM/D].The production well has a maximum rate of 500 STB/D [80 STM/D] and was limited by a maximum drawdown of 100 psi [6894.757KPa].This well was located at the far right column and perforated only in the bottom layer.Basic pressure-volume-temperature (PVT) data for the example are given in Table 2 (Firoozabadi and Thomas, 1990).Also, Figures 5 to 8 illustrate the saturationrelative permeability and capillary pressure-saturation curves for gas-oil and water-oil systems, respectively.
µ µ These curves apply to the matrix continuum only.Because of the limitation associated with the formulation, a pseudocapillary pressure curve, having a negligible effect on the depletion behavior, was assumed.The variation of gas/oil interfacial tension (IFT), σ, with pressure, as given in Table 2 has also been incorporated into the problem.For this purpose, the gas/oil capillary pressure is directly related to the IFT and therefore the gas/oil capillary pressure should be adjusted according to the ratio of IFT at reservoir pressure, divided by the value of the IFT at which the capillary pressure-saturation curve is specified, i.e. bubble point pressure, or in a mathematical form: (83) where the term p cgI corresponds to input capillary pressure values calculated using a surface tension σ 1 .Instead of an adaptive time step value, an equally spaced time step sizes of 0.1 year were used, and no attempt was made to optimize the time step size.The domain was assumed to be sealed and no flux boundary condition was prescribed at the surfaces.The reservoir was assumed to deform vertically and all displacements perpendicular to the lateral surfaces were fixed at zero.Also, the base of the reservoir was assumed to have no movement in all directions.The Young's modulus of elasticity for the solid skeleton was estimated in terms of the compressibility values of the porous skeleton using the relationship given by equation K = E/(3(1 -2v)) and by assuming a constant Poisson's ratio v = 0.2.Due to the lack of data on matrix bulk modulus, K m and pore fissured bulk modulus, K p , it was assumed that and such that the equivalent Biot-Willis parameter is equal to 0.8.Based on this, the values of K m and K p can be calculated using Equation (19).
Figure 9 shows the plot of oil production versus elapsed time.The results obtained by the present study show a significantly delayed reduction of the production rate when compared to the coupled model of Lewis and Ghafouri (1997).A negligibly small oscillation of the production rate was observed when the production rate was still high.The oscillation soon damped out after a few time steps.It is important to mention here that most of the previous uncoupled models predict a sharp decline of the production rate after nearly 2-6 years, the present coupled model predicts a much longer period with a smooth decline in total oil throughput, which finally ends after 9.2 years.
The same trend may be observed for the gas/oil ratio, GOR, in Figure 10.The results obtained by the SPE participants suggests a sharp increase in GOR after a few years due to the decline in the reservoir pressure, the present model, where coupling effects are also considered, this was postponed until the ninth year, see Figure 10.Note that the rate of GOR increases abruptly after nine years.We believe that the reason why the solution procedure halted at 9.3 years was mainly due to the extremely high GOR ratio calculated using Equation ( 80) during that period, thus causing the solution to diverge.Again, oscillations are observed in the producing GOR during the initial period due to the fluctuation in the oil production rate.
In order to investigate the influence of the reservoir pressure drawdown, the average pressure values in gridblock (5,1,1), are compared in Figure 11.The rate of pressure decline is significantly slower when the coupling effect is incorporated (Fig. 11), whereas all the SPE uncoupled models predict a similar trend with a rapidly decreasing pressure.This confirms the considerable impact of rock deformation on maintaining the initial reservoir pressure for a longer time, which can improve long term economical α ˜.   productivity of a reservoir.This effect, which is ignored in conventional uncoupled models, is particularly of great significance when the oil-bearing strata are quite deformable.An interesting comparison between the coupled models in Figure 11 shows that the present model predicted slower pressure dissipation in the reservoir.This observation is in accordance with the results reported in the literature (Khaled et al., 1984;Bai et al., 1993) where the matrix block pressure tends to dissipate much slower due to the presence of fissures.Another reason that leads to this result is the fact that in the present case, both the matrix and fissured blocks are assumed to be compressible (dictated by the value of K m and K p ), while in (Lewis and Ghafouri, 1997), these values were assumed to be zero due to their formulation.In addition to that, the present formulation which include the cross coupling terms can be another factor that slows down the pressure dissipation.Isobaric contours in the cross-section of the reservoir are shown in Figures 12-14 for 1, 5 and 9 years, respectively.The pressure declines towards the producing well with a relatively high-pressure gradient near the well bore.The pressure gradient in the region adjacent to the wellbore is usually very high compared to the far field.This is firstly due to a very high flow rate per unit area in this region.Moreover, in many producing wells the value of absolute permeability of the formation close to the well bore is greater than the values further back in the reservoir.This causes even more pressure drawdown around the well bore, the so-called skin effect.Due to the combined effect of the two foregoing factors, usually a very rapid drawdown occurs around the well.However, this effect is not predicted in the present simulation.The reason is that the size of elements used in this problem is larger than that of the region affected by the phenomenon, which is usually of the order of a few meters.This implies that the pressure value around the well   obtained in this study may be significantly different from the real values.More accurate results may be obtained by refining the mesh in this region.However, the pressure values in the other regions are by no means dramatically affected.
In all the solutions obtained by the participants, including Lewis and Ghafouri (1997), no solution for the reservoir displacements was shown.Figure 15 illustrates the vertical reservoir subsidence profile at the top surface for different time intervals.The reservoir subsides almost 1.0 ft over a period of 9.2 years, giving a subsidence rate of approximately 0.1 ft/year.This seemingly small value can however do considerable damage to the oil platform, well casing and pipeline, e.g.some old fields in the United States, i.e.Buena Vista (0.27 m), Fruitvale (0.04 m) and Santa Fe Spring (0.66 m). Figure 16 illustrates the displacement vs. time curve at the gridblock (5,1,1).Note that the displacement is almost linearly proportional to the oil rate and the pressure drawdown.

CONCLUSIONS
The mathematical formulation of a three-phase threedimensional fluid flow and rock deformation in fractured reservoir was introduced in this paper.The present formulation, consisting of both the equilibrium and multiphase mass conservation equations, accounts for the significant influence of coupling between the fluid flow and solid deformation.The final set of equations, which contain the additional cross coupling terms as compared to similar existing models, accounts for the displacements compatibility condition between the porous rock and the fissured block.A field scale example was employed as an "alpha" case to test the soundness of the current formulation and the robustness of the numerical scheme.The results illustrate a significantly different behavior for the case of a reservoir where the impact of coupling was also considered.
Figure 4Finite element mesh used for the SPE example (not in scale).
Figure 5 Figure 9Oil rate versus time for the SPE example.

Figure 10 Gas
Figure 10Gas/oil ratio versus time for the SPE example.

Figure 12
Figure 12 Pressure distribution (psi) contours in the cross-section of the reservoir at 1 year (1 psi = 6894.8Pa).

Figure 13
Figure 13 Pressure distribution (psi) contours in the cross-section of the reservoir at 5 years (1 psi = 6894.8Pa).

Figure 14
Figure 14 Pressure distribution (psi) contours in the cross-section of the reservoir at 9 years (1 psi = 6894.8Pa).