- Email: [email protected]
Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations
Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 81 (2011) 2270–2281
Original article
Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations V.V. Mourzenko a,∗ , I.I. Bogdanov a,b,1 , J.-F. Thovert a , P.M. Adler b a
Institut PPRIME, SP2MI, BP 30179, 86962 Futuroscope Chasseneuil cedex, France b UPMC - Sisyphe, 4 place Jussieu, 75252 Paris cedex 05, France
Received 15 November 2009; received in revised form 8 June 2010; accepted 8 December 2010 Available online 21 December 2010
Abstract A general three-dimensional numerical model for single phase, slightly compressible flow through fractured porous media is described. It is based on a discrete fracture representation. Three sets of applications are presented. In the first one, pressure drawdown well tests in closed oil reservoirs are simulated for complex model situations where the well intercepts a random fracture network with various fracture densities and conductivities. Then, the hydrodynamic response of a fractured aquifer is investigated by simulating on the field-scale single- or two-well pumping tests in the Poitiers Hydrological Experimental Site. Finally, a complete field-scale simulation of the production history in an oil reservoir with multiple wells is presented. © 2010 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Porous medium; Fracture; Reservoir; Well flow
1. Introduction This work addresses single phase, slightly compressible flow through fractured porous media. A very general, three-dimensional numerical model based on a discrete fracture representation is proposed, together with applications involving a single or multiple interacting wells. Due to the specific transport properties of fractures, the flow through a naturally fractured porous medium differs drastically from that in a conventional porous medium. The key feature is that the porous matrix provides the main storage for the fluids while transport takes place mainly through the fracture system. Furthermore, matrix/fracture flow interactions govern many of the medium transport properties. Because of the complexity and diversity of most natural fracture systems, the determination of fractured porous media transport properties remains an open issue of great practical importance. For instance, the present numerical tool can be applied for the interpretation of well test data, in order to quantify the characteristics of a reservoir, or conversely, to optimise the design of a producing well, given the reservoir characteristics.
∗
Corresponding author. Tel.: +33 5 49 49 67 72. E-mail addresses: [email protected] (V.V. Mourzenko), [email protected] (I.I. Bogdanov), [email protected] (J.-F. Thovert), [email protected] (P.M. Adler). 1 Present address: CHLOE - UPPA Avenue de l’Université, Bat. B1, 3ème étage droite, BP 1155, 64013 Pau cedex, France. 0378-4754/$36.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2010.12.014
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2271
Most existing models of transient compressible flow through fractured porous media are based on the concept of double-porosity media [2,18,8,14,10]. A discrete approach allows a more detailed description of the matrix/fracture flow interaction accounting for the flow behavior in individual elements of the medium, but it has been applied so far only to regular system of fractures [16,20,5,21,15]. We presented in [4] an extension to unsteady compressible flow of the numerical models of [12,3], which make use of a full three-dimensional description of random discrete fracture systems. Its main features are recalled here, as well as results of single-well tests simulations in medium-sized systems (hundreds of fractures), focused on the well response. Then, new applications on a much larger scale (thousands of fractures) are presented, involving many wells which can be simultaneously producing or simply observing. Attention is extended to the transient complex three-dimensional flow field in the fractured aquifer or reservoir. For instance, the possible disturbing effects of supposedly passive wells are demonstrated. 2. Governing equations The fractured porous medium can be represented as an arbitrary, generally random, set of fractures embedded into a solid porous matrix. At the large scale, each fracture can be viewed as a vanishingly thin layer, with singular transport properties. For instance, a steep pressure gradient may exist across a poorly conducting fracture, which results in an apparent pressure jump at the macroscopic scale. In addition, the typical time for pressure variations in transient flows is quite different in the matrix and in the fractures. The discrete fracture model used in the present work explicitly accounts for these features. Therefore, it allows detailed investigations of transitory flow regimes, and of the influence of the well intersections with the fracture network. 2.1. Flow equations Consider a porous matrix with a bulk permeability Km [L2 ] which may vary with space. Darcy law for the local seepage velocity v and the mass conservation for slightly compressible flow can be written in the matrix as v=−
Km ∇P μ
(a),
m Cm
∂P + ∇ · v = δw Jw ∂t
(b)
(1)
where μ is the viscosity, P is the pressure, m and Cm are the matrix porosity and total compressibility. Jw [L2 T−1 ] represents the exchanges with the well, at a location given by the Dirac function δw [L−2 ]. On this scale of description, the well appears as a line without thickness. If m Cm is constant, a matrix pressure diffusivity Dm can be defined and (1) can be written as a difusion equation for the pressure ∂P δ w Jw − ∇ · (Dm ∇P) = , ∂t m C m
Dm =
Km μm Cm
(2)
We assume that the hydraulic properties of a fracture can be characterized by two effective coefficients, namely a conductivity σ[L3 ] and a cross resistance ω[L−1 ] (see [3]), which relate the in-plane the flow rate j to the surface pressure gradient ∇ s P and the seepage velocity v⊥ of the net flow crossing the fracture to the pressure drop [P] across it by σ j = − ∇s P μ
(a),
v⊥ = −
1 [P] μω
(b).
(3)
For illustration, the fracture can be viewed as a plane channel of aperture b, filled with a porous material with permeability Kf , porosity f and total compressibility Cf . Then, σ and ω are given by σ = bKf and ω = b/Kf . A non zero resistance ω may exist even for totally open channels, if their walls are partially clogged by a deposited chemical [13]. Note that equations similar to (3) for the fracture hydraulic properties have been obtained by [7] by an upscaling procedure. The continuity equation for the flow through a fracture reads as bf Cf
∂P + ∇s · j = (v− − v+ ) · n + bδw Jw ∂t
(4)
2272
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
where v± is the seepage velocity in the matrix on either sides of the fracture, given by (1). Again, if bf Cf is constant in a fracture, a pressure diffusivity Df can be defined, and Eq. (4) can be written as ∂P (v− − v+ ) · n δ w Jw + , − ∇s · (Df ∇s P) = bf Cf f Cf ∂t
Df =
σ bμf Cf
(5)
The assumptions of uniform m Cm and bf Cf used to write down the diffusion equations (2) and (5) are by no means a requirement for the numerical model, since the discretized equations are really based on the general equations (1) and (4). The problem is governed by four dimensionless numbers relative to the fractures (b , σ , ω , η), in addition to the well properties, b =
b , R
σ =
σ , RKm
ω =
ωKm , R
η=
f Cf b m Cm R
(6)
In the present work, we only consider very conducting fractures and neglect pressure variations along the well. Unless otherwise stated, all the results are given for fractures without cross-resistance to flow (ω → 0), with large aspect ratio (b = η = 10−3 ), and a large fracture/matrix conductivity contrast (σ 1). The unit scale R is a characteristic lateral extent of the fractures. Other choices are possible. In particular, the well radius rw can be used as the characteristic length scale instead of the fracture size R. This is appropriate to describe the very early stage of a well test, when the flow is confined very near to the well [9]. Alternatively, the size of the drainage area A of a closed reservoir can be used as length unit (see Eq. (11)). 2.2. Overall boundary conditions The model can accommodate the usual boundary conditions in reservoir engineering, namely of Dirichlet type, with an imposed far-field pressure P∞ , or of Neumann type, with no-flux conditions at the boundaries of a closed reservoir. For systematic studies of the interactions of a single well with a fractured formation, a simple geometrical setting is used, with a parallelepipedic cell. The well is set along a vertical line throughout the cell, and spatial periodicity is applied in the x-, y-directions. Hence, this situation corresponds to a periodic square array of infinite vertical wells, or to a single well in a closed reservoir. All the simulations reported here correspond to a pressure drawdown test, i.e., to a well producing at constant rate Qw from a field initially at rest. The initial pressure is arbitrary, and can be taken as P0 = 0 without loss of generality. Of course, other settings can be devised to address any specific field case as illustrated in Section 6. 3. Numerical model We consider here both simple deterministic test configurations and more realistic and complex randomly fractured porous media. Whereas specific discretization techniques could be used in the first situation, for instance with local grid refinements, we are really interested in the latter case where geometrical randomness requires great flexibility. The numerical model is a generalization of [12,3] for fracture networks and fractured porous media. More details can be found in [4]. 3.1. Triangulation of fractured media and description of the well geometry The fracture network and the space located in between the fractures are successively triangulated by an advancing front technique. The scale of the discretization is controlled by a typical spacing δm between the mesh points. Threedimensional views of examples of triangulated fractured medium are displayed in Fig. 1. In a subsequent step, the well is superimposed on the three-dimensional unstructured mesh. Note that the initial triangulation is performed without considering the well. An alternative is to take the well into account from the start during the construction of the mesh, so that nodes would be distributed right on the well. However, the 3D triangulation of complex media is a time-consuming step. The present procedure allows to perform it only once for a given medium, and then to vary the well geometry or location at negligible cost.
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2273
Fig. 1. Three-dimensional views of triangulated periodic fractured media. The example in (a) contains Nfr = 32 hexagonal fractures with circumscribed radius R. The cell size is L = 4R, and δm = R/4. The tetrahedral volume elements with their centers in the cubic unit cell − 2R ≤ x, y, z ≤+2R are displayed. The protuding fractures sit astride the boundaries with the neighboring cells. The intersections of the fractures (thick lines) and of the tetrahedral volume elements with the cell boundaries are shown in the sample in (b), with L = 6R, Nfr = 166 and δm = R/4. This mesh contains 17,671 nodes and 111,556 tetrahedra.
The typical size δm of the volume and surface elements has to be small compared to the typical size R of the fractures, for a good rendering of the fractured medium geometry, and large compared to the well radius rw , since the well is considered as a line without thickness in the model.
3.2. Finite volume formulation The rock matrix is represented by tetrahedral volume elements, with given properties m , Cm and Km , and the fractures by triangular surface elements with given properties f , Cf , b, σ and ω. All these parameters can be set independently, on a per-element basis. The pressure is evaluated at the vertices of the surface and volume elements (◦, in Fig. 2c). The well is represented by a series of segments in the volume elements that it crosses. The pressure in the well is evaluated at the mid-points of these segments, and at the intersection of the well with fractures (• and , in Fig. 2c). The equations for the node pressures are obtained by integration of the balance Eqs. (1) and (4) over a control volume around each mesh point. Several types of control volumes exist according to the location of the mesh point, in the bulk of the rock matrix or on a fracture, with or without a well nearby. The simplest case of a node point in the matrix, without well nearby, is depicted in Fig. 2a. The continuity equation (1b) is integrated over a volume , made up of one fourth of all the incident tetrahedral volume elements
∂P ∂P m Cm m Cm dv + n · vds = 0 + ∇ · v dv = ∂t ∂t ∂
(7)
Since the matrix properties m , Cm and Km are taken as uniform over each tetrahedral volume element, the seepage velocity v in each tetrahedron is piecewise constant over ∂, and can be evaluated from (1), where the pressure gradient is deduced from the pressures at the four vertices.
2274
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
Pressure may take different values P± in the matrix on either sides of a fracture. A third value of the pressure Pf is introduced, in the middle plane of the fracture, which is used for Eq. (4). Balance equations are written over the three control volumes depicted in Fig. 2b, ± on either side of the fracture and f in the fracture, ±
m C m
∂P ± dv + ∂t
∂± −F
n · vds = q±
f
bf Cf
∂Pf ds + ∂t
∂f
n · jdl = q− − q+
(8)
The flux j in each triangle intercepting f is evaluated from (3a) where ∇ s P is a linear combination of the pressures at the triangle vertices, and the exchange terms q± are modeled according to (3b). The balance equations (8) yield three linear equations with the same form as for (7). More complex situations occur at the intersection of two or three fractures, but they are handled in a similar way. The well is treated as a one-dimensional highly conductive system superimposed on the medium. Accordingly, it can be viewed as a line source or sink in the fractured porous medium, decomposed into elementary parts corresponding to its intersections with the grid surface and volume elements. The source terms involving Jw in (1) and (4) have to be added to the balance equations for the control volumes (7) and (8). The flux Jw can be related to the well pressure Pw by decomposing the pressure field near the well into the superposition of a global trend, induced for instance by a large scale transverse flow and represented by a constant pressure gradient G in the mesh element, and of a disturbance induced by the presence of the well, modeled according to the analytical solution for a straight well in an unbounded medium. A finite volume formulation can also be implemented to describe the flow along the well, but in the present case, where pressure variations along the well are neglected, it reduces into an integral condition on the pumping rate Qw . 3.3. Solution of the discretized equations Time derivatives are discretized to first order, in a fully implicit formulation. Since the matrix for the set of linear equations to be solved at each time step is not symmetric, due to the exchange terms with the well, a biconjugate gradient algorithm is used. There is no stability restriction for the choice of the time step δt. However, it has to be set small enough in order to sample the transients of the flow that we are interested in, such as the early stages of the flow around a well, or the propagation of a wave pressure along a fracture which can be much faster than in the imbedding matrix. This implies criteria of the form Dδt < h2 where h and D are the relevant length scale and diffusivity. For instance, D = Df and h is the fracture extension when tracking the travel of a pressure wave along a fracture. The wide range of length scales and the contrast between the diffusivities of the fractures and matrix make a careful management of the time step necessary. In practice, different time steps have to be used at various stages of the simulations, depending on the relative contributions of the fractures and matrix to the flow and on its current transitory or quiescent character.
Fig. 2. Control volumes around a point in the matrix (a) and on a fracture (b). Two-dimensional illustration of the introduction of the well in the mesh (c).
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2
10
,Π
P
w,DA
w,DA [C] [D] Pseudo− "Hydraulic steady fracture" flow fracture flow
1
10
2275
[A] Early time fracture flow
1
[B] Transient fracture flow
P
w,DA
→
0
0.46
10
[E] Pseudo−steady global flow
← Πw,DA −1
10
1/25 −−−−−> 1/26 −−−−−−−−−−−−−>
0.78
tʹ
−2
10
−5
10
−3
10
−1
10
1
10
Fig. 3. The well pressure Pw,DA and its logarithmic time derivative w,DA versus the dimensionless time tDA , for a straight well intersecting a series of parallel plane fractures, with spacing H = L, in a closed reservoir with a square drainage area L × L. Data are for rw = 5 10−3 R, σ = 50, b = η = 5 10−4 , and δm = R/8. Curves correspond to Pw,DA and w,DA with (—) or without (−·− · −) well/matrix exchanges, and to the analytical prediction (10) for Pw,DA based on the raw (−−− −) or corrected (··· · ·) value of Df (see [3]). Straight lines indicate fitted slopes 0.78, 0.46 and 1. Horizontal arrows correspond to the values of the plateau of w,DA predicted by (9), with or without well/matrix exchanges.
4. Comparison with analytical solutions A few simple cases have analytical solutions, which have been used in [4] as references to assess the performance of the numerical model, such as the so-called drawdown test for a single vertical well in a uniform porous medium without any fracture, initially at rest, at a given constant pumping rate Jw per unit length of the well. The exact solution reduces after a short initial period to the well-known approximate logarithmic form [9] 4Dm t μJw 100r 2 ln 2 + 0.81 , (9) P(r, t) = ±1% for t > 4πKm r 4Dm Another basic solution is the pseudosteady long-time solution for a finite region with drainage area A [6,17] A Jw t μJw 2.2458 ln 2 + ln P(r, t) = , (10) + m Cm A 4πKm r CA where the shape factor CA depends on the reservoir geometry and on the location of the well. The term independent of time in (10) allows to evaluate the medium properties in the finite flow region. The results (9) and (10) can be straightforwardly transposed to the radial flow in a fracture in an impervious matrix. According to (10), the results can be expressed in terms of the dimensionless time tDA based on the drainage area A = L2 and of the dimensionless well pressure Pw,DA based on the imposed flow rate Jw tDA =
4Dm t , L2
Pw,DA =
4πKm Pw μJw
(11)
The results of the numerical simulations agree with the theoretical predictions (9) and (10) within a few percents for time and spatial steps δtDA ∼ 0.01 and δm /L ∼ 0.1. A well intersecting a single fracture has also given rise to analytical results [6,19] and can constitute the basis for modeling more complex situations. Our simulations, with or without taking into account the well/matrix exchanges, account for the various well-known successive regimes (Fig. 3). At early times, flow is mainly due to fracture compressibility. After a very brief transient (region [A]), a plateau in the curves for w,DA [B] corresponds to the prediction of (9); it is followed by a steeper increase [C], where Pw,DA grows as t0.78 . If the porous matrix were impervious, the pressure would increase linearly
2276
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
with time at this stage, according to (10). Later on, flow in the matrix becomes significant, and another regime is reached [D]. The slope of the curve of w,DA is about 0.46, up to tDA = 1/4, which is of the order of the typical time tDA ∼ 1 (t ∼ H2 /4Dm ) for the pressure wave to travel from the fracture to the middle of the matrix layer. This regime corresponds to a power dependence of Pw on time as t0.46 , as observed with hydraulic fractures, rather than to the halfed logarithmic law of [6], again because of the impervious domain boundaries. Finally, another pseudo steady regime [E] is reached for tDA > 1/4, with a law of the form of (10). Two hydraulic conductivities can be measured from the well pressure curve; the plateau [B] of w,DA provides the fracture conductivity (see (9)), and an extrapolation at t = 0 of the linear increase of Pw in the pseudo steady regime [E] provides an estimate of a global effective permeability Keff for the fractured medium (see (10)). 5. Well through a randomly fractured porous medium 5.1. Fracture network model We consider here randomly fractured media, as modeled by [11]. The fractures are plane polygonal objects, randomly located and oriented in space (see Fig. 1). The density of fractures is quantified by the dimensionless parameter ρ , which is the number of fractures per excluded volume Vex [11,1]. It also quantifies the connectivity of the fracture network, since ρ is the mean number of intersections per fracture. In particular, ρ determines the percolation of the fracture network, with a percolation threshold ρc ≈ 2.3, which has a dramatic influence on the macroscopic permeability of the fractured medium [3]. For very conducting fractures, it is expected that the number ni of intersections of the well with the network is a key parameter for the flow, since it conditions the fluid exchanges. Since the fractures are located according to a Poisson process, their intersections with the well are also Poissonian. Thus, their spacings si obey a negative exponential law f (si ) =
1 −si /Si e , Si
Si =
2 S
(12)
where S is the volumetric surface area of fractures and Si is the mean spacing. For a well of length L, the mean expected number of intersections is Ni = L/Si . 5.2. Examples of pressure drawdown curves As a rule, all the typical flow regimes present in Fig. 3 can be observed in the drawdown curves for wells in randomly fractured media. In particular, the two latest regimes, namely, the “hydraulic fracture” and the pseudo-steady flow, can always be clearly identified. However, the early and intermediate regimes are often blurred or do not show up at all, and additional features can appear in the curves. First, during the period when flow takes place mainly in the fractures, the pressure wave starts propagating in the fractures which cross the well, then spreads out in all the fractures connected to these. Of course, different curves are obtained according to whether the well intersects less or more fractures than the expected number Ni . Also, fluid exchanges between the well and different parts of the fracture network may take place in parallel, but over different time scales. Furthermore, the drainage areas corresponding to the various parts of the network connected with the well can be quite different, which yields different onset times and durations for the pseudo steady fracture flow. Typical results for a network with density ρ = 6 are given in Fig. 4, for rw = 0.01R and two values of the conductivity σ = 102 and 103 . In this case, the actual number of intersections is equal to the expected value Ni = 6. The four flow regimes identified in Fig. 3 can be observed, i.e., transient and pseudo steady fracture flows, followed by transient and pseudo steady global flows. The decrease of the derivative w,DA at tDA ∼ 10−4 is commonly observed, and corresponds to the spread of the pressure wave in the ramified network, with decreasing hydraulic resistance. 5.3. Effective properties of the fractured porous medium Systematic simulations of drawdown tests have been run over a wide range of network density ρ and fracture permeability σ [4]. In the late pseudo steady stage, the well pressure always increases linearly with time. The results
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2277
Pw,DA Π
Pw,DA Π
w,DA
w,DA
P
→
Pw,DA→
w,DA
−1
10 −1
10
←Π
← Πw,DA
w,DA
−2
10 −2
10
−3
nI = 6, σ = 102
10
3
nI = 6, σ = 10
−3
−4
10 −6 10
−4
−2
10
0
10
10
10 −6 10
tʹ
−4
−2
10
10
a
0
10
tʹ
b
Fig. 4. Pressure Pw,DA (—) and its derivative w,DA = dPw,DA /d ln tDA (−·− · −) for a pressure drawdown test, for a fracture network with ρ = 6 in a unit cell with size L = 6R. The well intersects 6 fractures. Data are for rw = 0.01R, b = 10−3 ; σ = 102 (a) and σ = 103 (b). The analytical solution by [17] for a spacing equal to the mean value Si = R is given for the theoretical (−−− −) and corrected (··· · ·) Df (see [15]).
can be analyzed by comparison with the pressure response of a well in a homogeneous medium with permeability Kr and storage coefficient g Cg (see (10)) A Jw 2.2458 μJw ln 2 + ln (13) Pw (t) = t+ g C g A 4πKr rw CA It has been shown in [4] that as a rule, Kr is significantly larger than the fractured medium effective permeability Keff and strongly depends on the number of well/fractures intersections. Hence, Kr cannot be regarded as an intrinsic property of the fractured medium, since it depends on boundary conditions. On physical grounds, it can be expected that the well pressure response should depend on the far-field effective permeability Keff of the medium, and on the particular interactions of the well with the fracture network in a given situation. This can be modeled in a first approximation by introducing a skin factor S [9] in Eq. (13), A Jw 2.2458 μJw Pw,DA (t) = ln 2 + ln + 2S (14) t+ g Cg A 4πKeff rw CA S is positive if permeability near the well is smaller than in the surrounding medium.
15
S ρ’=6, σ’=10 ρ’=6, σ’=100 ρ’=6, σ’=1000 ρ’=4, σ’=10 ρ’=4, σ’=100 ρ’=4, σ’=1000 ρ’=1.5, σ’=10 ρ’=1.5, σ’=100 ρ’=1.5, σ’=1000
10
5
0
−5 0
1
2
3
4
5
6
sI / SI = NI / nI
Fig. 5. Skin factor S as a function of the spacing si /Si . Data are for a well with a radius rw = 0.01R in a randomly fractured medium, with various densities ρ and fracture conductivities σ . The broken lines correspond to the fit (15).
2278
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
Fig. 6. Three-dimensional view of the simulation domain for the SEH experimental site. The central part is a 500 × 500 × 100 m3 fractured zone (yellow); the margin zone is green, the black lines are wells; red, green and blue lines are fracture traces. The mesh contains more than 3 × 106 volume elements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
The skin factor is plotted in Fig. 5 as a function of the spacing between the well/fracture intersections normalized by its expected value, si /Si . A linear fit is possible, which provides a guide line to predict the performances of a well, given its connectivity with the fracture network, for highly conducting and well connected fractures si S = 2.5 − 1 − 1.6 (r = 0.963) (15) Si 6. Field simulations We present in this section two examples of applications on the field scale, which are numerically more demanding owing to the problem size. In the first example, we study the behaviour of the fractured aquifer at the Experimental Hydrological Site (SEH) of Poitiers. The computational domain covers 8 × 8 × 0.1 km3 (Fig. 6). A central 500 × 500 × 100 m3 region contains a discrete fracture network, which is meshed as described in Section 3. It is surrounded by a wide area where the fractured medium is replaced by an equivalent homogeneous medium, with the same effective properties as the central region. More than 30 vertical wells are located in the central region. The fracture network is generated stochastically. It contains three families of sub-vertical fractures, with their densities and orientations set according to the regional trends, and one family of sub-horizontal fractures, devised so that the statistics of the well fracture intersections comply with the actual field observations. Many simulations have been run, including single-well drawdown pumping tests and two-well injection/production tests. The other wells can be regarded as non-conducting, in which case they do not play any role, or as conducting
Fig. 7. Pressure distribution during a steady pumping at one of central wells, in the domain of Fig. 6.
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2279
Fig. 8. Drawdown curves for 16 observation wells during a steady pumping at one of the central wells. Colored curves are typical responses of observation wells recorded on the experimental site. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
conduits, which is more realistic. They correspond then to the observation wells which provide on the site an insight in the piezometric variations but also actively contribute and disturb the flow by creating additional man-made flow-paths. Fig. 7 presents an example of pressure field for a steady state pumping test, where the influence of the fracturation is clearly visible, in particular in the anisotropy that it introduces. Different kinds of response can be obtained depending on the well selected for the test. As shown in Section 5.1, it depends first on the degree of connection of the well with the fracture network. It also depends on the connectivity of the intercepted fractures with the rest of the network. Similarly, the response of an observation well also depends on its connectivity, and on the characteristics of the aquifer properties in its surrounding. Fig. 8 displays an example of the drawdown curves at all the observation wells when a pumping with a steady flow rate is performed in a central well. The numerical model reproduces well the variety of behaviors actually recorded on the site, where three kinds of typical responses can be obtained. The vertical fluxes in the pumping well and in some observation wells are detailed in Fig. 9, at successive instants during a steady drawdown test. The discontinuities in the profiles of Fig. 9a correspond to exchanges between the well and intercepting fractures. Direct transfers between the well and the matrix rock is negligible. At the pumping well the steady profile is reached almost instantaneously.
Fig. 9. Evolution with time of the vertical distribution of the vertical flow rate in the production well (a) and two observation wells at 70 m (b) and 160 m (c) from the pumping well. Variations of the colors from blue through cyan to red correspond to the increase of time. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
2280
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
Fig. 10. Pressure fields in the fractured oil reservoir after 2 (a) and 4 (b) years of production.
This is not true in distant observation wells. Note first that a significant flow rate can develop along the wells, which demonstrates their disturbing influence. For instance, the well in Fig. 9c transfers fluid from the top half of the aquifer to a fracture, and from its middle to its lower parts. In addition, the vertical flow profiles evolve with time, both in amplitude and in shape, before they settle at their final steady state, as seen in Fig. 9b. The drainage of the aquifer by the fracture which crosses the well near its top extremity is stronger at early time than in the final state, as a result of an evolution of the overall flow pattern in the whole domain. This illustrates the variety of time scales associated with the flow paths in the matrix rock and in the fracture network, and the four-dimensional character of the flow field. In the second example, a whole oil reservoir was meshed up to its boundaries where mixed conditions were imposed: no flow at the ceiling and part of the lateral boundaries, pressure conditions at the interface with an underlying aquifer and interaction with another neibouring reservoir at other lateral boundaries.
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2281
Several vertical, horizontal and oblique wells are included, and the fracturation is introduced according to actual field data for the main events on the kilometric scale, or from stochastic generation for more than 700 smaller scale ones, with varying hydraulic properties. The simulations allow to follow the evolution of the oil pressure field, when the pumping history is prescribed at each of the wells, either from actual chronicles or according to other scenarios for optimisation purposes. Fig. 10 shows an example of such fields after two and four years of production. 7. Concluding remarks The numerical model presented in this paper for 3D single phase compressible flow through fractured porous media was validated in simple reference cases and performs well in large scale complex situations. The flexibility of the numerical description of the discrete fracture network makes possible to directly simulate the operation of wells in a real setting, given field data relative to the fracture locations and characteristics. In addition, the generation of the numerical mesh prior to the introduction of the well allows for an easy optimisation of the well design, according to any relevant criterion. Although it was not illustrated here, the numerical model is ready to handle heterogeneous matrix or fracture properties, far from or in the absence of a well. This is an important feature, since spatial variations of the fracture aperture may have a significant influence on the exchanges with the matrix. Additional work is needed to accurately describe the fluid exchanges between a well and its immediate heterogeneous surroundings. Work is underway to improve the description of matrix/fracture exchanges during the fast transients that take place at the early stages of pumping tests. Another desirable extension of the model is the incorporation of two-phase flow, which will be addressed in the near future. References [1] P.M. Adler, J.-F. Thovert, Fractures and Fracture Networks, Kluwer Academic Publishers, Dordrecht, 1999. [2] G.I. Barenblatt, Y.P. Zheltov, Fundamental equations of filtration of homogeneous liquids in fissured rocks, Soviet Dokl. Akad. Nauk 13 (1960) 545–548. [3] I.I. Bogdanov, V.V. Mourzenko, J.-F. Thovert, P.M. Adler, Effective permeability of fractured porous media in steady state flow, Water Resour. Res. 39 (2003) 1023. [4] I.I. Bogdanov, V.V. Mourzenko, J.-F. Thovert, P.M. Adler, Pressure drawdown well tests in fracture porous media, Water Resour. Res. 39 (2003) 1021. [5] N.S. Boulton, T.D. Streltsova, Unsteady flow to a pumped well in a two-layered water-bearing formation, J. Hydrol. 35 (1977) 245–256. [6] N.S. Boulton, T.D. Streltsova, Unsteady flow to a pumped well in a fissured water-bearing formation, J. Hydrol. 35 (1977) 257–270. [7] I. Canamon, F.J. Elorza, R. Ababou, Thermo-hydro-mechanical simulation of a 3d fractured porous rock: preliminary study of coupled matrix-fracture hydraulics, Eurotherm Seminar 81, ET81–11, Albi, 2007. [8] Z.-X. Chen, Transport flow of slightly compressible fluids through double-porosity, double-permeability systems – a state-of-the-art review, Transport Porous Media 4 (1989) 97–116. [9] R.C. Earlougher, Advances in Well Test Analysis, SPE Monograph Series, New York, 1977. [10] J.P. Gwo, R. O’Brien, P.M. Jardine, Mass transfer in structured porous media: embedding mesoscale structure and microscale hydrodynamics in a two-region model, J. Hydrol. 208 (1998) 204–222. [11] O. Huseby, J.-F. Thovert, P.M. Adler, Geometry and topology of fracture systems, J. Phys. A30 (1997) 1415–1444. [12] N. Koudina, R. Gonzalez-Garcia, J.-F. Thovert, P.M. Adler, Permeability of three-dimensional fracture networks, Phys. Rev. E57 (1998) 4466–4479. [13] M. Onur, A. Satman, Interpretation of single-well pressure transient data from naturally fractured reservoirs, In Situ 22 (1998) 181–237. [14] J.C. Reis, Effect of fracture spacing distribution on pressure transient response in naturally fractured reservoirs, J. Petrol. Sci. Eng. 20 (1998) 31–47. [15] S.G. Shikaze, E.A. Sudicky, F.W. Schwartz, Density dependent solute transport in discretely-fractured geological media: is prediction possible? J. Contam. Hydrol. 34 (1998) 273–291. [16] D.T. Snow, Anisotropic permeability of fractured media, Water Resour. Res. 5 (1969) 1273–1289. [17] T.D. Streltsova, Well Testing in Heterogeneous Formations, An Exxon Monograph, John Wiley & Sons, New York, 1988. [18] J.R. Warren, P.J. Root, The behaviour of naturally fractured reservoirs, Soc. Pet. Eng. J. 228 (1963) 245–255. [19] G.J. Weir, Single-phase flow regimes in a discrete fracture model, Water Resour. Res. 35 (1999) 65–73. [20] C.R. Wilson, P.A. Witherspoon, Steady state flow in rigid networks of fractures, Water Resour. Res. 10 (1974) 328–335. [21] R. Young, Pressure transient in a double-porosity medium, Water Resour. Res. 28 (1992) 1261–1270.
Mathematics and Computers in Simulation 81 (2011) 2270–2281
Original article
Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations V.V. Mourzenko a,∗ , I.I. Bogdanov a,b,1 , J.-F. Thovert a , P.M. Adler b a
Institut PPRIME, SP2MI, BP 30179, 86962 Futuroscope Chasseneuil cedex, France b UPMC - Sisyphe, 4 place Jussieu, 75252 Paris cedex 05, France
Received 15 November 2009; received in revised form 8 June 2010; accepted 8 December 2010 Available online 21 December 2010
Abstract A general three-dimensional numerical model for single phase, slightly compressible flow through fractured porous media is described. It is based on a discrete fracture representation. Three sets of applications are presented. In the first one, pressure drawdown well tests in closed oil reservoirs are simulated for complex model situations where the well intercepts a random fracture network with various fracture densities and conductivities. Then, the hydrodynamic response of a fractured aquifer is investigated by simulating on the field-scale single- or two-well pumping tests in the Poitiers Hydrological Experimental Site. Finally, a complete field-scale simulation of the production history in an oil reservoir with multiple wells is presented. © 2010 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Porous medium; Fracture; Reservoir; Well flow
1. Introduction This work addresses single phase, slightly compressible flow through fractured porous media. A very general, three-dimensional numerical model based on a discrete fracture representation is proposed, together with applications involving a single or multiple interacting wells. Due to the specific transport properties of fractures, the flow through a naturally fractured porous medium differs drastically from that in a conventional porous medium. The key feature is that the porous matrix provides the main storage for the fluids while transport takes place mainly through the fracture system. Furthermore, matrix/fracture flow interactions govern many of the medium transport properties. Because of the complexity and diversity of most natural fracture systems, the determination of fractured porous media transport properties remains an open issue of great practical importance. For instance, the present numerical tool can be applied for the interpretation of well test data, in order to quantify the characteristics of a reservoir, or conversely, to optimise the design of a producing well, given the reservoir characteristics.
∗
Corresponding author. Tel.: +33 5 49 49 67 72. E-mail addresses: [email protected] (V.V. Mourzenko), [email protected] (I.I. Bogdanov), [email protected] (J.-F. Thovert), [email protected] (P.M. Adler). 1 Present address: CHLOE - UPPA Avenue de l’Université, Bat. B1, 3ème étage droite, BP 1155, 64013 Pau cedex, France. 0378-4754/$36.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2010.12.014
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2271
Most existing models of transient compressible flow through fractured porous media are based on the concept of double-porosity media [2,18,8,14,10]. A discrete approach allows a more detailed description of the matrix/fracture flow interaction accounting for the flow behavior in individual elements of the medium, but it has been applied so far only to regular system of fractures [16,20,5,21,15]. We presented in [4] an extension to unsteady compressible flow of the numerical models of [12,3], which make use of a full three-dimensional description of random discrete fracture systems. Its main features are recalled here, as well as results of single-well tests simulations in medium-sized systems (hundreds of fractures), focused on the well response. Then, new applications on a much larger scale (thousands of fractures) are presented, involving many wells which can be simultaneously producing or simply observing. Attention is extended to the transient complex three-dimensional flow field in the fractured aquifer or reservoir. For instance, the possible disturbing effects of supposedly passive wells are demonstrated. 2. Governing equations The fractured porous medium can be represented as an arbitrary, generally random, set of fractures embedded into a solid porous matrix. At the large scale, each fracture can be viewed as a vanishingly thin layer, with singular transport properties. For instance, a steep pressure gradient may exist across a poorly conducting fracture, which results in an apparent pressure jump at the macroscopic scale. In addition, the typical time for pressure variations in transient flows is quite different in the matrix and in the fractures. The discrete fracture model used in the present work explicitly accounts for these features. Therefore, it allows detailed investigations of transitory flow regimes, and of the influence of the well intersections with the fracture network. 2.1. Flow equations Consider a porous matrix with a bulk permeability Km [L2 ] which may vary with space. Darcy law for the local seepage velocity v and the mass conservation for slightly compressible flow can be written in the matrix as v=−
Km ∇P μ
(a),
m Cm
∂P + ∇ · v = δw Jw ∂t
(b)
(1)
where μ is the viscosity, P is the pressure, m and Cm are the matrix porosity and total compressibility. Jw [L2 T−1 ] represents the exchanges with the well, at a location given by the Dirac function δw [L−2 ]. On this scale of description, the well appears as a line without thickness. If m Cm is constant, a matrix pressure diffusivity Dm can be defined and (1) can be written as a difusion equation for the pressure ∂P δ w Jw − ∇ · (Dm ∇P) = , ∂t m C m
Dm =
Km μm Cm
(2)
We assume that the hydraulic properties of a fracture can be characterized by two effective coefficients, namely a conductivity σ[L3 ] and a cross resistance ω[L−1 ] (see [3]), which relate the in-plane the flow rate j to the surface pressure gradient ∇ s P and the seepage velocity v⊥ of the net flow crossing the fracture to the pressure drop [P] across it by σ j = − ∇s P μ
(a),
v⊥ = −
1 [P] μω
(b).
(3)
For illustration, the fracture can be viewed as a plane channel of aperture b, filled with a porous material with permeability Kf , porosity f and total compressibility Cf . Then, σ and ω are given by σ = bKf and ω = b/Kf . A non zero resistance ω may exist even for totally open channels, if their walls are partially clogged by a deposited chemical [13]. Note that equations similar to (3) for the fracture hydraulic properties have been obtained by [7] by an upscaling procedure. The continuity equation for the flow through a fracture reads as bf Cf
∂P + ∇s · j = (v− − v+ ) · n + bδw Jw ∂t
(4)
2272
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
where v± is the seepage velocity in the matrix on either sides of the fracture, given by (1). Again, if bf Cf is constant in a fracture, a pressure diffusivity Df can be defined, and Eq. (4) can be written as ∂P (v− − v+ ) · n δ w Jw + , − ∇s · (Df ∇s P) = bf Cf f Cf ∂t
Df =
σ bμf Cf
(5)
The assumptions of uniform m Cm and bf Cf used to write down the diffusion equations (2) and (5) are by no means a requirement for the numerical model, since the discretized equations are really based on the general equations (1) and (4). The problem is governed by four dimensionless numbers relative to the fractures (b , σ , ω , η), in addition to the well properties, b =
b , R
σ =
σ , RKm
ω =
ωKm , R
η=
f Cf b m Cm R
(6)
In the present work, we only consider very conducting fractures and neglect pressure variations along the well. Unless otherwise stated, all the results are given for fractures without cross-resistance to flow (ω → 0), with large aspect ratio (b = η = 10−3 ), and a large fracture/matrix conductivity contrast (σ 1). The unit scale R is a characteristic lateral extent of the fractures. Other choices are possible. In particular, the well radius rw can be used as the characteristic length scale instead of the fracture size R. This is appropriate to describe the very early stage of a well test, when the flow is confined very near to the well [9]. Alternatively, the size of the drainage area A of a closed reservoir can be used as length unit (see Eq. (11)). 2.2. Overall boundary conditions The model can accommodate the usual boundary conditions in reservoir engineering, namely of Dirichlet type, with an imposed far-field pressure P∞ , or of Neumann type, with no-flux conditions at the boundaries of a closed reservoir. For systematic studies of the interactions of a single well with a fractured formation, a simple geometrical setting is used, with a parallelepipedic cell. The well is set along a vertical line throughout the cell, and spatial periodicity is applied in the x-, y-directions. Hence, this situation corresponds to a periodic square array of infinite vertical wells, or to a single well in a closed reservoir. All the simulations reported here correspond to a pressure drawdown test, i.e., to a well producing at constant rate Qw from a field initially at rest. The initial pressure is arbitrary, and can be taken as P0 = 0 without loss of generality. Of course, other settings can be devised to address any specific field case as illustrated in Section 6. 3. Numerical model We consider here both simple deterministic test configurations and more realistic and complex randomly fractured porous media. Whereas specific discretization techniques could be used in the first situation, for instance with local grid refinements, we are really interested in the latter case where geometrical randomness requires great flexibility. The numerical model is a generalization of [12,3] for fracture networks and fractured porous media. More details can be found in [4]. 3.1. Triangulation of fractured media and description of the well geometry The fracture network and the space located in between the fractures are successively triangulated by an advancing front technique. The scale of the discretization is controlled by a typical spacing δm between the mesh points. Threedimensional views of examples of triangulated fractured medium are displayed in Fig. 1. In a subsequent step, the well is superimposed on the three-dimensional unstructured mesh. Note that the initial triangulation is performed without considering the well. An alternative is to take the well into account from the start during the construction of the mesh, so that nodes would be distributed right on the well. However, the 3D triangulation of complex media is a time-consuming step. The present procedure allows to perform it only once for a given medium, and then to vary the well geometry or location at negligible cost.
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2273
Fig. 1. Three-dimensional views of triangulated periodic fractured media. The example in (a) contains Nfr = 32 hexagonal fractures with circumscribed radius R. The cell size is L = 4R, and δm = R/4. The tetrahedral volume elements with their centers in the cubic unit cell − 2R ≤ x, y, z ≤+2R are displayed. The protuding fractures sit astride the boundaries with the neighboring cells. The intersections of the fractures (thick lines) and of the tetrahedral volume elements with the cell boundaries are shown in the sample in (b), with L = 6R, Nfr = 166 and δm = R/4. This mesh contains 17,671 nodes and 111,556 tetrahedra.
The typical size δm of the volume and surface elements has to be small compared to the typical size R of the fractures, for a good rendering of the fractured medium geometry, and large compared to the well radius rw , since the well is considered as a line without thickness in the model.
3.2. Finite volume formulation The rock matrix is represented by tetrahedral volume elements, with given properties m , Cm and Km , and the fractures by triangular surface elements with given properties f , Cf , b, σ and ω. All these parameters can be set independently, on a per-element basis. The pressure is evaluated at the vertices of the surface and volume elements (◦, in Fig. 2c). The well is represented by a series of segments in the volume elements that it crosses. The pressure in the well is evaluated at the mid-points of these segments, and at the intersection of the well with fractures (• and , in Fig. 2c). The equations for the node pressures are obtained by integration of the balance Eqs. (1) and (4) over a control volume around each mesh point. Several types of control volumes exist according to the location of the mesh point, in the bulk of the rock matrix or on a fracture, with or without a well nearby. The simplest case of a node point in the matrix, without well nearby, is depicted in Fig. 2a. The continuity equation (1b) is integrated over a volume , made up of one fourth of all the incident tetrahedral volume elements
∂P ∂P m Cm m Cm dv + n · vds = 0 + ∇ · v dv = ∂t ∂t ∂
(7)
Since the matrix properties m , Cm and Km are taken as uniform over each tetrahedral volume element, the seepage velocity v in each tetrahedron is piecewise constant over ∂, and can be evaluated from (1), where the pressure gradient is deduced from the pressures at the four vertices.
2274
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
Pressure may take different values P± in the matrix on either sides of a fracture. A third value of the pressure Pf is introduced, in the middle plane of the fracture, which is used for Eq. (4). Balance equations are written over the three control volumes depicted in Fig. 2b, ± on either side of the fracture and f in the fracture, ±
m C m
∂P ± dv + ∂t
∂± −F
n · vds = q±
f
bf Cf
∂Pf ds + ∂t
∂f
n · jdl = q− − q+
(8)
The flux j in each triangle intercepting f is evaluated from (3a) where ∇ s P is a linear combination of the pressures at the triangle vertices, and the exchange terms q± are modeled according to (3b). The balance equations (8) yield three linear equations with the same form as for (7). More complex situations occur at the intersection of two or three fractures, but they are handled in a similar way. The well is treated as a one-dimensional highly conductive system superimposed on the medium. Accordingly, it can be viewed as a line source or sink in the fractured porous medium, decomposed into elementary parts corresponding to its intersections with the grid surface and volume elements. The source terms involving Jw in (1) and (4) have to be added to the balance equations for the control volumes (7) and (8). The flux Jw can be related to the well pressure Pw by decomposing the pressure field near the well into the superposition of a global trend, induced for instance by a large scale transverse flow and represented by a constant pressure gradient G in the mesh element, and of a disturbance induced by the presence of the well, modeled according to the analytical solution for a straight well in an unbounded medium. A finite volume formulation can also be implemented to describe the flow along the well, but in the present case, where pressure variations along the well are neglected, it reduces into an integral condition on the pumping rate Qw . 3.3. Solution of the discretized equations Time derivatives are discretized to first order, in a fully implicit formulation. Since the matrix for the set of linear equations to be solved at each time step is not symmetric, due to the exchange terms with the well, a biconjugate gradient algorithm is used. There is no stability restriction for the choice of the time step δt. However, it has to be set small enough in order to sample the transients of the flow that we are interested in, such as the early stages of the flow around a well, or the propagation of a wave pressure along a fracture which can be much faster than in the imbedding matrix. This implies criteria of the form Dδt < h2 where h and D are the relevant length scale and diffusivity. For instance, D = Df and h is the fracture extension when tracking the travel of a pressure wave along a fracture. The wide range of length scales and the contrast between the diffusivities of the fractures and matrix make a careful management of the time step necessary. In practice, different time steps have to be used at various stages of the simulations, depending on the relative contributions of the fractures and matrix to the flow and on its current transitory or quiescent character.
Fig. 2. Control volumes around a point in the matrix (a) and on a fracture (b). Two-dimensional illustration of the introduction of the well in the mesh (c).
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2
10
,Π
P
w,DA
w,DA [C] [D] Pseudo− "Hydraulic steady fracture" flow fracture flow
1
10
2275
[A] Early time fracture flow
1
[B] Transient fracture flow
P
w,DA
→
0
0.46
10
[E] Pseudo−steady global flow
← Πw,DA −1
10
1/25 −−−−−> 1/26 −−−−−−−−−−−−−>
0.78
tʹ
−2
10
−5
10
−3
10
−1
10
1
10
Fig. 3. The well pressure Pw,DA and its logarithmic time derivative w,DA versus the dimensionless time tDA , for a straight well intersecting a series of parallel plane fractures, with spacing H = L, in a closed reservoir with a square drainage area L × L. Data are for rw = 5 10−3 R, σ = 50, b = η = 5 10−4 , and δm = R/8. Curves correspond to Pw,DA and w,DA with (—) or without (−·− · −) well/matrix exchanges, and to the analytical prediction (10) for Pw,DA based on the raw (−−− −) or corrected (··· · ·) value of Df (see [3]). Straight lines indicate fitted slopes 0.78, 0.46 and 1. Horizontal arrows correspond to the values of the plateau of w,DA predicted by (9), with or without well/matrix exchanges.
4. Comparison with analytical solutions A few simple cases have analytical solutions, which have been used in [4] as references to assess the performance of the numerical model, such as the so-called drawdown test for a single vertical well in a uniform porous medium without any fracture, initially at rest, at a given constant pumping rate Jw per unit length of the well. The exact solution reduces after a short initial period to the well-known approximate logarithmic form [9] 4Dm t μJw 100r 2 ln 2 + 0.81 , (9) P(r, t) = ±1% for t > 4πKm r 4Dm Another basic solution is the pseudosteady long-time solution for a finite region with drainage area A [6,17] A Jw t μJw 2.2458 ln 2 + ln P(r, t) = , (10) + m Cm A 4πKm r CA where the shape factor CA depends on the reservoir geometry and on the location of the well. The term independent of time in (10) allows to evaluate the medium properties in the finite flow region. The results (9) and (10) can be straightforwardly transposed to the radial flow in a fracture in an impervious matrix. According to (10), the results can be expressed in terms of the dimensionless time tDA based on the drainage area A = L2 and of the dimensionless well pressure Pw,DA based on the imposed flow rate Jw tDA =
4Dm t , L2
Pw,DA =
4πKm Pw μJw
(11)
The results of the numerical simulations agree with the theoretical predictions (9) and (10) within a few percents for time and spatial steps δtDA ∼ 0.01 and δm /L ∼ 0.1. A well intersecting a single fracture has also given rise to analytical results [6,19] and can constitute the basis for modeling more complex situations. Our simulations, with or without taking into account the well/matrix exchanges, account for the various well-known successive regimes (Fig. 3). At early times, flow is mainly due to fracture compressibility. After a very brief transient (region [A]), a plateau in the curves for w,DA [B] corresponds to the prediction of (9); it is followed by a steeper increase [C], where Pw,DA grows as t0.78 . If the porous matrix were impervious, the pressure would increase linearly
2276
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
with time at this stage, according to (10). Later on, flow in the matrix becomes significant, and another regime is reached [D]. The slope of the curve of w,DA is about 0.46, up to tDA = 1/4, which is of the order of the typical time tDA ∼ 1 (t ∼ H2 /4Dm ) for the pressure wave to travel from the fracture to the middle of the matrix layer. This regime corresponds to a power dependence of Pw on time as t0.46 , as observed with hydraulic fractures, rather than to the halfed logarithmic law of [6], again because of the impervious domain boundaries. Finally, another pseudo steady regime [E] is reached for tDA > 1/4, with a law of the form of (10). Two hydraulic conductivities can be measured from the well pressure curve; the plateau [B] of w,DA provides the fracture conductivity (see (9)), and an extrapolation at t = 0 of the linear increase of Pw in the pseudo steady regime [E] provides an estimate of a global effective permeability Keff for the fractured medium (see (10)). 5. Well through a randomly fractured porous medium 5.1. Fracture network model We consider here randomly fractured media, as modeled by [11]. The fractures are plane polygonal objects, randomly located and oriented in space (see Fig. 1). The density of fractures is quantified by the dimensionless parameter ρ , which is the number of fractures per excluded volume Vex [11,1]. It also quantifies the connectivity of the fracture network, since ρ is the mean number of intersections per fracture. In particular, ρ determines the percolation of the fracture network, with a percolation threshold ρc ≈ 2.3, which has a dramatic influence on the macroscopic permeability of the fractured medium [3]. For very conducting fractures, it is expected that the number ni of intersections of the well with the network is a key parameter for the flow, since it conditions the fluid exchanges. Since the fractures are located according to a Poisson process, their intersections with the well are also Poissonian. Thus, their spacings si obey a negative exponential law f (si ) =
1 −si /Si e , Si
Si =
2 S
(12)
where S is the volumetric surface area of fractures and Si is the mean spacing. For a well of length L, the mean expected number of intersections is Ni = L/Si . 5.2. Examples of pressure drawdown curves As a rule, all the typical flow regimes present in Fig. 3 can be observed in the drawdown curves for wells in randomly fractured media. In particular, the two latest regimes, namely, the “hydraulic fracture” and the pseudo-steady flow, can always be clearly identified. However, the early and intermediate regimes are often blurred or do not show up at all, and additional features can appear in the curves. First, during the period when flow takes place mainly in the fractures, the pressure wave starts propagating in the fractures which cross the well, then spreads out in all the fractures connected to these. Of course, different curves are obtained according to whether the well intersects less or more fractures than the expected number Ni . Also, fluid exchanges between the well and different parts of the fracture network may take place in parallel, but over different time scales. Furthermore, the drainage areas corresponding to the various parts of the network connected with the well can be quite different, which yields different onset times and durations for the pseudo steady fracture flow. Typical results for a network with density ρ = 6 are given in Fig. 4, for rw = 0.01R and two values of the conductivity σ = 102 and 103 . In this case, the actual number of intersections is equal to the expected value Ni = 6. The four flow regimes identified in Fig. 3 can be observed, i.e., transient and pseudo steady fracture flows, followed by transient and pseudo steady global flows. The decrease of the derivative w,DA at tDA ∼ 10−4 is commonly observed, and corresponds to the spread of the pressure wave in the ramified network, with decreasing hydraulic resistance. 5.3. Effective properties of the fractured porous medium Systematic simulations of drawdown tests have been run over a wide range of network density ρ and fracture permeability σ [4]. In the late pseudo steady stage, the well pressure always increases linearly with time. The results
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2277
Pw,DA Π
Pw,DA Π
w,DA
w,DA
P
→
Pw,DA→
w,DA
−1
10 −1
10
←Π
← Πw,DA
w,DA
−2
10 −2
10
−3
nI = 6, σ = 102
10
3
nI = 6, σ = 10
−3
−4
10 −6 10
−4
−2
10
0
10
10
10 −6 10
tʹ
−4
−2
10
10
a
0
10
tʹ
b
Fig. 4. Pressure Pw,DA (—) and its derivative w,DA = dPw,DA /d ln tDA (−·− · −) for a pressure drawdown test, for a fracture network with ρ = 6 in a unit cell with size L = 6R. The well intersects 6 fractures. Data are for rw = 0.01R, b = 10−3 ; σ = 102 (a) and σ = 103 (b). The analytical solution by [17] for a spacing equal to the mean value Si = R is given for the theoretical (−−− −) and corrected (··· · ·) Df (see [15]).
can be analyzed by comparison with the pressure response of a well in a homogeneous medium with permeability Kr and storage coefficient g Cg (see (10)) A Jw 2.2458 μJw ln 2 + ln (13) Pw (t) = t+ g C g A 4πKr rw CA It has been shown in [4] that as a rule, Kr is significantly larger than the fractured medium effective permeability Keff and strongly depends on the number of well/fractures intersections. Hence, Kr cannot be regarded as an intrinsic property of the fractured medium, since it depends on boundary conditions. On physical grounds, it can be expected that the well pressure response should depend on the far-field effective permeability Keff of the medium, and on the particular interactions of the well with the fracture network in a given situation. This can be modeled in a first approximation by introducing a skin factor S [9] in Eq. (13), A Jw 2.2458 μJw Pw,DA (t) = ln 2 + ln + 2S (14) t+ g Cg A 4πKeff rw CA S is positive if permeability near the well is smaller than in the surrounding medium.
15
S ρ’=6, σ’=10 ρ’=6, σ’=100 ρ’=6, σ’=1000 ρ’=4, σ’=10 ρ’=4, σ’=100 ρ’=4, σ’=1000 ρ’=1.5, σ’=10 ρ’=1.5, σ’=100 ρ’=1.5, σ’=1000
10
5
0
−5 0
1
2
3
4
5
6
sI / SI = NI / nI
Fig. 5. Skin factor S as a function of the spacing si /Si . Data are for a well with a radius rw = 0.01R in a randomly fractured medium, with various densities ρ and fracture conductivities σ . The broken lines correspond to the fit (15).
2278
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
Fig. 6. Three-dimensional view of the simulation domain for the SEH experimental site. The central part is a 500 × 500 × 100 m3 fractured zone (yellow); the margin zone is green, the black lines are wells; red, green and blue lines are fracture traces. The mesh contains more than 3 × 106 volume elements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
The skin factor is plotted in Fig. 5 as a function of the spacing between the well/fracture intersections normalized by its expected value, si /Si . A linear fit is possible, which provides a guide line to predict the performances of a well, given its connectivity with the fracture network, for highly conducting and well connected fractures si S = 2.5 − 1 − 1.6 (r = 0.963) (15) Si 6. Field simulations We present in this section two examples of applications on the field scale, which are numerically more demanding owing to the problem size. In the first example, we study the behaviour of the fractured aquifer at the Experimental Hydrological Site (SEH) of Poitiers. The computational domain covers 8 × 8 × 0.1 km3 (Fig. 6). A central 500 × 500 × 100 m3 region contains a discrete fracture network, which is meshed as described in Section 3. It is surrounded by a wide area where the fractured medium is replaced by an equivalent homogeneous medium, with the same effective properties as the central region. More than 30 vertical wells are located in the central region. The fracture network is generated stochastically. It contains three families of sub-vertical fractures, with their densities and orientations set according to the regional trends, and one family of sub-horizontal fractures, devised so that the statistics of the well fracture intersections comply with the actual field observations. Many simulations have been run, including single-well drawdown pumping tests and two-well injection/production tests. The other wells can be regarded as non-conducting, in which case they do not play any role, or as conducting
Fig. 7. Pressure distribution during a steady pumping at one of central wells, in the domain of Fig. 6.
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2279
Fig. 8. Drawdown curves for 16 observation wells during a steady pumping at one of the central wells. Colored curves are typical responses of observation wells recorded on the experimental site. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
conduits, which is more realistic. They correspond then to the observation wells which provide on the site an insight in the piezometric variations but also actively contribute and disturb the flow by creating additional man-made flow-paths. Fig. 7 presents an example of pressure field for a steady state pumping test, where the influence of the fracturation is clearly visible, in particular in the anisotropy that it introduces. Different kinds of response can be obtained depending on the well selected for the test. As shown in Section 5.1, it depends first on the degree of connection of the well with the fracture network. It also depends on the connectivity of the intercepted fractures with the rest of the network. Similarly, the response of an observation well also depends on its connectivity, and on the characteristics of the aquifer properties in its surrounding. Fig. 8 displays an example of the drawdown curves at all the observation wells when a pumping with a steady flow rate is performed in a central well. The numerical model reproduces well the variety of behaviors actually recorded on the site, where three kinds of typical responses can be obtained. The vertical fluxes in the pumping well and in some observation wells are detailed in Fig. 9, at successive instants during a steady drawdown test. The discontinuities in the profiles of Fig. 9a correspond to exchanges between the well and intercepting fractures. Direct transfers between the well and the matrix rock is negligible. At the pumping well the steady profile is reached almost instantaneously.
Fig. 9. Evolution with time of the vertical distribution of the vertical flow rate in the production well (a) and two observation wells at 70 m (b) and 160 m (c) from the pumping well. Variations of the colors from blue through cyan to red correspond to the increase of time. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
2280
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
Fig. 10. Pressure fields in the fractured oil reservoir after 2 (a) and 4 (b) years of production.
This is not true in distant observation wells. Note first that a significant flow rate can develop along the wells, which demonstrates their disturbing influence. For instance, the well in Fig. 9c transfers fluid from the top half of the aquifer to a fracture, and from its middle to its lower parts. In addition, the vertical flow profiles evolve with time, both in amplitude and in shape, before they settle at their final steady state, as seen in Fig. 9b. The drainage of the aquifer by the fracture which crosses the well near its top extremity is stronger at early time than in the final state, as a result of an evolution of the overall flow pattern in the whole domain. This illustrates the variety of time scales associated with the flow paths in the matrix rock and in the fracture network, and the four-dimensional character of the flow field. In the second example, a whole oil reservoir was meshed up to its boundaries where mixed conditions were imposed: no flow at the ceiling and part of the lateral boundaries, pressure conditions at the interface with an underlying aquifer and interaction with another neibouring reservoir at other lateral boundaries.
V.V. Mourzenko et al. / Mathematics and Computers in Simulation 81 (2011) 2270–2281
2281
Several vertical, horizontal and oblique wells are included, and the fracturation is introduced according to actual field data for the main events on the kilometric scale, or from stochastic generation for more than 700 smaller scale ones, with varying hydraulic properties. The simulations allow to follow the evolution of the oil pressure field, when the pumping history is prescribed at each of the wells, either from actual chronicles or according to other scenarios for optimisation purposes. Fig. 10 shows an example of such fields after two and four years of production. 7. Concluding remarks The numerical model presented in this paper for 3D single phase compressible flow through fractured porous media was validated in simple reference cases and performs well in large scale complex situations. The flexibility of the numerical description of the discrete fracture network makes possible to directly simulate the operation of wells in a real setting, given field data relative to the fracture locations and characteristics. In addition, the generation of the numerical mesh prior to the introduction of the well allows for an easy optimisation of the well design, according to any relevant criterion. Although it was not illustrated here, the numerical model is ready to handle heterogeneous matrix or fracture properties, far from or in the absence of a well. This is an important feature, since spatial variations of the fracture aperture may have a significant influence on the exchanges with the matrix. Additional work is needed to accurately describe the fluid exchanges between a well and its immediate heterogeneous surroundings. Work is underway to improve the description of matrix/fracture exchanges during the fast transients that take place at the early stages of pumping tests. Another desirable extension of the model is the incorporation of two-phase flow, which will be addressed in the near future. References [1] P.M. Adler, J.-F. Thovert, Fractures and Fracture Networks, Kluwer Academic Publishers, Dordrecht, 1999. [2] G.I. Barenblatt, Y.P. Zheltov, Fundamental equations of filtration of homogeneous liquids in fissured rocks, Soviet Dokl. Akad. Nauk 13 (1960) 545–548. [3] I.I. Bogdanov, V.V. Mourzenko, J.-F. Thovert, P.M. Adler, Effective permeability of fractured porous media in steady state flow, Water Resour. Res. 39 (2003) 1023. [4] I.I. Bogdanov, V.V. Mourzenko, J.-F. Thovert, P.M. Adler, Pressure drawdown well tests in fracture porous media, Water Resour. Res. 39 (2003) 1021. [5] N.S. Boulton, T.D. Streltsova, Unsteady flow to a pumped well in a two-layered water-bearing formation, J. Hydrol. 35 (1977) 245–256. [6] N.S. Boulton, T.D. Streltsova, Unsteady flow to a pumped well in a fissured water-bearing formation, J. Hydrol. 35 (1977) 257–270. [7] I. Canamon, F.J. Elorza, R. Ababou, Thermo-hydro-mechanical simulation of a 3d fractured porous rock: preliminary study of coupled matrix-fracture hydraulics, Eurotherm Seminar 81, ET81–11, Albi, 2007. [8] Z.-X. Chen, Transport flow of slightly compressible fluids through double-porosity, double-permeability systems – a state-of-the-art review, Transport Porous Media 4 (1989) 97–116. [9] R.C. Earlougher, Advances in Well Test Analysis, SPE Monograph Series, New York, 1977. [10] J.P. Gwo, R. O’Brien, P.M. Jardine, Mass transfer in structured porous media: embedding mesoscale structure and microscale hydrodynamics in a two-region model, J. Hydrol. 208 (1998) 204–222. [11] O. Huseby, J.-F. Thovert, P.M. Adler, Geometry and topology of fracture systems, J. Phys. A30 (1997) 1415–1444. [12] N. Koudina, R. Gonzalez-Garcia, J.-F. Thovert, P.M. Adler, Permeability of three-dimensional fracture networks, Phys. Rev. E57 (1998) 4466–4479. [13] M. Onur, A. Satman, Interpretation of single-well pressure transient data from naturally fractured reservoirs, In Situ 22 (1998) 181–237. [14] J.C. Reis, Effect of fracture spacing distribution on pressure transient response in naturally fractured reservoirs, J. Petrol. Sci. Eng. 20 (1998) 31–47. [15] S.G. Shikaze, E.A. Sudicky, F.W. Schwartz, Density dependent solute transport in discretely-fractured geological media: is prediction possible? J. Contam. Hydrol. 34 (1998) 273–291. [16] D.T. Snow, Anisotropic permeability of fractured media, Water Resour. Res. 5 (1969) 1273–1289. [17] T.D. Streltsova, Well Testing in Heterogeneous Formations, An Exxon Monograph, John Wiley & Sons, New York, 1988. [18] J.R. Warren, P.J. Root, The behaviour of naturally fractured reservoirs, Soc. Pet. Eng. J. 228 (1963) 245–255. [19] G.J. Weir, Single-phase flow regimes in a discrete fracture model, Water Resour. Res. 35 (1999) 65–73. [20] C.R. Wilson, P.A. Witherspoon, Steady state flow in rigid networks of fractures, Water Resour. Res. 10 (1974) 328–335. [21] R. Young, Pressure transient in a double-porosity medium, Water Resour. Res. 28 (1992) 1261–1270.