Fully Coupled Flow-Deformation Model for a Naturally Fractured Gas Reservoir

The Lacq gas field, a large scale fractured reservoir subjected to a significant depletion in its production life, has been investigated. A fully coupled model two phase fluid flow in deformable fractured porous media has been utilised in the investigation. Spatial discretisation of the governing equations has been achieved using the standard Galerkin method while the finite difference technique is employed for the discretisation of the time domain. A three dimensional finite element model representing a quarter of the reservoir and surrounding rock layers has been created and subjected to the actual course of field production from 1957 to 1989. A good agreement has been found between simulation results and field data on reservoir pressure evolution, maximum ground settlement and subsidence profile along the northwest-southeast direction.


INTRODUCTION
Fluid flow and deformation are naturally coupled phenomena in studying hydrocarbon reservoirs.As a result of this coupled behaviour, subsidence and seismic activities are more common ground's responses to fluid extraction/injection from/into underground reservoirs.For example, in the production life of Ekofisk field in the North Sea, seabed subsidence has been experienced over an area of 50 km2 [5].Accurate modelling of depletion induced subsidence requires an understanding of the mechanical behaviour of the reservoir and surrounding rock in response to the depletion.Triggering seismic events in case of fluid injection has been proven to be a straightforward phenomenon as shown by observations and tests performed in Colorado [6].Also, in some cases, production has been found to initiate seismic activity [16,18,21,23].In Lacq gas field (France), while seismic activities first time felt in 1969, 12 years after start of production, the surface subsidence has been observed before first earthquake during the levelling in 1967 and repeated in 1979 and 1989.However, the surface subsidence records show a linear relationship between reservoir depletion and ground settlement in spite of seismic activities, showing an elastic behaviour of ground to reservoir depletion.Segall et al [23] suggested that the induced earthquakes are caused by poroelastic stressing associated with fluid withdrawal.In their model, the coupling between fluid flow and deformation realised using effective stress concept while the reservoir considered free of constraints on its boundaries i.e. the change in total stress assumed to be zero (implicit model, [22]), and they directly used the field pressure data to compute the deformations and stresses in the vicinity of the reservoir.
In Lacq field, a highly fractured reservoir is located in an anticline structure of the French Pyrenean foreland.The reservoir consists of Jurassic dolomite with levels of low porosities (4-6%) alternating with levels of very low porosities (less than 1%).The matrix permeability is extremely low (0.001 mD) while the well tests give permeability values from 5 to 200mD.This gives evidence of "double porosity" behaviour of the reservoir.The migration of gas is the result of fracturing which is responsible for high permeability of the rock [20].
In this paper, a fully coupled flow-deformation three-dimensional finite element model has been utilised to analyse the behaviour of Lacq gas field during its production life till 1989 corresponding to a significant pressure drop of about 60 MPa.The model is based on double porosity concept with two phase fluid flow (gas and water).The coupling between flow and deformation has been realised through effective stress concept [10].Good agreement has been found between model data and field data.

LACQ GAS FIELD CHARACTERISATION
Lacq gas field is a large west-east oriented reservoir located in an anticline structure of the French Pyrenean foreland.The depth of the reservoir changes from 3.2 km to about 5 km due to the anticline structure with an average depth of 4 km and effective thickness of 250m.The structural traps are formed by impermeable marly levels which bound the reservoir.The initial pressure of Lacq Profond gas field was 66 MPa when production began in 1957 and reduced to 9 MPa in 1989, a significant drop which caused more than 900 seismic events with magnitude larger than 1 between 1974 and 1989 [3].
The reservoir consists of Jurassic dolomite with levels of low porosities (4-6%) alternating with levels of very low porosities (less than 1%).The overall permeability of blocks is less than 2%.The matrix permeability is extremely low (0.001 mD) while the well tests give permeability values from 5 to 200mD.The migration of gas is the result of fracturing which is responsible for high permeability of the rock.This gives evidence of "double porosity" behaviour of the reservoir.The reservoir responded to the depletion as a homogenous matrix.No compartments with different relative pressures have been observed during gas extraction [4,17,20].
Geological formation of the field consists of 5 major layers: below levels (BL), reservoir, marls, reef and upper levels (UL).The mechanical properties of these layers are given in table 1 [2,3].The appearing of fractures in reservoir rock will affect the deformability modulus, so the overall deformability modulus calculated based on given intact value (54 GPa) and relationship proposed by Hoek and Diederichs [7] assuming a very blocky structure.With an OGIP of 2.6x10 11 Sm 3 , production from Lacq field began in 1957 and reached to the rate of 2.1x10 7 Sm 3 /d in 1961 and lasted for 24 years until the first decrease in plateau rate occurred in 1985 followed by another decrease in 1986 reaching the plateau rate of 8x10 6 Sm 3 /d.The gas was initially highly over-pressured (66 MPa), and production has caused the gas pressure to decline dramatically to 9 MPa in 1989.The temporal evolution of gas pressure has been shown in figure 1.
Figure 1.Gas pressure evolution  Subsidence associated with gas extraction has been observed and measured along a northwest-southeast trending line in four different times.This data provide valuable information on the nature of coupling between fluid flow and deformation in Lacq field.First levelling has been performed in 1887, far before any production, followed by three other levelling in 1967, 1979 and 1989.The measured subsidence has been shown in figure 2 with a maximum subsidence located in centre of reservoir.Plotting the maximum observed surface subsidence against the reservoir depletion (figure 3) reveals an almost linear relationship between deformation and pressure drop proving the linear response of the ground to the applied loads due to gas extraction [23].A shallow oil field at 640-700 m depth was discovered in 1949, and by 1986 a total volume of 4x10 6 m 3 oil was produced.The oil pore pressure has changed between 6.5 MPa and 5.9 MPa which is not comparable to the pressure drop in gas reservoir, so the surface subsidence is mainly due to gas extraction from lower reservoir.

COUPLED FORMULATION
The coupled formulation consists of two separate, yet overlapping models: the deformation model, and the flow model.The deformation model is based on the theory of elasticity, the effective stress principle and the equations of equilibrium.The flow model is based on the theory of double porosity [1,24].Two interacting fluid flow regions are identified: one representing the flow in the porous blocks, and the other representing the flow in the fracture network.The two flow regions are coupled through a leakage term controlling the transfer of fluids from the porous blocks into the fracture network, and vice versa.The rate of fluid transfer at any point is assumed to depend on the fluid pressure difference in the porous blocks and the fracture network, the phase permeability of the porous blocks, phase viscosity and the fracture spacing.
The coupling between the fluid flow and deformation models is established through the effective stress principle.In its most general form, the effective stress for a fractured porous medium, saturated by two fluid phases, is expressed as [10]: in which σ is the total stress, β is the incremental effective stress parameter, p is the fluid pressure, and I is the second order identity tensor.Incremental effective stress parameters are quantified as following [22]: is the compressibility of porous block and c is overall compressibility of fractured porous medium.Coefficient β represents the contribution of pores to the overall deformability of medium [11], while ψ represents the contribution of wetting fluid into deformation of pore α.
represents the compressibility of solid grains.The effective stress parameters are restricted to the following constraint: The α 1 is the well-known Biot coefficient.Segall et al. [23] proposed that Biot coefficient for Lacq field is equal to 0.25 using the experimental data reported by Laurent et al. [12] for different limestones.We modified the quantity of α for fractured medium, assuming that the α 0.25 is still valid for porous blocks: Combining equation of equilibrium and the effective stress equation yields the differential equation describing the deformation field of an elastic fractured porous medium saturated with two immiscible fluids The flow model is developed by combining the equation of linear momentum balance for the fluid phases (Darcy's law) with the mass balance equation of the fluids (neglecting the inertial and viscous effects): sign α for α p for α f is the formation volume factor, k and µ are permeability and dynamic viscosity of phase , c is the coefficient of fluid compressibility and d V /V represents the change in the volume of the fluid constituent β in pore space, α, over the volume of the fractured porous medium.Γ is the leakage term controlling the exchange rate of constituent (β) between the pores and fractures.
The leakage term has an important role in simulation of fractured reservoirs.The fluid transfer between pores and fractures is produced by the leakage term.Basically the leakage term defines the early, intermediate and late time response of a fractured reservoir.In two phase fluid flow it controls the imbibitions of matrix blocks and counter current of flow from matrix blocks to fractures.The leakage term can be explicitly evaluated by employing the Darcy law between physical matrix block and fractures, or implicitly by introducing the shape factor σ that takes into account the shape and size of virtual matrix blocks.However, many simulators use the latter to evaluate the leakage term in the form of: Both leakage and pressure difference in pore and fracture are time-dependent and also depend on the boundary conditions.However the shape factor converges to a constant value in the later time.The dimension of shape factor parameter is the inverse of a squared length i.e. shape factor parameter is a constant divided by squared length.So far several workers have used different approaches to evaluate shape factor constant analytically and numerically.For a cubic matrix block (three sets of fractures), the shape factor constant varies from 12 for Kazemi et al [9] to 60 for Warren and Root [24].Among them, the values calculated by Lim and Aziz [13] and Quintard and Whitaker [19], which are the analytical solutions to the diffusion equation for different boundary conditions, seem to be more appropriated.Here the value calculated by Lim and Aziz [13] has been used, (7) In which, , and are matrix block's dimensions in x, y and z directions.The fully coupled equations are obtained by combining the governing equations for the deformation and flow models: The coefficients Λ , and Λ , will be determined from the capillary pressuresaturation characteristic curves for pores and fractures: In which, n and n are the pores and fractures porosity, S and S are the pores and fractures saturation, and, P and P are the capillary pressure in pores and fractures respectively.

NUMERICAL MODEL AND RESULTS
Finite element method has been invoked to solve the governing equations.Spatial discretisation has been performed using the Galerkin method.Displacements and phase pressures are considered as primary variables and approximated in space by means of the shape functions and their nodal values, and , The shape function vector or matrix consisting from nodal shape functions N are defined as, A weak form of the governing equations is needed for their transformation into a set of algebraic equations in space.For this purpose, using Galerkin method which is based on minimizing the weighted average of residual error of the solution over the entire element domain, the governing equations are multiplied by vector T as the weighting factors, and then integrated over the whole element domain, Ω with boundary Γ, in which, is a differential operator corresponding to the definition of engineering strains, is the gradient vector, is the elastic stiffness matrix, is the identity vector, and is the permeability tensor.By introducing, , T and , the derivatives of the shape function for the displacements, coupling and phase pressure terms, and expanding the equations 11 and 12 by application of Green theorem, the final form of the governing equations may be expressed as, 13 ∑ , , , sign α 14 in which, is the element stiffness matrix, is the coupling matrix, is flow matrix corresponding to the mobility of phase αβ, , and are mass matrices corresponding to cross coupling between phases αβ and ζξ, and leakage terms respectively, is the vector of nodal forces, is the vector of nodal flux of phase αβ and is the gravity vector, Temporal discretisation is performed using the standard finite difference technique such that for an arbitrary function f, over a time interval Δt one can write [22], with the parameter θ (0 θ 1) controlling the type of approximation used for time discretisation.For different values of θ, we obtain different well-known numerical integration schemes such as forward (θ 0) and backward (θ 1) schemes.The Galerkin scheme θ 2/3 provides unconditional stability and first order of accuracy.Experience has shown that Galerkin scheme provides better stability than the Crank-Nicolson scheme θ 0.5 in spite of the lower accuracy.With Crank-Nicolson scheme the profile of the solution displays some oscillations, particularly in the vicinity of the perturbation.Having discretised spatially and temporal, the final set of linear algebraic equations are: Equations ( 17) and ( 18) provide a set of linear algebraic equations and can be written in matrix form of , in which is the element's general stiffness matrix, is the unknowns' vector and is the generalized force vector.Since the coefficients in stiffness matrix are dependent on the unknowns, an iterative procedure is employed within each time step.At each iteration level, the nonlinear coefficients are updated using the values of unknowns calculated in the previous iteration.The generalized force vector is updated in each increment.The eight-node prism element has been used to interpolate all unknowns.Two-point integration scheme (a total of eight points) has been used.The finite element code has been developed as part of this research.
A uniform flat reservoir with effective dimensions of 12km x 6km x 250m was assumed.A three dimensional model of a quarter of the reservoir and surrounding rocks has been created.The size of the model is 20km x 20km x 6.5km.The reservoir and four rock layers, listed in table 1, has been shown in figure 4. The matrix and fracture porosity were set to 1.8% and 0.5% respectively.The intrinsic permeability of matrix is considered 0.001 mD (1x10 -18 m 2 ) while the fracture's is 200mD (2x10 -13 m 2 ).The production simulated by applying uniform extraction rate all over the reservoir, starting from zero in 1957 and increasing linearly to 2.1x10 7 Sm 3 in 1961.The production rate remains unchanged till 1985, and then reduces to 8x10 6 Sm 3 until 1989, end of simulation.The simulated reservoir pressure evolution has been shown in figure 5 which shows good agreement with field data.data The maximum subsidence which is occurring in the middle of reservoir versus gas pressure depletion has been shown in figure 6 which shows a linear behaviour as mentioned before.The surface subsidence from this simulation along the northwest-southeast direction corresponding to the times that the levelling was performed (1967, 1979 and 1989) has been shown in figure 7 showing good agreement with the field data.In figure 8, the surface subsidence, along the northwest-southeast direction, from this simulation in 1989 was compared with that of Segall et al. [23] and it was shown that present model predict better subsidence.

CONCLUSION
A fully coupled two phase fluid flow in fractured porous media has been utilised to investigate a gas filed case.Lacq gas field is a huge fractured reservoir subjected to a significant depletion in its production life during 1957 to 1989.As a result of coupling nature of flow and deformation in reservoir engineering, surface subsidence and seismic activities have been observed in Lacq field.A three dimensional finite element model of quarter of the reservoir and surrounding rock layers has been created and subjected to the actual course of production.A good agreement has been found between simulation results and field data on reservoir pressure evolution, maximum ground settlement and subsidence profile along the northwest-southeast direction.

Figure 4 .Figure 5 .
Figure 4. Three dimensional model of Lacq gas reservoir and surrounding rock layers

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Maximum surface subsidence from this simulation and field data versus reservoir depletion

Table 1 .
Mechanical properties of geological layers of Lacq gas field