Numerical simulations on the influence of matrix diffusion to carbon sequestration in double porosity fissured aquifers

International Journal of Greenhouse Gas Control 3 (2009) 431–443

Contents lists available at ScienceDirect

International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc

Numerical simulations on the influence of matrix diffusion to carbon sequestration in double porosity fissured aquifers Ju´lio F. Carneiro * Geophysical Centre of E´vora, Geosciences Department, University of E´vora, Rua Roma˜o Ramalho 59, 7000 E´vora, Portugal

A R T I C L E I N F O

A B S T R A C T

Article history: Received 24 August 2008 Received in revised form 14 December 2008 Accepted 13 February 2009 Available online 27 March 2009

The double porosity model for fissured rocks, such as limestones and dolomites, has some features that may be relevant for carbon sequestration. Numerical simulations were conducted to study the influence of matrix diffusion on the trapping mechanisms relevant for the long-term fate of CO2 injected in fissured rocks. The simulations show that, due to molecular diffusion of CO2 into the rock matrix, dissolution trapping and hydrodynamic trapping are more effective in double porosity aquifers than in an equivalent porous media. Mineral trapping, although assessed indirectly, is also probably more relevant in double porosity aquifers due to the larger contact surface and longer contact time between dissolved CO2 and rock minerals. However, stratigraphic/structural trapping is less efficient in double porosity media, because at short times CO2 is stored only in the fissures, requiring large aquifer volumes and increasing the risk associated to the occurrence of imperfections in the cap-rock through which leakage can occur. This increased risk is also a reality when considering storage in aquifers with a regional flow gradient, since the CO2 free-phase will move faster due to the higher flow velocities in fissured media and discharge zones may be reached sooner. ß 2009 Elsevier Ltd. All rights reserved.

Keywords: Carbon sequestration Fissured aquifers Double porosity Matrix diffusion Numerical simulations

1. Introduction The capture of CO2 from point sources and its storage in geological formations has the potential for becoming a key strategy in decreasing the greenhouse gases emissions and mitigate the climate change effects (IPCC, 2005). This strategy, often referred as Carbon Capture and Storage (CCS), relies on the identification of geological formations that are able to guarantee a long-term safe and reliable storage of the CO2. Suitable geological formations are coal seams, depleted oil and gas reservoirs and deep saline aquifers (Celia et al., 2007; IPCC, 2005; Lokhorst and Wildenborg, 2005). Of these, saline aquifers provide undoubtedly the largest storage volume, estimated around 400–10,000 Gtons of CO2 (Castello and Bouchi-Lamontagne, 2007). The screening criteria for selecting reservoirs or saline aquifers include minimum values of permeability and porosity (Bachu, 2003; Bentham and Kirby, 2005) to ensure that injectivity and storage capacity are adequate. Those criteria are usually met by porous sedimentary rocks, such as sandstones, generally regarded as the most suitable formations for carbon storage.

* Tel.: +351 266745301; fax: +351 266745397. E-mail address: [email protected]. 1750-5836/$ – see front matter ß 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijggc.2009.02.006

Nevertheless, there are some commercial and pilot injection sites in fissured carbonate rocks, such as limestones and dolomites, where CO2 is being injected or is considered as part of Enhanced Oil Recovery (EOR) or Enhanced Gas Recovery (EGR) programs. Such is the case of the Weyburn-Midale field in Canada, where 15–20 millions tons of CO2 are expected to be injected as part of an EOR program (Lombardi et al., 2006), or of the Casablanca field in Spain, where REPSOL is studying the possibility of injecting CO2 (Le Gallo, 2007). Limestones, dolomites and other intensively fissured rocks can be described as a double porosity media (Barenblatt et al., 1960), where flow is mainly along fissures but with a potentially large storage volume in the rock matrix. This study is a preliminary assessment of the relevance of the double porosity behaviour to the long-term fate of CO2 injected in a fissured saline aquifer. The study is based on numerical simulations conducted with the TOUGH2 generic purpose simulator (Pruess et al., 1999). The paper starts with a general description of the main features of the double porosity media, followed by the details about the methodology and scenarios adopted. The following section presents the results of the numerical simulation and finally a discussion is made about the relevance of the double porosity behaviour to the long-term fate of the injected CO2. In the simulations presented, and throughout this paper, it is considered matrix diffusion as fickian transport of components in

432

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

the rock matrix due to molecular diffusion. Therefore, the simulations do not intend to represent the behaviour in fissured media where adsorption or surface diffusion (transport through physically adsorbed layers) dominates the CO2 trapping mechanisms, such as in coal seams (Shi and Durucan, 2005), and rather aim to represent the behaviour of highly fissured sedimentary rocks, such as limestones and dolomites. 2. Background on the double porosity model The concept of the double porosity media was first introduced by Barenblatt et al. (1960) and by Warren and Root (1963), and it is based on considering that fissured rock aquifers can be regarded as two overlapping continua: the fissure system where flow takes place; and the rock matrix, much less permeable than fissures, but providing most of the storage capacity (Fig. 1). Since both fissures and matrix are regarded as continua, several authors were able to find more or less complex analytical solutions for flow and mass transport (e.g., Barker, 1982; Barker et al., 2000; Moench, 1984; Sudicky and Frind, 1982). The main feature of the double porosity media with respect to phases, components or heat transport is the matrix diffusion process, in which a phase, component or heat can migrate by molecular diffusion from/to the fissures into/from the matrix (Fig. 2). The effects of the matrix diffusion are several: (a) matrix porosity becomes accessible to dissolved components and storage increases; (b) apparent retardation with respect to other components in the fissures, due to migration to the matrix and access to a larger volume of reaction sites; (c) tailing on breakthrough curves, where components take a long time to come out of the matrix because diffusion coefficients tend to be small (Barker, 1985b; Carrera et al., 1998; Maloszewski and Zuber, 1993). The double porosity model ultimately became of widespread use in oil reservoir engineering (e.g., Di Donato and Blunt, 2004; Kazemi and Gilman, 1993; Sarma and Aziz, 2006), where it was first enunciated, but it is also applied in groundwater flow studies (e.g., Barker, 1985a; Moench, 1984), contaminant transport analysis (e.g., Maloszewski and Zuber, 1985, 1993), underground waste repositories (e.g., Hadermann and Ro¨sel, 1985; Lever and Bradbury, 1985) and in geothermal energy studies (Kolditz, 1995). The model has been tested in those branches of science and technology and has proved, for the proper geologic conditions, remarkably useful and reliable. Nevertheless, the specific features of double porosity, most notably matrix diffusion, have not been addressed with respect to the study of carbon storage, although they may play a significant role in the long-term fate of CO2.

Fig. 1. Schematic illustration of the double porosity concept. (a) Naturally fissured rock, (b) Warren-Root’s idealized three-dimensional cubic fissure system; (c) idealized slab fissure system. From Kruseman and de Ridder (1994). Fissures provide any significant permeability, while most of the storage capacity is provided by the intact blocks of rock (the matrix).

Fig. 2. Schematic illustration of the matrix diffusion process. After Carrera et al. (1998). In the double porosity model transfer across the interface fissures/matrix is diffusive, imposed by concentration gradients (reflected in a molecular diffusion coefficient) and advective transport is considered negligible in the matrix.

2.1. Relevant physical processes The main physical processes active in the CO2 flow and transport in a double porosity aquifer are: (a) the movement of the injected CO2, and displacement of water along the fissures (with negligible free-phase CO2 transfer between fissures and matrix). This process is described by the mass balance equation (written in terms of the CO2 free-phase):

@ðf f rc Sc Þ f ¼ rðrc uc Þ f þ q fm @t

(1)

where Sc and rc represent the CO2 saturation and density, respectively, f is the porosity and uc is the Darcy velocity of the CO2 phase. The subscripts f and m stand for the fissure and matrix flow components, and qfm represents the free-phase CO2 transfer between the fissures and the matrix. Notice that qfm implies the existence of Darcian flow of CO2 in the rock matrix, but since the permeability of the rock matrix is very small when compared to the fissures permeability, the term qfm will in practice tend to zero. Significant values of qfm can only be imposed by large injection pressures, comparable to those that would fracture the cap-rock. These processes result in the development of a dry-out zone around the injection well, from which all but residual water is displaced, and in a partially CO2 saturated zone, to which we shall refer as the free-phase zone. (b) the dissolution of the CO2 in contact with water. The phase partitioning behaviour of CO2–H2O mixtures has been described by several equations of state. One such equation is given by Spycher et al. (2003) as:   Fc Ptot ð1  Y w Þ ðP  P0 ÞV c XC ¼ (2) exp  0 RT 55:508KcðgÞ where Xc is the CO2 mass fraction in the groundwater, F is the fugacity coefficient, P is pressure (and the subscripts 0 and tot refer to a reference and total pressure), Yw is the water mass fraction in the carbon dioxide phase, Vc is the average partial molar volume of CO2 and K0 is the gas equilibrium constant.

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

(c) the advection along fissures of the dissolved CO2 enriched water, and the molecular diffusion of dissolved CO2 into the matrix. This process can be described (for slab shape matrix blocks) as: @ðX c Þ f @ð~ vw X c Þ f D @ðX c Þm  ¼ þ þR @t @x a @z z¼a (3) 2 @ðX c Þm @ ðX c Þm ¼ D þR @t @z2 where vw is the water velocity in the fissures, D* is the apparent matrix diffusion coefficient and a refers to a characteristic length of the matrix block (Barker, 1982). R refers to a sink/ source term due to chemical reaction of the dissolved CO2. Matrix diffusion depends mainly on the ratio between fissure porosity and matrix porosity, s = fm/ff, on the characteristic time to diffuse across a matrix block, tcb = a2/D*, where a is a characteristic length, and on the matrix block shape. Hydrodynamic dispersion occurs also along the fissures, but is usually less important than advection. Eqs. (1) and (3) are coupled by the joint variation in CO2 and water saturation, Sc þ Sw ¼ 1, and by phase partitioning behaviour of CO2–H2O, Eq. (2). In porous aquifers (Fig. 3a) the CO2 dissolution and reaction are highly constrained by the groundwater displacement, which creates a relatively narrow contact surface between CO2 and groundwater. Dissolution and reaction will occur mainly at the outer rim of the free-phase zone. What is hypothesized in this paper is that in double porosity aquifers, dissolution and reaction

Fig. 3. Schematic illustration of the physical processes acting on CO2 transport in: (a) porous aquifer; (b) double porosity fissured aquifer. Free-phase CO2 can occur in the matrix closer to the injection well, but according to the double porosity model, matrix permeability is too small to allow for significant advection.

433

may proceed faster than in porous aquifers because the injected CO2 primarily displaces groundwater along fissures, with a large volume of groundwater remaining in the rock matrix. The nearly stagnant matrix water provides a large contact surface between water and CO2, increasing dissolution and reaction (Fig. 3b). Due to the matrix diffusion process and given the time scales involved, all the matrix porosity will eventually become accessible to CO2, providing the volume necessary for carbon storage to be practical (carbonate rocks frequently show matrix porosities between 10% and 20%, and aquifers such as the Chalk in SE England, can show a matrix porosity up to 45% (Price et al., 1993a,b)). 3. Methodology This study compares the efficiency of the CO2 trapping mechanisms by simulating the long-term fate of CO2 injected in single porosity and double porosity aquifers with equivalent hydraulic properties and total porosity. Numerical simulations were conducted using the generic purpose simulator TOUGH2 (Pruess et al., 1999) and the Equation of State (EOS) module applied was ECO2N (Pruess, 2005; Pruess and Spycher, 2007), developed specifically for simulating the behaviour of CO2 injected in a saline aquifers. TOUGH2 includes an option to model multiple porosity media through a process designated as MINC (Multiple Interacting Continua), in which resolution of the physical processes is achieved by appropriate subgridding of the matrix blocks. In the MINC procedure matrix and fractures may exchange fluid (or heat) locally by means of ‘‘interporosity flow,’’ which is driven by the difference in pressures (or temperatures) between matrix and fractures. In order to accurately describe such flows TOUGH2 resolves the driving pressure, temperature, and mass fraction gradients at the matrix/fracture interface and within the matrix (Pruess and Narasimhan, 1982, 1985). The simulated scenarios are conceptual, with no attempt to reproduce experimental data. The simulations are concerned only with the matrix diffusion process, and ignore other features relevant for CO2 injection but that would mask the matrix diffusion effect. That is the case of: (a) chemical reactions between dissolved CO2 and the rock matrix (including sorption); (b) the mechanical effects associated to the injection of a pressurized fluid, and; (c) the changes in porosity and permeability that may be associated with (a) and (b). Still the simulations are representative of the freephase and dissolved phase exchange between fissures and matrix. The layout of the simulations is based on an ECO2N test case (Pruess, 2005), with a CO2 injection well fully penetrating a homogeneous, isotropic, saline aquifer of 50 m thickness (Fig. 4a). The numerical grid is extended to a distance of 20 km from the existing well. Vertical discretization considers 10 layers of uniform thickness with CO2 injection at the lowermost layer. Two flow conditions are considered: (a) radial flow—so that the finite difference grid represented in Fig. 4a is radial and the injection well is at origin of the grid; (b) regional flow—aimed at simulating conditions under which there is groundwater flow before and after the CO2 injection, so that imbibition with groundwater will eventually occur. This was achieved by imposing a constant, uniform, hydraulic gradient from left to right. The CO2 injection well is located at R = 6000 m. The finite difference grid (Fig. 4b) is no longer radial, and instead it represents a vertical cross-section parallel to the regional flow direction. This setup, although not representing the 3-D flow conditions imposed by an injection well, reduces the computational demand, and is representative of the CO2 fate along a flow path that remains directionally invariant during and after the CO2 injection. The simulations consider a fissured double porosity aquifer and its equivalent single porous media, hereafter designated as the Equivalent Porous Media (EPM). The EPM results are used as base

434

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

Fig. 4. Layout of the numerical simulations: (a) radial flow conditions, indicating the main hydraulic parameters for the double porosity aquifer and its Equivalent Porous Media (EPM); (b) regional flow conditions, groundwater flow from left to right. The smaller image shows the alignment represented in the numerical model.

cases against which the double porosity simulations are compared. Hydraulic parameters for the simulations were adopted according to one of the ECO2N examples (Pruess, 2005) and are given in Table 1. EPM equivalence to the fissured aquifer was guaranteed by maintaining the same total porosity (12%) and the same permeability (1013 m2, fissure permeability in the double porosity aquifer). Simulations were conducted for different ratios between matrix and fracture porosity, s, always making sure that the total porosity remains equal to 12%. Cubic, rectangular and slabs matrix block geometry were used, defined by fissure spacing (S) of 1, 5 and 10 m. 3.1. MINC discretization The secondary mesh representing the matrix discretization, according to the MINC approach, was tested to define the number of subgridding levels necessary to accurately represent the influence of the matrix. For this purpose the setup in Fig. 4a was simplified by removing vertical discretization and setting the aquifer thickness (H) to 10 m, so that the effect of ignoring buoyancy is not too unrealistic. The flow pattern is then 1-D radial. Tests were made considering discretization of the matrix into 4, 6, 8 and 10 levels. The results of those tests are illustrated in Fig. 5,

in terms of dissolved CO2 in the rock matrix, and of free-phase CO2 in the fissures. Although differences exists, the results with 4, 6, 8 and 10 MINC levels are very similar, and so it was decided to use four MINC levels, each taking 25% of the matrix volume. Notice that in every case the free-phase CO2 migration into the matrix was zero, which is according to expectations since the matrix permeability is of the same order of magnitude as usually expected for the cap-rock. 4. Results The simpler setup used for studying the MINC subgridding, without vertical discretization, also permits to plot the variation of CO2 components within the rock matrix, which is very difficult to achieve when vertical discretization exists (in which case it is only possible to plot the mean value of each parameter). Such plots are useful in illustrating the matrix diffusion effect on the injected CO2 and are helpful for explaining the CO2 behaviour in more complex simulations. The CO2 saturation profiles (Fig. 6) show that the free-phase zone can be considerably larger in the double porosity environments than in the EPM simulations. This indicate that for short times most CO2 is accommodated only in the fissures. In the limit,

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

435

Table 1 Hydraulic and modeling parameters. Parameter

Equivalent porous media

Double porosity media

Permeability, k (m2) Porosity, f

1013 0.12

Fissures: 1013; matrix: 1016 s = fm/ff = 10–20; ff + fm = 0.12

Pore compressibility, c (Pa1)

4.5  1010 water—Van Genuchten function with l = 0.457 CO2—Corey curve with l = 0.457

Relative permeability, kr Residual saturation, Sr

water - 0.3; CO2 - 0.05

Capillary pressure, Pcap

Van Genuchten function with l = 0.457, P0 = 19.61 kPa

Fracture spacing, S

–

1, 5 and 10 m, cubic, rectangular and slab shape matrix blocks

Temperature, T (8C) Initial pressure, P0 (bar)

45, isothermal simulations Hydrostatic pressure with P0 = 120 bar at the elevation of the injection node 2  1010

Apparent molecular diffusion coefficient, D* (m2/s)

Fig. 5. Sensibility to discretization of the matrix according to the Multiple Interacting Continua (MINC) method. (a) Dissolved CO2 in the matrix; (b) free-phase CO2 in the fissures. Results are very similar and sometimes indistinguishable.

and assuming that reaction and dissolution are irrelevant during the injection period, the free-phase zone radius in the fissured media (Rf) relates to the same radius in an EPM according to: R f ½t ffi Rs

pffiffiffiffiffiffiffiffiffiffiffiffi s þ 1 for small t

(4)

where Rs and s = fm/ff stand for the free-phase radius and the ratio of matrix (fm) and fissure (p fffiffiffiffiffiffiffiffiffiffiffiffi f) porosities, respectively. However, Rf can become lower than Rs s þ 1 due to the diffusion of CO2 into the matrix. For the radial flow conditions in Fig. 6, s = 0.10 and thus one would expect Rf ffi 3.3Rs. In Fig. 6a, we obtain Rf = 3.2Rs at the end of injection period and Rf = 1.1Rs at the end of the simulation, reflecting the increasing importance of CO2 diffusion into the matrix. Notice that unlike in the EPM aquifer in which the freephase radius increases slightly after the injection period, in the double porosity media the free-phase radius decreases with time, eventually becoming of the same order of magnitude as in the EPM situation (Fig. 6b). The dry-out zone also shows significant differences (Fig. 6a). At the end of the injection period the radius of the dry-out zone follows the proportion given by Eq. (4), but as soon as the injection stops the influence of matrix diffusion becomes obvious. In the EPM environment the dry-out zone decreases very slowly since it can only interact with the small volume of residual groundwater in and around the dry-out zone. As for the double porosity, the dryout zone is only in the fissures; a larger volume of groundwater remains trapped in the rock matrix and starts interacting with the

CO2, dissolving it and diffusing it along the matrix blocks. Consequently water starts to move into the dry-out zone, with the CO2 saturation decreasing in time. The mass fraction of dissolved CO2 in the nearly stagnant matrix water is depicted in Fig. 7a, which corresponds to a matrix block (5 m thick) bounded by fissures in the lower (z = 0 m) and upper (z = 5 m) limit. The x-axis is the distance from the injection well. Fig. 7b shows the same section for the EPM simulation. The EPM results (Fig. 7b) are fairly simple, with the water becoming saturated in dissolved CO2 in the area where there is free-phase CO2, and with smaller mass fractions in a very narrow fringe at the contact between the free-phase zone and the displaced groundwater. The mass fraction of dissolved CO2 varies slowly with time since dissolution can only occur in the outer contact between the CO2 free-phase zone and the groundwater. In the rock matrix of the double porosity aquifer (Fig. 7a) as the CO2 displaces the water in the fissures, the groundwater in the matrix gets trapped and remains as isolated pockets of water surrounded by free-phase CO2 in the fissures. Molecular diffusion acts to transport the dissolved CO2 to the inner matrix, allowing for further CO2 to dissolve in the water at the interface. Therefore, in time, the dissolved CO2 in the rock matrix moves toward the centre of the matrix, increasing the concentration in the matrix, until equilibrium is reached across the matrix (visible for the 1000-year simulation in Fig. 7a). Fig. 8 compares the free-phase mass and the dissolved phase mass. The decrease in the free-phase CO2 in the double porosity

436

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

Fig. 6. Saturation profiles for simplified radial flow conditions. (a) CO2 saturations versus distance; (b) free-phase radius variation with time. A CO2 saturation of 0.01 was adopted as the outer limit of the free-phase zone. Double porosity results considering s = 10, slab matrix geometry, and fissures spaced 5 m.

media is remarkable when compared to the EPM, with the freephase CO2 almost disappearing at the end of the simulation, as an effect of the enhanced dissolution of CO2 caused by matrix diffusion (Fig. 8b). The total mass of dissolved CO2 is not only much bigger in the double porosity scenario (Fig. 8a), but also the rate of dissolution is much faster. 4.1. Radial flow conditions The above 1-D radial simulations, although with limited validity to represent a real CO2 injection scenario since it ignores buoyancy, set the framework for understanding the relevance of matrix diffusion in carbon storage and its importance in a more realistic setup can be studied with the model layout illustrated in Fig. 4. Fig. 9 shows the free-phase distribution for the EPM and for fissures of the double porosity model under radial flow conditions. Notice that the horizontal scales in Fig. 9a and b are different, reflecting the need for the double porosity media to take a larger aquifer volume to accommodate the same amount of CO2. The relation in Eq. (4) does hold during the injection period, but while the plume remains relatively stable in the double porosity environment, it continues to grow in the EPM, decreasing the ratio between Rf and Rs. Buoyancy of the free-phase is faster in the double porosity media than in the EPM, reflecting the larger flow

velocity imposed by the small fissure porosity along which the free-phase moves. Consequently, the CO2 free-phase plume in the double porosity will reach the top of the aquifer (base of the caprock) rapidly and will then move predominantly in the horizontal, reaching faster its maximum spread. In the EPM, and because vertical flow is slower, the plume will tend to spread in the horizontal for longer times. For the case illustrated in Fig. 9 the CO2 free-phase has attained its maximum spread at the end of the injection period (ffi8 km radius), while in the EPM it continues to spread beyond the end of the injection (ffi2.2 km radius) and up to the end of the simulation period (ffi3.6 km). Fig. 10 compares the free-phase mass and the dissolved phase mass for several double porosity scenarios. In every case, even for relatively small s ratios and block sizes, the decrease in the freephase CO2 in the double porosity media is remarkable when compared to the EPM. For the case illustrated in Fig. 9 (s = 10, cubes, S = 5 m), 100 years after the beginning of injection more than 41% of the total injected CO2 had already dissolved into the groundwater, against only 10% in the EPM. In the case with a large s ratio (s = 20, not uncommon in double porosity environments), that percentage raises to 63%. Therefore, even considering a buoyancy effect, which reduces the volume of matrix in contact with CO2, matrix diffusion can play an important role in the longterm fate of CO2 in fissured aquifers.

Fig. 7. Mass fraction of dissolved CO2 (X CO2 aq ) for simplified radial flow conditions: (a) matrix of double porosity aquifer; (b) EPM. Double porosity results considering s = 10, slab matrix geometry and fissures spaced 5 m.

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

437

Fig. 8. Variation with time of: (a) dissolved CO2 (b) and free-phase CO2 with time, for simplified radial flow conditions. Double porosity results considering s = 10, slab matrix geometry and fissures spaced 5 m.

Fig. 11 shows the dissolved CO2 profiles 100 years from the start of injection, for the fissures and for the matrix (mean value across the matrix) when compared to the EPM situation (Fig. 11a). Notice that although the mean mass fraction of CO2 dissolved in the matrix is lower than that in the EPM, in fact a larger mass of CO2 has dissolved in the matrix (3.1 Mton, see Fig. 10), since the matrix is completely saturated with water into which the CO2 can dissolve. On the contrary, in the EPM the larger mass fraction of dissolved CO2 occurs in an area where there is partial saturation with water, since the free-phase CO2 exists in almost the same area (compare with Fig. 9a), and thus only 0. 9 Mton have been dissolved in the

EPM situation. In the fissures 0.3 Mton of CO2 have been dissolved, which is comparatively more favourable than in the EPM which has a porosity of 12% against a fissure porosity of only 1%. 4.2. Regional flow conditions In order to test the influence of re-imbibition of the fissures and matrix with water, a set of simulations were conducted in which a regional flow was imposed, so that there is always an advective movement. The free-phase CO2 is displaced and the groundwater moves back into the fissures and matrix. Obviously, some of the

Fig. 9. CO2 saturation (Sc) profiles for radial flow conditions: (a) EPM; (b) fissures of double porosity, for 30 years (top images) and 100 years (bottom images). Double porosity results considering s = 10, cubic matrix geometry, and fissures spaced 5 m. Scale varies.

438

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

Fig. 10. Variation with time of: (a) dissolved CO2 and (b) free-phase CO2 for radial flow conditions.

injected CO2 will remain as a residual component, but most of it will be displaced by the imbibition water. The regional flow was forced from left to right of the model domain (see Fig. 4b) by imposing Dirichlet boundary conditions at the cells in the left and right limits of the domain, with the pressures set as shown in Fig. 4b. The double porosity simulations results presented for this scenario refer only to the situation when s = 10, with cubic matrix geometry defined by three sets of orthogonal fissures spaced 5 m. Figs. 12 and 13 show the CO2 saturation and dissolved CO2 profiles, respectively, at the end of the injection period (1 year) and after 500 years. The influence of the regional flow is more obvious in the double porosity model due to the larger advective velocity,

which spreads the plume faster down-gradient from the injection well. The variation of the total mass of free-phase CO2 (Fig. 14b), shows very similar behaviour between the EPM and the double porosity media up to 40 years, after which the free-phase CO2 in the double porosity media starts decreasing at a faster rate than in the EPM, so that at the end of the simulation the free-phase CO2 in the double porosity aquifer is only 14% of the total injected, while in the EPM 60% of the CO2 remains in free-phase. Compared to the situation without regional flow (see Fig. 10), the most striking is that the dissolution of CO2 proceeds initially at a similar rate in the EPM and in the double porosity media. Tests were made by applying a more intensive discretization of the

Fig. 11. Mass fraction of dissolved CO2 (X CO2 aq ) for radial flow conditions. (a) EPM; (b) double porosity (top – fissures; bottom – matrix). Double porosity results considering s = 10, cubic matrix geometry, and fissures spaced 5 m. Scale varies.

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

439

Fig. 12. CO2 saturation (Sc) profiles for regional flow conditions: (a) EPM; (b) fissure of double porosity, for 1 year (top images) and 500 years (bottom images). Double porosity results considering s = 10, cubic matrix geometry, and fissures spaced 5 m. Scale varies.

matrix (10 levels) to check whether it could be a result of the MINC discretization, but the same behaviour was noticed. That behaviour is possibly due to enhanced dissolution of CO2 by the incoming groundwater, and that effect is only overcome by the influence of matrix diffusion after a sufficiently large volume of matrix is in contact with the CO2. That is, for short times, the dominant factor for promoting CO2 dissolution is the regional flux. However, for longer times the matrix diffusion has a dominant influence, and this is further enhanced by the regional flow. As the free-phase and dissolved CO2 are moving along fissures, CO2 starts diffusing out of the matrix and into the fissures as soon as these are filled with incoming groundwater. Downstream, as the displaced groundwater moves along the fissures, the CO2 diffuses again into the matrix, and thus dissolution proceeds at a faster rate than in the situation without regional flow.

discussion, since in fissured rocks it is likely to be a function of the roughness of the fissures, and thus it may be quite difficult to assess using the double porosity concept in which fissures are seen as smooth surfaces. 5.1. Stratigraphic/structural trapping Stratigraphic and structural trapping refers to the process by which a fluid in gas, liquid or supercritical phase is contained in a static position beneath impermeable layers (Bachu et al., 2007). The traps include anticlines, faults and structural and stratigraphic pinch-outs. Bachu et al. (2007) indicates the following formulae for quantifying the theoretical storage volume (V C;t ) due to stratigraphic and structural trapping: V C;t ¼ Ahfð1  Swirr Þ

(5)

5. Discussion The CO2 trapping mechanisms active in a double porosity fissured media are not different from those in a porous media: (i) stratigraphic/structural trapping; (ii) dissolution trapping; (iii) mineral trapping; (iv) residual trapping; and (v) hydrodynamic trapping (Bachu et al., 2007; Bradshaw et al., 2007; Chadwick et al., 2008). The results of the simulations described in the previous sections make it possible to establish some considerations on the relevance of those trapping mechanism in double porosity reservoirs. Residual trapping is not addressed in this

where A and h are the trap area and average thickness, respectively, f is the effective porosity and Swirr is the irreducible water saturation. For any given fixed area trap, the volume of CO2 that can be stored is proportional to the porosity, implying that if fissure porosity was the only one accessible to CO2, the storage volume would be necessarily small. The free-phase zone radius as studied from the gas saturation profiles show that for short times (say, the injection period), and as expected, the free-phase zone radius in the double porosity environment will be larger than in the EPM, since at short times

440

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

Fig. 13. Mass fraction of dissolved CO2 (X CO2 aq ) for regional flow conditions. (a) EPM; (b) double porosity(top – fissures; bottom – matrix). Double porosity results considering s = 10, cubic matrix geometry, and fissures spaced 5 m. Scale varies.

only the fissure porosity is available for storage. For later times, and as the CO2 dissolves and diffuses into the rock matrix, the freephase zone in the double porosity reservoir will remain stationary or will tend to decrease faster than in single porosity media. For relatively thick aquifers, in which buoyancy effects play an important role, and because flow velocities are higher in double porosity aquifers, the maximum extent of the free-phase plume is reached sooner than in single porosity aquifers, where the plume may grow in extent for longer periods. The larger area taken by the free-phase CO2 in double porosity aquifers is associated with a higher risk of intercepting imperfec-

tions in the cap-rock that may lead to leakage of CO2. Therefore, because the structural and stratigraphic trapping is relevant mainly during the injection period, it is obvious that this trapping mechanism is less efficient in double porosity aquifers than in single porosity aquifers. 5.2. Dissolution trapping The dissolution trapping occurs when the CO2 moving through the reservoir porous space dissolves into the groundwater which it comes into contact with. The CO2 will dissolve into the

Fig. 14. Variation with time of: (a) dissolved CO2 and (b) free-phase CO2 for the regional flow conditions. Double porosity results considering s = 10, cubic matrix geometry, and fissures spaced 5 m.

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

441

groundwater fully or partially, depending on time, CO2 saturation rate of the water, water chemistry and the rate of contact between CO2 and water. Bachu et al. (2007) indicates the following formulae for quantifying the theoretical storage mass (M C;t ) due to dissolution trapping: M C;t ¼ Ahfðrs X C;s  r0 X C Þ

(6)

where r and XC are the density and mass fraction on the groundwater, respectively, and the subscripts s and 0 stand for initial CO2 content and CO2 content at saturation. As seen from the numerical simulations, the solubility trapping in double porosity fissured aquifers is much more efficient than in single porosity environments, because the effective volume of water (i.e., Ahf in Eq. (6)) is much larger in double porosity environments since the water in the rock matrix will not be displaced (at least not to high degree) and will provide a larger volume of water into which the CO2 can diffuse and dissolve. Figs. 10a and 14a are particularly elucidative of this effect, with the mass of dissolved CO2 increasing faster in the double porosity environment than in the EPM. Thus, the matrix diffusion process plays a major role in the evolution of the CO2 plume and is a very effective process for definitive sequestration of the CO2. 5.3. Mineral trapping

2

> 0:001

 RSCO

2

(7)

< 0:5 Þ

where RX CO > 0:001 is the radius in which the dissolved CO2 2 concentration is above 0.001 and RSCO < 0:5 is the radius for which 2 the free-phase CO2 saturation is below 0.5. X CO2 > 0:001 and SCO2 < 0:5 , were fixed arbitrarily, in the first case small enough to guarantee that all relevant amount of dissolved CO2 is accounted for, and in the later case, to ensure that the water phase is the predominant one. As for the double porosity simulations, the approximation used was: V c ¼ 2ph½f f ðR f ;xCO

2

Additionally because the flow rate in the matrix is extremely slow, nearly stagnant, there will be a longer contact time during which reactions may occur. Therefore, although evaluated in an indirect manner, it is likely that mineral trapping in double porosity environments has an increased importance when compared to the same process in porous media. 5.4. Hydrodynamic trapping

Mineral trapping refers to the mechanism by which the CO2 reacts with the host rock and the groundwater and precipitates in mineral form. The scale of occurrence is of hundreds to thousands of years and it depends strongly on the mineralogy and chemistry of the groundwater. Other controlling factors are the pressure and temperature, the contact surface between the mineral grains and the groundwater containing dissolved CO2, and the flow rate of the fluids past that contact surface (Bachu et al., 2007). The code used for the numerical simulations does not have reaction modelling capabilities, nor the possibility to compute the contact volume between rock and the dissolved CO2. However, the aquifer volume, under radial flow conditions, with dissolved CO2 was approximated for the EPM, by: V c ¼ 2phfðRX CO

Fig. 15. Volume of aquifer with dissolved CO2, for the simplified flow conditions (buoyancy effects ignored).

> 0:001

 R f ;SCO

2

< 0:5 Þ

þ fm Rm;xCO

2

> 0:001 

(8)

where the first term in the square brackets accounts for the volume in the fissures and the second term accounts for the volume in the matrix. Notice that in the case of the matrix there is no need to consider an upper limit to the CO2 saturation, since there is no freephase CO2. The evaluation of Vc was made for the simpler 1-D radial simulations (illustrated in Fig. 7), no buoyancy effects considered, for several double porosity scenarios (Fig. 15). At the end of the injection period, the volume in contact between rock and CO2 is already higher than that in the EPM, due to the diffusion of dissolved CO2 into the rock matrix, and it continues to increase for hundreds of years, suggesting that the potential for reaction is much higher than in the EPM.

Hydrodynamic trapping exists when CO2 is injected in strata in which there is flow imposed by a regional hydraulic gradient, so that the groundwater migrates in long, regional and basin scale flow systems. Even if no structural or stratigraphic traps prevent the CO2 from moving laterally, the CO2 migrates together with the groundwater at such low velocities that it would take thousands of years to reach the discharge areas (Bachu et al., 2007). The numerical simulations considering a regional flow show that the free-phase CO2 decreases rapidly in the double porosity domain (Fig. 14), eventually vanishing. This process is not different from that observed without the regional flow, and again it results from the dissolution and diffusion into the stagnant matrix water. However, when the CO2 is moving together with the groundwater, matrix diffusion plays an even more important role, since the freephase is permanently becoming into contact with matrix water into which no CO2 has yet dissolved. The possibility to dissolve and diffuse in the matrix is renewed continuously, leading to a faster disappearance of the free-phase CO2. Enhanced CO2 dissolution also occurs in single porosity media, since the regional flow promotes the contact between the free-phase CO2 and the groundwater. However, the dissolution rate is slower than in the double porosity media, given that in the EPM the free-phase CO2 is largely isolated from contact with the flowing water. Hence, hydrodynamic trapping may be a more effective trapping mechanism in double porosity environments than in single porosity reservoirs, due to the enhanced dissolution. However, the larger flow velocity in fissured media will cause the free-phase CO2 to move faster down-gradient, which can be a serious disadvantage due to the increased risk of encountering caprock imperfections. 6. Conclusions The double porosity conceptualization of fissured aquifers and reservoirs is a well established model in the hydrogeology and oil engineering fields, with many theoretical, laboratory and field studies confirming its reliability to describe fluid flow and transport in fissured rocks. Nevertheless, in what regards carbon storage in geological media, the properties of the double porosity fissured media have usually not been considered. And yet, since the matrix

442

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443

porosity can be very high (for instance, the Chalk aquifer in England can show matrix porosities of up to 45% (Price et al., 1993a,b)), double porosity aquifers may be suitable formations for CO2 storage. The numerical simulations presented in this paper show that the matrix diffusion process may play a relevant part in carbon sequestration in fissured aquifers. According to the simulations dissolution trapping is more important in double porosity fissured aquifers, because water in the rock matrix is not displaced (at least not in a large volume) by the injected CO2, which moves along fissures and only displaces the water on them. This situation leads to enhanced dissolution of CO2 at the fissures/matrix interface, and to the diffusion of the dissolved CO2 across the matrix block. The amount of CO2 that is able to be trapped by dissolution in this manner is much larger then expected in porous media. Additionally, the matrix diffusion process also allows for a larger contact surface and longer contact time between dissolved CO2 and the rock minerals promoting a more favourable environment for chemical reactions. Even if the simulations do not account directly for mineral reactions, the indirect evidences is that the mineral trapping in double porosity aquifers will also be more important than in porous aquifers. Hydrodynamic trapping was addressed in the simulations, and again it proved a more effective means for the sequestration of the CO2 in double porosity reservoirs than in the equivalent single porosity media, with the matrix diffusion being enhanced by the CO2 movement along the fissures, where it is permanently coming into contact with water in the rock matrix that is not yet saturated with dissolved CO2. However, faster movement would imply that the free-phase CO2 could reach discharge zones sooner. Therefore, hydrodynamic trapping will only be more favourable in double porosity media depending on the balance between rate of dissolution and the transport velocity. Structural and stratigraphic trapping play a distinct role, being less effective in double porosity systems because only the fissures are available for storage of the free-phase CO2 which will spread through a larger area than in porous aquifers. If the storage strategy was aiming at keeping the free-phase CO2 under structures such as anticlines or domes, a smaller storage volume would be possible until the spill-out limit had been reached. Furthermore, this need for larger aquifer volumes to accommodate a similar amount of CO2 results in an increased risk of encountering imperfections in the cap-rock through which the CO2 could leak. Thus, as far as structural/stratigraphic effectiveness and safety, double porosity aquifers are less suitable for carbons storage than equivalent single porosity media. Several other issues need to be addressed in order to have a full picture of the suitability of fissured saline aquifers for storing carbon dioxide. Apart from the obvious need to conduct laboratory and field works to which the numerical simulations can be confronted, two issues are particularly relevant and should be addressed:  The influence of the mineral reaction in the matrix diffusion and on the permeability of fissures; carbonate rocks such as limestones and dolomites show double porosity features and injection of CO2 into those rocks could lead to deposition of calcite, decreasing both porosity and permeability. Laboratory and numerical investigations have been conducted that study the reaction of CO2 injected in carbonate rocks (Izgec et al., 2008a,b; Madland et al., 2006), but not accounting for double porosity behaviour nor the matrix diffusion effect;  Channelling, i.e., transport of solutes and components along preferential fissures, is a well-known phenomena that can be responsible for very high advective velocities (e.g., Berkowitz and Braester, 1991; Moreno and Neretneiks, 1993; Tsang, 1987). The occurrence of channelling would obviously be problematic for

CO2 storage in open basins, and the studies addressing this factor should be conducted. Much knowledge about this subject has been collected in contaminant hydrogeology and in studies for radioactive waste repositories, but there is a need to adapt the conclusions of those studies to the particular subject of CO2 storage. In any case, the matrix diffusion effect of double porosity reservoirs may have interesting consequences for the definitive sequestration of CO2 in fissured aquifers, mainly due to the enhanced dissolution, mineral and hydrodynamic trapping associated to the mass transfer between fissures and rock matrix. However, the larger aquifer volumes necessary to accommodate the free-phase CO2 and the higher flow velocities in double porosity aquifers imply increased risks of intercepting imperfections in the cap-rock and of reaching discharge zones, leading to CO2 leakage. In that respect, and despite matrix diffusion promoting a long-term safe sequestration of CO2, double porosity aquifers should be considered as potential carbon storage reservoirs only if robust characterization of the cap-rock is possible and the existing flow and discharge conditions are well known. References Bachu, S., 2003. Screening and ranking of sedimentary basins for sequestration of CO2 in geological media in response to climate change. Environmental Geology 44 (3), 277–289. Bachu, S., Bonijoly, D., Bradshaw, J., Burruss, R., Holloway, S., Christensen, N.P., Mathiassen, O.M., 2007. CO2 storage capacity estimation: methodology and gaps. International Journal of Greenhouse Gas Control 1 (4), 430–443. Barenblatt, G.E., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. Journal of Applied Mathematics and Mechanics 24 (5), 1286–1303. Barker, J.A., 1982. Laplace transform solutions for solute transport in fissured aquifers. Advances in Water Resources 5, 98–104. Barker, J.A., 1985a. Generalised well-function evaluation for homogeneous and fissured aquifers. Journal of Hydrology 76, 143–154. Barker, J.A., 1985b. Modeling the effects of matrix diffusion on transport in densely fissured media. In: Day, J.B. (Ed.), Hydrogeology in the Service of Man—18th Congress of the Intl. Assoc. Hydrogeologists. IAHS Publication, Cambridge, UK, pp. 250–269. Barker, J.A., Wright, T.E.J., Fretwell, B.A., 2000. A pulsed-velocity method of doubleporosity solute transport modelling. In: Dassargues, A. (Ed.), Tracers and Modelling in Hydrogeology. IAHS Pub. 262, Liege, Belgium. Bentham, M., Kirby, G., 2005. CO2 storage in saline aquifers. Oil & Gas Science and Technology-Revue de L’Institut Francais du Petrole 60 (3), 559–567. Berkowitz, B., Braester, C., 1991. Solute transport in fracture channel and parallel plate models. Geophysical Research Letters 18 (2), 227–230. Bradshaw, J., Bachu, S., Bonijoly, D., Burruss, R., Holloway, S., Christensen, N.P., Mathiassen, O.M., 2007. CO2 storage capacity estimation: issues and development of standards. International Journal of Greenhouse Gas Control 1 (1), 62–68. Carrera, J., Sanchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guimera, J., 1998. On matrix diffusion: formulations, solution methods and qualitative effects. Hydrogeology Journal 6 (1), 178–190. Castello, M., Bouchi-Lamontagne, M., 2007. CO2 capture and storage in the subsurface. In: A Technological Pathway for Combating Climate Change. Geoscience Issues, Ifp, Ademe, Brgm, Rueil-Malmaison, 64 pp. Celia, M.A., Nordbotten, J.M., Gasda, S.E., Kavetski, D., Bachu, S., 2007. Geological storage as a carbon mitigation option. Geochimica et Cosmochimica Acta 71 (15), A153–A1153. Chadwick, R.A., Arts, R., Bernstone, C., May, F., Thibeau, S., Zweigel, P. (Eds.), 2008. Best Practice for the Storage of CO2 in Saline Aquifers: Observations and Guidelines from the SACS and CO2STORE Projects, vol. iv. British Geological Survey Occasional Publication. British Geological Survey, Keyworth, pp. 267. Di Donato, G., Blunt, M.J., 2004. Streamline-based dual-porosity simulation of reactive transport and flow in fractured reservoirs. Water Resources Research 40 (4). Hadermann, J., Ro¨sel, F., 1985. Radionuclide chain transport in inhomogeneous crystalline rocks: limited matrix diffusion and effective surface sorption. Nagra Technical Report, NTB 85-40. IPCC, 2005. IPCC Special Report on Carbon Dioxide Capture and Storage, Vol. x. Cambridge University Press, Cambridge, 431 pp. Izgec, O., Demiral, B., Bertin, H., Akin, S., 2008a. CO2 injection into saline carbonate aquifer formations I: laboratory investigation. Transport in Porous Media 72 (1), 1–24. Izgec, O., Demiral, B., Bertin, H., Akin, S., 2008b. CO2 injection into saline carbonate aquifer formations II: comparison of numerical simulations to experiments. Transport in Porous Media 73 (1), 57–74.

J.F. Carneiro / International Journal of Greenhouse Gas Control 3 (2009) 431–443 Kazemi, H., Gilman, J.R., 1993. Multiphase flow in fractured petroleum reservoirs. In: Bear, J., Tsang, C.F., de Marsily, G. (Eds.), Flow and Contaminant Transport in Fractured Rock. Academic Press, Inc., San Diego, pp. 267– 323. Kolditz, O., 1995. Modeling flow and heat-transfer in fractured rocks—dimensional effect of matrix heat diffusion. Geothermics 24 (3), 421–437. Kruseman, G.P., de Ridder, N.A., 1994. Analysis and Evaluation of Pumping Test Data. International Institute for Land Reclamation and Improvement, Wageningen, The Netherlands, 377 pp. Le Gallo, Y., 2007. Field case Casablanca—reservoir model and regional model for simulation of CO2 injection and leakage scenario. In: CO2NET Seminar 2008, Warsaw, Poland. Lever, D.A., Bradbury, M.H., 1985. Rock matrix diffusion and its implications for radionuclide migration. Mineralogical Magazine 49 (351), 245–254. Lokhorst, A., Wildenborg, I., 2005. Introduction on CO2 geological storage. Classification of storage options. Oil & Gas Science and Technology-Revue De L Institut Francais Du Petrole 60 (3), 513–515. Lombardi, S., Altunina, L.K., Beaubien, S.E., 2006. Advances in the geological storage of carbon dioxide: international approaches to reduce anthropogenic greenhouse gas emissions. In: NATO Science Series. Series IV Earth and Environmental Sciences, Vol. xiv, Springer, Dordrecht [Great Britain], 362 pp. Madland, M.V., Finsnes, A., Alkafadgi, A., Risnes, R., Austad, T., 2006. The influence of CO2 gas and carbonate water on the mechanical stability of chalk. Journal of Petroleum Science and Engineering 51 (3–4), 149–168. Maloszewski, P., Zuber, A., 1985. On the theory of tracer experiments in fissured rocks with a porous matrix. Journal of Hydrology 79, 333–358. Maloszewski, P., Zuber, A., 1993. Tracer experiments in fractured rocks: matrix diffusion and the validity of models. Water Resources Research 29 (8), 2723– 2735. Moench, A.F., 1984. Double-porosity models for a fissured groundwater reservoir with fracture skin. Water Resources Research 20, 831–846. Moreno, L., Neretneiks, I., 1993. Fluid flow and solute transport in a network of channels. Journal of Contaminant Hydrology 14, 163–192.

443

Price, M., Downing, R.A., Edmunds, W.M., 1993a. The Chalk as an aquifer. In: Downing, R.A., Price, M., Jones, G.P. (Eds.), The Hydrogeology of the Chalk of North West Europe. Oxford Scientific Publications, Chapter 3. Price, M., Jones, G.P., Downing, R.A., 1993b. The Hydrogeology of the Chalk of NorthWest Europe, Vol. viii. Clarendon Press, Oxford, 300 pp. Pruess, K., 2005. ECO2N: a TOUGH2 fluid property module for mixtures of water, NaCl, and CO2. LBNL-57952. Lawrence Berkeley National Laboratory, Berkeley, CA. Pruess, K., Narasimhan, T.N., 1982. On fluid reserves and the production of superheated steam from fractured, vapor-dominated geothermal-reservoirs. Journal of Geophysical Research 87 (Nb11), 9329–9339. Pruess, K., Narasimhan, T.N., 1985. A practical method for modeling fluid and heat-flow in fractured porous-media. Society of Petroleum Engineers Journal 25 (1), 14–26. Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User’s Guide, Version 2. 0. LBNL-43134. Lawrence Berkeley National Laboratory, Berkeley, CA. Pruess, K., Spycher, N., 2007. ECO2N—a fluid property module for the TOUGH2 code for studies of CO2 storage in saline aquifers. Energy Conversion and Management 48 (6), 1761–1767. Sarma, P., Aziz, K., 2006. New transfer functions for simulation of naturally fractured reservoirs with dual-porosity models. SPE Journal 11 (3), 328–340. Shi, J.Q., Durucan, S., 2005. CO2 storage in deep unminable coal seams. Oil & Gas Science and Technology-Revue de L’Institut Francais du Petrole 60 (3), 547–558. Spycher, N., Pruess, K., Ennis-King, J., 2003. CO2–H2O mixtures in the geological sequestration of CO2 I. Assessment and calculation of mutual solubilities from 12 to 100 degrees C and up to 600 bar. Geochimica et Cosmochimica Acta 67 (16), 3015–3031. Sudicky, E., Frind, E., 1982. Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures. Water Resources Research 18 (6), 1634–1642. Tsang, Y.W.a.T.C.-F., 1987. Channel model of flow through fractured media. Water Resources Research 23, 467–479. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. Society of Petroleum Engineers Journal 3, 245–255.