CO2 flow through a fractured rock volume: Insights from field data, 3D fractures representation and fluid flow modeling

International Journal of Greenhouse Gas Control 18 (2013) 183–199

Contents lists available at ScienceDirect

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

Review

CO2 flow through a fractured rock volume: Insights from field data, 3D fractures representation and fluid flow modeling S. Bigi a,∗ , M. Battaglia a , A. Alemanni b , S. Lombardi a , A. Campana c , E. Borisova c , M. Loizzo c a

Department of Earth Sciences, Sapienza – University of Rome, Italy Upstream Gas, Enel Trade Spa, Rome, Italy c Etudes et Productions Schlumberger, Clamart, France b

a r t i c l e

i n f o

Article history: Received 7 September 2012 Received in revised form 17 July 2013 Accepted 19 July 2013 Available online 15 August 2013 Keywords: Fluid flow Fracture network Numerical modeling Field work

a b s t r a c t Carbon capture and storage (CCS) projects require an accurate evaluation of the sealing potential of faults and highly fractured zones to minimize the potential for CO2 leakage. A study on the control exerted by fracture and fault networks on fluid flow, and in particular on CO2 leakage, should be based upon a representation of discrete fracture networks (DFN) that is as close as possible to that observed in the field. The present research integrates field work analysis, digital data representation, and fluid flow modeling to build a discrete fracture network (defined here as an Analogue Model; AM) that preserves the field-observed fracture geometry and relative proportion of the associated petrophysical parameters (aperture, length, direction and dip). Our study area is an outcrop in the caldera of Latera (Central Italy) where CO2 is naturally released and gas discharge in the atmosphere can be directly observed. We then compare the AM results to those generated by inputting the same fracture parameters into commercial DFN models. Our results highlight that these latter generally underestimate permeability values (by about two orders of magnitude) and hide fault-related anisotropies observed in the field, which instead are very well defined by the AM. The models were applied to a study site in the Latera caldera (Central Italy), where geologically produced CO2 leaks to the atmosphere along an exposed fault, and to simulate gas release through a fractured reservoir. Simulated leakage correlates well with field measurements that show CO2 spot anomalies at fault and fracture intersections and indicate how gas migration pathways are controlled by discontinuity permeability, complex fracture orientations, and fracture positions within the analyzed rock volume. © 2013 Elsevier Ltd. All rights reserved.

Contents 1. 2. 3.

4. 5.

6. 7.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The Latera caldera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Outcrop description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Fieldwork analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Development of the Analogue Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discrete Fracture Network (DFN) models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Upscaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. Analogue Model upscaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. DFNs upscaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3. Comparison of the upscaling process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical models of CO2 leakage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

184 184 185 186 187 188 188 190 190 192 192 193 196

∗ Corresponding author at: Department of Earth Sciences, Sapienza – University of Rome, P.le A. Moro 5, 00185 Roma, Italy. Tel.: +39 06 4991 4992; fax: +39 06 445 4729. E-mail addresses: [email protected], [email protected] (S. Bigi). 1750-5836/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijggc.2013.07.011

184

8.

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix A. Supplementary data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1. Introduction Field and laboratory evidence suggests that we can rarely model flow and transport in a fractured reservoir by treating it as a uniform or mildly non-uniform isotropic continuum. Instead, we must generally account for the highly heterogeneous, directional dependence, dual or multicomponent nature, and multi-scale behavior of fractured rocks. One can either depict the reservoir as a network of discrete fractures (with permeable or impermeable matrix blocks) or as a stochastic (single, dual or multiple) continuum. Modeling discrete fractures follows two basic approaches: probabilistic and deterministic. 3D probabilistic models of fracture networks are based on field data used as input parameters (i.e. orientation, spatial distribution, and density), whereas deterministic models are based on kinematics and dynamic constraints (e.g. van Dijk et al., 2000; Masini et al., 2010). However, accuracy may be limited by the high number of parameters controlling fracturing processes (lithology, stresses, etc.), statistical bias, and the inherent difficulty in accurately estimating uncertainty. Instead, the probabilistic approach creates a statistical depiction of a fractured rock volume, which results in a continuous medium representation of the petrophysical parameters. A third, hybrid approach is to create a stochastic continuous medium that contains a relatively small number of discrete features (see Neuman, 2005, and references therein). In this work, we simulate the 3D fluid flow through a fracture network by integrating digitized field data and numerical modeling to simulate CO2 flow within a fractured medium. Similar procedures have been proposed by several authors, with the difference being that fluid flow simulation in this study used a discrete fracture network (DFN) generated by statistical elaboration of field data (e.g. Gilmour et al., 1986; Barton and Zoback, 1992; Walsh et al., 1998; Caine and Forster, 1999). Data from outcrops representative of reservoir rocks (similar lithology, age, tectonic evolution) can provide information on: (a) fracture size, aperture, length, density and connectivity; (b) the opening and propagation mode of fractures and the control exerted on fracture distribution; (c) fault zone and fracture network architecture and porosity/permeability distribution; and (d) fracture chronology to reconstruct the deformation history. We test our methodology on a fractured outcrop within the Latera caldera (Central Italy), a volcanic area where CO2 and other gases leak naturally. This site has been selected a natural field laboratory for the study of natural CO2 release by several European research projects focused on Carbon Geological Storage. The chosen outcrop, within a small abandoned quarry, is in one of the areas with maximum gas emission. The outcrop includes a main fault plane and its associated damage zone for a total volume of about 26 m3 (3.5 m × 2.5 m × 3 m). The size, orientation, and position of each single fracture plane were measured in the field. 3D representation of the fracture network is built by digitizing the field data AutoCAD (Autodesk, 2012). The use of AutoCAD for the construction of a fracture network in a 3D representation is a relatively new approach, although already proposed (see La Bella, 2006; van Dijk, 2002, 2013). The obtained Analogue Model of the Latera quarry (based on unrefined field data only) is then compared against Discrete Fracture Network models generated by commercial software used in the oil and gas industry (Petrel 2009 by Schlumberger and

198 199 199 199

Move 2011 by Midland Valley). Finally, our 3D Analogue Model is used to simulate the leakage of CO2 using the software COMSOL 4.2 Multiphysics (by COMSOL AB) and the numerical model is compared against field maps of gas emission in the area. Notwithstanding its limitations (e.g. outcrops are exposed to weathering and are not in the same physical condition of a reservoir), the Analogue Model represents a hybrid depiction of a fractured rock volume, comparable to Discrete Fracture Networks, build with a stochastic approach, and fluid flow simulation results are largely comparable to the result of soil gas survey data from the same area. Although we stress the major role played by field work and manual data analysis in our work, we make use of different software and modules to create the DFNs and run the numerical experiments: (a) Move 2011, Module fracture modeling, by Midland Valley; (b) Petrel 2009 with Module fracture modeling, by Schlumberger and Golder Associates; (c) COMSOL Multiphysics with Earth Science module by COMSOL.

2. The Latera caldera The Latera geothermal field, close to the Tyrrhenian coast of Italy, is characterized by a high geothermal gradient and an extensional volcano-tectonic regime (Fig. 1). The Latera caldera is the result of three main collapse phases. First, the formation of an old caldera, now buried, related to the eruption of ignimbrites from the magma chamber (located about 2 km below the present-day surface); second, the sinking of the eastern sector due to the formation of the nearby Bolsena caldera; and third, the multistage formation of the present caldera (∼0.16 My Palladino and Simei, 2005 and references therein). The area is a good natural test site to study CO2 flow in fractured media, and to verify gas sampling methods and monitoring techniques. The production of CO2 is due to thermo-metamorphic reactions at depth. CO2 is the carrier for trace (i.e. He) and other minor gas species (i.e. H2 S, CH4 and H2 ) (Fig. 2). This geological setting allows to measure the fluid flow from fractures systems and fault planes, using soil gas and CO2 flux surveys (e.g. Annunziatellis et al., 2008; Pettinelli et al., 2010). Good exposures are available in abandoned quarries of kaolinite and sulfur, which are products of alteration of the lava deposits (Fig. 3). The young age of the volcanic deposits (0.3–0.1 My), crossed by fractures and faults, allows us to refer them to a recent, single extensional volcano-tectonic event in the area (Palladino and Simei, 2005 and references therein). Because of this, exposed outcrops may be considered as qualitatively representative of the geometrical distribution of fractures at depth. The main regional faults and fractures trend N10◦ E and N50◦ E, suggesting the occurrence of a normal fault system connected to the latest extensional tectonics of the Tyrrhenian margin (Figs. 1 and 2; Zito et al., 2003; Mazzoli et al., 2012 and reference therein). The correlation between the alignments of gas vents and the main normal faults suggests that these faults control the gas emission in the area (Annunziatellis et al., 2008 and references therein). The local (outcrop scale) density of fractures and faults can vary dramatically from about 1 fracture every 150 m to about one fracture every 0.5–0.4 cm. Fractures and faults cluster around the two main directions of N50◦ E and N10◦ E, even at the local scale (Fig. 2).

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

185

Fig. 1. Geological sketch of the Italian peninsula. The caldera of Latera is close to the Tyrrhenian coast (indicated by a box), an area characterized by significant volcanic activity and a high geothermal gradient.

2.1. Outcrop description The studied outcrop is located in the northwestern sector of the caldera, in an abandoned kaolinite and sulfur quarry. Exposed lithologies are lavas and pyroclastic products related to the caldera volcanic activity. Different fault zone architectures, observed throughout the caldera, are well exposed in the quarry. In this extensional volcano-tectonic regime, fault architectures represent different stages of the faulting process (i.e. from younger to mature) (Figs. 2 and 3). Younger Fault Zones are characterized by thick damage zones (see Caine et al., 1996) that are formed by conjugate systems of normal faults and shear fractures, and a thin (often absent) core zone. Damage zones are 4–5 m wide and are repeated every 8–10 m. Cross cutting relationships among fractures and striated slip surfaces within the damage zone indicates that all of them belong to the same deformation event. Such synthetic, antithetic faults and shear fractures form a network where planes provide a preferred structure for fluid circulation. Although generally the preferred

orientation for fluid circulation in fracture networks is parallel to the faults and fractures intersection (Sibson, 1996), in this area, because of its volcanic origin, a strong degassing process pushes CO2 rich fluids to the surface from bottom to top (i.e. as evinced by the deposits of sulfur crystals along the fault plane; Fig. 3b). The so called Mature Fault Zones are characterized by spacing values of the order of tens of meters (see Fig. 3c) and well developed core zones. The Mature Fault Zone cropping out in the quarry consists of two main slip planes which separate a 5 m wide fault core from the two surrounding damage zones that are 5–6 m wide each. These planar–slip surfaces are approximately parallel and sub-vertical, dip toward 150◦ , and have abrasive slickensides that indicate a normal sense of shear. The core zone, which consists of mainly gray, clayey material mixed with small lithic fragments, can be considered a low-permeability zone as evidenced by its low CO2 flux rate compared to the surrounding damage zones (Fig. 3c). The greater permeability in the Younger Fault Zone, as suggested by the higher CO2 flux rates measured across this interval, strongly advocated the use of data from one of these damage zones to build

186

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

Fig. 2. Geology of the Latera caldera field. (a) Geological map; (b) regional soil gas survey in the caldera; (c) faults and fractures measured in the area. Note that the main discontinuities have patterns similar to the gas distribution in the soil (after Annunziatellis et al., 2008, modified).

a 3D depiction of a fractured outcrop and simulate fluid flow across a fractured rock volume. 3. Methodology Data acquisition techniques and methodologies for studies of fractures can be distinguished by the scale of observation and dimension of fractures (van Dijk, 1998; van Dijk et al., 2000). In many cases, experimental data are collected in one or two dimensions from outcrops (cm to m scale) or boreholes to obtain information on fracture density, spacing, orientation, aperture,

nature of filling, and alteration in order to characterize a fractured rock volumes. A standard, one dimensional technique is the scan line survey (Priest and Hudson, 1981), where scans along vertical or horizontal lines are employed to describe the cumulative distribution of fracture spacing, fracture length and fracture aperture (La Pointe and Hudson, 1985; Marrett et al., 1999; Ortega et al., 2006; among many others). In outcrop-based studies, quantitative parameters for fracture analysis (e.g. fracture spacing/fracture length, fracture aperture) can also be acquired by two dimensional analysis, considering distribution with respect to an area instead of a single line. Area sampling involves mapping the fracture trace

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

187

Fig. 3. Feldspar and sulfurs quarry in the Latera Caldera with the location of two different types of fault architectures recognized in the area. (a) Young fault zone, a high permeability fault zone; note the yellow patch, evidenced with the squared symbol, corresponding to sulfur deposits deposited along the fracture planes; (b) mature fault zone showing a clay rich core as product of mechanical and geochemical reworking of fault rocks, suggesting very low permeability values; CO2 flux measurements across the fault zone show that flux values are very low in the zone corresponding to the core of the fault (after Annunziatellis et al., 2008) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article).

patterns and recording fracture characteristics at locations in the area mapped (Wu and Pollard, 1995). The 3D reconstruction of a fractured volume, instead, is usually carried out using computer software, able to generate fractures networks in synthetic rock volumes (Caine and Forster, 1999). In this work, we create a 3D model of the outcrop described in Section 2.1 (Fig. 4). The construction of this model, called an Analogue Model (AM), is carried out in two steps: (1) field work analysis and (2) data processing. The name Analogue Model indicates that the approach to the three-dimensional construction of the fracture network adopted here is deterministic and based on values measured directly in the field; therefore it can also be described as a Data Driven Fracture Model, as already proposed by van Dijk (2002). 3.1. Fieldwork analysis We performed a 2D mapping of the chosen fracture network on two adjacent, orthogonal surfaces of the outcrop that are 3.5 m × 3 m and 2.5 m × 3 m in size. On each surface we used c. 0.5 m × 0.5 m sampling windows to identify and map fractures planes, first on rectified photographs of the rock surface, then numbered, classified and drawn using one or more lines (Fig. 5 and Fig. 1 of Supplementary Data). For each mapped fracture we measured (a) fracture position and trace (length), orientation, aperture and fill (distinguished in classes); (b) opening mode; (c) roughness and curvature (Table 1). In addition, six scan lines, two for each surface

(running parallel to the directions of the main outcrop walls) and two on the quarry walls, were also measured to have data for the stochastic elaboration and DFN reconstruction using Petrel 2009 and Move 2011 (Figs. 3 and 5 for location). Along these scan lines we collected data on the main orientation, length, and aperture of Table 1 Data measured in the field at wall 1 and to wall 2, separated in four fractures groups (attitude is expressed as dip direction – dip). Closed

Open

Shear fracture

Small fault

334–83 295–85 120–60 140–89 190–65 154–65 340–88 315–85 107–85 190–15 320–75 203–85 190–65

350–55 220–70 222–88 328–89 331–89 130–50 329–80 100–80 198–80 140–80 200–85 150–89 150–89 150–89 350–82 100–80 219–85 225–88

145–80 155–85 110–88 210–88 100–60 375–70 150–88 315–88 300–85 350–78 330–80 350–55 300–55 310–83

310–75 40–80 40–80 312–89 312–89

188

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

negligible confining pressure; for simplicity, we assign a single aperture value to each open fracture (Figs. 6 and 7). 3. shear fractures include mode II and III discontinuities, filled by white gouge material. The gouge material homogeneous in grain size and is constituted by clays created by the mechanical and chemical alteration of the lavas cropping out in this area. The aperture values range from a few mm to a few cm, while fracture planes strike around N50◦ and have dips near 70◦ (Figs. 6 and 7); 4. small faults include mode II and III faults. Kinematic indicators on the fault planes suggest a normal sense of motion with a small horizontal slip component (pitch angles around 60◦ ). Unfortunately, the lack of markers did not allow us to evaluate the range of offset, except for a few cases where it was found to be on the order of a few centimeters. Faulting zones range in width from a few cm to a few tens of cm, are characterized by arrays of anastomosed slip planes developed inside the fault zone, and are filled by white gouge material with grain size ranging from micro breccia to clay. Small faults have trends and dips similar to those of shear fractures (Figs. 6 and 7). 3.2. Development of the Analogue Model

Fig. 4. Workflow followed in this study. The field work is at the base of the reconstruction with AutoCAD of the three-dimensional fractured volume, and the DFNs built using Petrel 2009 and Move 2011. Only the Analogue Model was used for fluid flow simulation (see text for detailed explanation).

fractures, the morphology and frequency of fault plane discontinuities, and the number of fractures per unit length and their average spacing (Fig. 5 and Table 1 of Supplementary Data). Considering that the effect of surface orientation on the distribution of fractures is especially important in the case of sub vertical outcrop surfaces, we chose two outcrop walls perpendicular to each other and performed (as described above) a number of additional scan lines measured in the same quarry and oriented parallel and normal to the main tectonic direction in the area (Figs. 3 and 5; Supplementary Data, Table 1 and Table 2; see van Dijk, 2002). This approach allows us to reduce bias due to surface orientation, since we directly measured the fracture attitudes. We identified four different fracture groups, based on morphological, kinematic and genetic features: open, closed, shear fractures and small faults (Fig. 6). This classification was performed in an effort to correlate the observed structural features and the permeability of the measured fractures and fault planes. Fracture aperture and roughness are in fact highly elusive parameters, so we tried to constrain them using field data. We assigned four different aperture values to the four classes recognized. When possible, minimum and maximum fracture aperture values were measured in the field: 1. closed fractures include discontinuities belonging to the opening mode I, with apertures up to 0.5 mm, orientations around N150◦ , and dip angles around 80◦ (Figs. 6 and 7); 2. open fractures include mode I discontinuities with no filling, apertures higher than 0.5 mm, and orientations mainly around N50◦ and N140◦ . Some fractures in this class have quite large apertures, up to a few centimeters, probably because of

The 2D fracture lines were drawn from rectified photographs as 3D planes in AutoCAD, using the orientation data of each fracture measured in the field. The lines identifying each fracture were picked and correlated one by one manually. Although timeconsuming, the manual approach made it possible to image the fractures as warped surfaces. Fracture planes crossed or cut each other within the AutoCAD three-dimensional space, ending at external faces of the considered volume. Although this may not be exactly what we observed in the field, this choice allowed us to keep the fracture density closer to that obtained from the scan lines analysis, as can be deduced by comparing the values of P32 and P10 (respectively, P32 = 4.32; P10 (scanlines 1, 2, 3, 4) = 5.1; 6.9; 3.8; 5.6) according to the relationship that spacing is equal to the inverse of P32 for sub parallel fractures (Fig. 8a and Fig. 1 of Supplementary Data). Generally, fractures in a 3D reconstruction are disk shaped, especially if the model will be used to develop flow or transport models (Baecher, 1983; Zhang et al., 2002). We decided to assign each plane the shape resulting from the extension in the third dimension of the intersection between the fracture plane and the outcrop surface (Fig. 8a). This was done because some of the measured fractures were longer than the facets of the studied volume. Adapting them to a disk would have significantly reduced the surface extent of the fractures. The Analogue Model shown in Fig. 8 is the digital 3D representation of the rock volume analyzed in the field. Note that we use only the fracture data measured and mapped in the field, without any correction or adjustment. The resultant Analogue Model can now be evaluated against fracture models obtained with more common procedures in terms of number of fractures (discrete features) and permeability distribution (model upscaling). We focus on two procedures: (a) Discrete Fracture Network models (DFN) generated by Petrel 2009 and Move 2011 – see Section 4; (b) a representative elementary volume (REV) depiction of the fractured medium, see Bear (1972), created by upscaling the Analogue Model using Petrel 2009 – see Section 5. 4. Discrete Fracture Network (DFN) models We built several DFN models using statistical tools implemented in Petrel 2009 (a software package by Schlumberger employed to study oil/gas reservoirs) and Move 2011 (a software package by Midland Valley employed in structural geology), and compared

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

189

Fig. 5. Mapping of the two sides of the outcrop. (a) The two walls with the orientation of the outcrop and the location of the four scanlines; (b) mapping of faults and fractures, numbered in the field on the two walls (each grid cell is a square with a side of 0.5 m; the shape is irregular due to the photo distortion); (c) all measured and mapped data, symbols indicate the four different fracture groups identified in the field; (d) data for faults and fractures along the four scanlines; the orientation of the two outcrops is signed in all the projections (Schmidt projection, lower hemisphere). All data are available in Tables 1 and 2, and Supplementary Data.

and validated the DFNs’ number of fractures, fracture size distribution, and fracture density against the Analogue Model. Petrel 2009 and Move 2011 require the following initial input parameters: the fracture intensity, mean orientation for each set of fractures (obtained by the fracture sets orientation analysis), and the statistical fracture-size distribution. The required input data was taken from the scan lines and the two scan areas (the two walls of the same outcrop) measured during our field work in Latera (see Supplementary Data, Tables 1 and 2). The fracture intensity is defined as the amount of fracturing per unit volume, and is given by the number of fractures per volume (parameter P30 in Petrel 2009), the fracture length per volume (P31), or the fracture area per volume (P32). The reference volume used for all the DFN models is that measured in the field (∼26 m3 ; Table 2). Petrel 2009 generates DFN models using one

of these intensity parameters. In the present work, we first calculate the values of P30, P31, and P32 for each of the four different fracture sets recognized in the field (Open, Closed, Shear Fractures and Small Faults; see Section 3.1). These parameters are employed to build the DFNs shown in the Fracture Set Model of Fig. 9. Then we use the mean density values of the parameters to build a Fracture Network Model where all fractures are considered belonging to a single set with the same aperture and permeability (Fig. 9). This is unlike borehole intensity fracture analysis, where fracture count can be biased by the relative orientations of well and fractures. The recognized fracture groups are based on aperture values and not on different directions; in general for equally spaced families, fractures oriented at a high angle (e.g. perpendicular) to the trace will be measured more frequently along the scan line, whereas low angle fractures (sub-parallel) will be measured much

190

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

P32 – were first created independently for each of the four different fracture sets (Open, Closed, Shear Fractures and Small Faults) and then combined together in a single DFN model. The other three models were based again on either P30, P31 or P32, but the DFN was created merging all fracture sets in a single group. Four DFN were created using Move 2011, based either on P30 or P32, following the same procedure described in the paragraph above (Fig. 9). Comparison with the Analogue Model indicates that the number, area, attitude and length of the fractures are different for each DFN generated by Petrel 2009 and Move 2011, as a function of the fracture density values (P30, P31 and P32). In the case of P32, the number of fractures in the DFN models is equivalent to that of the Analogue Model. 5. Upscaling Upscaling converts a DFN (and its previously defined properties, such as facture aperture or length) in new properties required for simulation of a dual porosity or dual permeability model within a fluid flow simulator (e.g. Eclipse by Schlumberger, or any other equivalent). New properties, such as porosity, permeability, and sigma coefficient (representing connectivity) are assigned to a reference grid. The optimal dimension of this grid is also evaluated during the upscaling. Our Analogue Model can be considered as a simple fracture set created deterministically as the result of fault planes extraction from seismic images. When deterministic methods are used, as in this work, Petrel 2009 performs an approximation of the input fracture plane when it deviates from planar, simplifying the shape by calculating a planar best fit surface based on the input data. In this way all the input fracture planes were converted to planar objects. This process resulted in a simplification of the Analogue Model, although the spatial relationships were preserved. 5.1. Analogue Model upscaling

Fig. 6. Examples of each group of fractures recognized in the outcrop of Fig. 5. (a) Closed fracture; (b) open fracture; (c) shear fracture; (d) small fault.

more rarely. To minimize these aspects, we performed scan lines along two different orthogonal directions. The data (localization, length, orientation, and classification) acquired for the scan lines are summarized in Table 2 and in Tables 1 and 2 in the Supplementary Data. Fracture length distributions and fracture orientations are the other two important input parameters for the generation of the DNFs. The Fisher model provided average dip, dip azimuth values and concentration for fractures orientation (Table 2). Fracture lengths measured in the field follow a log-normal cumulative size distribution. The log-normal distribution is implemented in Petrel 2009 but not in Move 2011. Only the normal distribution is implements in Move 2011. Six DFNs were built using Petrel 2009. The first three DFN models, each based on one single parameter – either P30, P31 or

Our Analogue Model was imported into Petrel 2009 as a deterministic fracture model. The spatial coordinates defining each fracture in the CAD model were exported as a text file and imported into Petrel 2009 (Fig. 8b and c), where they were used as a DFN. Geological parameters (i.e. field measured aperture, length, dip and dip azimuth) were assigned to each fracture of the imported network. Finally, permeability values were allocated to each group of fractures. The Petrel 2009 code applies information on where and how the fractures behave in 3D space, using high confidence fracture patches from the initial Analogue Model. To maintain the spatial relationships between properties in adjacent cells, we modeled the fractures explicitly (deterministic method). When the deterministic method is used, Petrel 2009 performs an approximation of the input fracture plane when it deviates from being planar. Therefore, if complex fault surfaces are used as input they can be slightly altered to better fit a planar surface (Caine and Forster, 1999). As described above, the software simplifies the shape by calculating a planar best fit surface based on the input data. The physical and geometrical characteristics of each fracture and their spatial distribution are the main parameters needed to describe the permeability distribution within a volume of fractured rock. Although roughness is considered one of the main features that control permeability in many laboratory studies (e.g. Maerz et al., 1990; Lanaro, 2000; Lapcevic et al., 1999), its effect decreases as the fracture aperture increases and can therefore be ignored in many practical applications (Bandis et al., 1985). Aperture can be reasonably assumed to be the main controlling parameter for permeability, neglecting roughness for all the fractures having

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

191

Fig. 7. Some examples of the brittle deformation in the area of the caldera of Latera. (a) Small normal fault dissecting ash deposits; (b) lava deposits deformed by extensional shear fractures system; note the isolated patch of sulfur on the outcrop wall, suggesting a main horizontal flow direction; (c) small, closed normal fault offsetting a tuff bed; (d) open fractures in lava deposits; (e) high angle strike slip fault.

an aperture different from zero, especially if we assume that the porosity of the host rock is negligible (Long et al., 1996; Aydin, 2000). We assigned four different values of aperture and permeability to each fracture group recognized in the field (Table 2). For the closed and open groups of fractures we used an aperture corresponding to the maximum value measured in the field, which is 0.5 mm and 10 mm respectively. For the other two groups, the fault rock material filling the Shear Fractures and Small Faults complicate the aperture – permeability relationship. For the Shear Fracture group, the average aperture values measured in the field are similar to those of the Open group, but they are systematically filled by gouge material. To simulate the effect of the low permeability fill, we assigned the intermediate aperture of 5 mm, smaller than the values measured in the field. The same criterion was applied for the Small Fault group, to whom we assigned an aperture of 15 mm. These faults are characterized by high aperture values in the field and, for this reason, were considered the most permeable structures in the reconstructed volume.

If we approximate a fracture by parallel plates, the laminar flow in the fracture can be described by the following cubic equation (Snow, 1969; Witherspoon et al., 1980): Q =

b3 P P = b 12  

(1)

where Q is the flow within the fracture, b is the aperture (distance between the two walls of the fracture), P is the pressure gradient, and  is the fluid dynamic viscosity. The fracture permeability  is approximated by =

b2 12

(2)

Permeability values, expressed in mD, were calculated in Petrel 2009 for the four groups of fractures using (2) (see Table 2), and assigned to each fracture of the Analogue Model imported from AutoCAD. We employed the flow based analysis, implemented in Petrel 2009, for the upscaling process (Dershowitz and Fidelibus, 1999).

192

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

volume equal to or larger than 1.5 m × 1.5 m × 2.5 m. The outcrop volume sampled in the field (3.5 m × 2.5 m × 3 m) should thus be sufficient to estimate the permeability distribution (Fig. 10). 5.2. DFNs upscaling During the upscaling process we converted the Discrete Fracture Network models calculated by Petrel 2009 and Move 2011 to equivalent porous medium properties over a grid. The porous medium volume, considered equivalent to that measured in the field (∼26 m3 ), is formed by 208 cubic cells with a side of 0.5 m (Fig. 11). The fracture permeability in an upscaled cell is the effective fracture permeability based on all the fractures in that cell. The result of the upscaling process is a grid with equivalent properties for porosity and permeability. The assignment of the hydraulic properties (we use as input parameters P30, P31 and P32 values, although for simplicity we discuss only results for P32) to the grid nodes was carried out in Move 2011 using the Oda method (Oda, 1985) and in Petrel 2009 using the Flow Based Method (Dershowitz and Fidelibus, 1999). The permeability tensor is calculated by assuming that the fracture network is well connected; note that it does not take the fracture size and the connectivity into account. This means that a potential permeability is generated for each grid cell assuming that all fractures are connected and that permeability is independent of fracture size. Oda’s fracture tensor represents the volume weighted average product of fracture areas and transmissivity. The more accurate Flow Based Method is computationally challenging compared to the Oda method, being designed for a maximum of a few thousands of cells. This flow simulation is actually practical in our model given its relatively small number of cells and fractures. 5.3. Comparison of the upscaling process

Fig. 8. Analogue Model. (a) Model obtained by the generation of each plane mapped in the field (the grid corresponds to 0.5 m) in AutoCAD; note that plane intersected each other, are warp shaped and cut at the end of the model volume; (b) the same model with the four groups of fractures evidenced by different colors; (c) Analogue Model imported in Petrel 2009. Here the colors correspond to different aperture/permeability values (defined in Petrel 2009).

This method creates a finite element grid for each grid cell and simulates flow under a pressure gradient to calculate the permeability in all directions. This simulation, done separately in each of the 6 directions for a full tensor calculation, allows us to define the representative grid dimension for our Analogue Model and permeability distribution. The grid-block permeability may be strongly influenced by the boundary conditions imposed on the flow equations and the size of the grid-blocks. We used four different cell sizes (see Fig. 10a–d) and four different boundary conditions (see Fig. 10e–h). The results indicate that the Analogue Fracture Network is representative of a permeability distribution within a rock

We applied either the Oda or the Flow Based method to all the obtained fracture network models: (i) the Analogue Model; (ii) the Fracture Set Model, the DFN created by combining the stochastic models for Open, Closed, Shear Fractures and Small Faults groups, and (iii) the Fracture Network model, the DFN obtained merging and modeling all fractures together without any difference. The results of the upscaling process for the analogue and stochastic models are shown in Figs. 11 and 12. Permeability values were calculated independently for each cell. These calculations can contain bias and high approximations, especially when a cell has a small number of, or sparse, data (e.g. the Fracture Network Model in Fig. 12a). In general, the upscaled permeability distribution for the Fracture Set Model is closer to the permeability distribution of the Analogue Model. Upscaled permeability values from all the stochastic DFNs are always smaller than those from the Analogue Model (Fig. 12). In this case, upscaled permeability values estimated from Fracture Network models do not compare well with the field data from the Analogue Network model: permeability estimates are two orders of magnitude smaller and do not show the strong anisotropy of the Analogue Model. The Fracture Set Model shows a more homogeneous distribution of permeability (see Fig. 12b). In the case of the Flow Based method, the permeability distribution of the two DFN (Fracture Set and Fracture Network models) reproduce the strong anisotropy seen from the field data of the Analogue Model (Fig. 12c). Permeability values along the horizontal direction i are higher than those along the horizontal j and vertical k directions, indicating a preferred path for gas migration that is parallel to the direction of intersection of the main conjugate faults in the volume (e.g. Sibson, 1996).

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

193

Table 2 Parameters calculated for the four Fractures Sets (one for each fractures group) and for the Fracture Network. Mean Vector is the average values of dip and dip direction obtained with the Fisher method implemented in Move 2011 (Midland Valley). Concentration values obtained by contour analysis performed with Georient 9.5.0 for each Fracture Set (Open, Closed, Shear Fractures and Small Faults) and the Fracture Network (all fractures together).

Aperture (mm) Permeability (mD) P 30 P 31 P 32 Mean vector (◦ ) Density (%)

Closed

Open

Shear fractures

Small faults

Fracture network

0.5 2.08 × 108 0.433 0.320 0.595 28–213 31

10 8.33 × 109 0.633 0.434 1.011 58–185 39

5 2.08 × 109 0.533 0.505 1.827 38–317 36

15 1.87 × 1010 0.166 0.230 0.9993 80–345 60

– 5.33 × 109 1.766 1.489 4.432 22–285 30

6. Numerical models of CO2 leakage We developed a numerical model of CO2 flow through the fractured outcrop using COMSOL Multiphysics (COMSOL AB), a Finite Element Analysis software package that simulates flow in fractured reservoirs using a dual porosity approximation. COMSOL represents the fractures as boundaries between adjacent matrix blocks. Darcy’s law governs velocities in the matrix blocks, while the flow along the blocks boundaries is determined taking into account fracture aperture. Representing the fracture as an interior boundary is particularly efficient, because it eliminates the need to create a geometry with a high aspect ratio. A very long and narrow fracture domain would otherwise require a dense mesh consisting of great number of tiny elements (COMSOL, 2011). The time-dependent fluid flow in the porous block matrix is controlled by Darcy’s equation [f  + s (1 − )]

∂p −∇ · ∂t



m





∇p ;

ub =

m ∇p 

(3)

where P is the fluid pressure gradient in the pore space,  is the porosity of the matrix blocks, f and s are the compressibility of the fluid and solid respectively,  gives the permeability of the

matrix blocks and  is the fluid’s dynamic viscosity. The fluid velocity ub gives the Darcy’s velocity, which is a volume flow rate per unit area of porous media. The fractures in the model are a sequence of interior boundaries. For this type of analysis to be valid, the equation for the fracture must solve for the same dependent variable p as the equation for the matrix block. The equation for velocity in the fracture follows a modified form of Darcy’s law, where the coefficients take into account the relatively small resistance to flow along the fractures and the fracture apertures, which gives dimensional consistency between the fractures and the porous matrix: Sfr

∂p −∇ · ∂t



fr





dfr ∇ p ;

uf = −

fr 

∇p

(4)

where Sfr is the fracture-storage coefficient, fr describes fracture permeability and dfr is the fracture aperture (Table 3). Initially, numerical tests of the flow properties of the fractures were performed on a homogeneous porous block whose permeability fr in the flow direction was set to 5 × 10−3 mD (Table 3). A homogeneous model, while unrealistic, enables the model predictions to be checked against the expected theoretical results. A pressure gradient, equal to the lithostatic pressure difference between the

Fig. 9. Discrete fracture networks built by Move 2011 and Petrel 2009, using as input parameters the field data for parameters P30, P31 and P32 (see text for parameters description). Five models were generated by combining all fractures together (Fracture Network models). Five more by building a DFN for each fracture group and then combining the 4 DFNs in a single one (Fracture Set Models).

194

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

Fig. 10. Upscaling analysis performed in Petrel 2009. Grid dimension is based on the permeability distribution calculation. Used cell size: (a) 0.5 m × 0.5 m × 2.5 m; (b) 1 m × 1 m × 2.5 m; (c) 1.5 m × 1.5 m × 2.5 m; (d) 3.35 m × 0.85 m × 2.5 m. Permeability values obtained for the different cell sizes: (e) linear pressure drop; (f) constant pressure drop; (g) periodic pressure drop; (h) no flow on lateral boundaries. Note that the permeability values are coincident for a cell size of 1.5 m, the representative grid dimension.

bottom and the top of the block, was applied in the vertical direction (z-direction). This pressure gradient is the minimum requested to obtain a flow through the model itself. It is not representative of the real pressure gradient of the gas in the field, which is likely to Table 3 Parameters used in the COMSOL Models. (1) Zinszner and Pellerin (2007); (2) Bear (1972); (3) laboratory data; (4) Brunton (1988); (5) measured values. Permeability values are expressed in millidarcys (mD). Matrix porosity (1) Matrix permeability (2) Matrix density (3) Matrix compressibility (4) Aperture group “Closed” Aperture group “Open” Aperture group “Shear Fracture” Aperture group “small Faults” Permeability group “Closed” Permeability group “Open” Permeability group “Shear Fracture” Permeability group “small Faults”

0.45 5 × 10−3 mD 1.133 g/cm3 10−13 Pa−1 0.5 mm 10 mm 5 mm 15 mm 833,333,333 mD 8,333,333,333 mD 2,083,333,333 mD 18,750,000,000 mD

be higher. The CO2 moved through the block with a homogeneous, constant velocity. It is important to keep in mind that it is always possible to employ more complex equations describing flow behavior. For example, we could consider the permeability as a function of time, or assume that CO2 is compressible. Future development of our research could provide the opportunity to employ a more realistic mathematical description when modeling the fluid flow (e.g. Putra et al., 1999; Konzuk and Kueper, 2004). The architecture of COMSOL Multiphysics is quite flexible, allowing one to build more complex models in a step by step fashion that employs progressively more complex physics while maintaining the initial integration volume. We imported each fracture of the Analogue Fracture model from the AutoCAD file into COMSOL Multiphysics, maintaining the geometry developed from the field data (Fig. 4). We carried out five simulations to measure the different flow properties of these heterogeneous fracture distributions (Fig. 13). The first simulation shows the flow of CO2 in a fractured medium composed exclusively by closed fractures (Fig. 13c). These are the fractures with

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

195

 (a) Analogue model; (b) fracture set (P32) generated in Petrel 2009. Note that the higher Fig. 11. Permeability distribution in the horizontal (i, j) and vertical direction k. permeability zone (in red) seems to be similar in the two models, but it is more evidenced in the Analogue Model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

the smaller permeability (∼10−8 m2 = 10,132.5 D; 1 D is equivalent to 9.869233 × 10−13 m2 ) and a aperture of 0.5 mm (see Section 3.1 on petrophysical characterization for additional details). These fractures slightly change the pattern and rates of the CO2 release

at the top of the block (Fig. 12a). The second simulation shows the flow of CO2 in a fractured medium where all the fractures belongs to the Shear Fractures group (thickness of 5 mm and permeability around 10−6 m2 = 1,013,250 D). The CO2 leakage is now

196

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

7. Discussion

Fig. 12. Different permeability distributions the Analogue and Stochastic models in  (a) Results obtained using Oda method the horizontal (i, j) and vertical direction k. applied to the Analogue Model and to the DFN Fracture network (all fractures as one set and the P32 density) in Move 2011 and Petrel 2009; (b) results obtained using Oda Method applied to the Analogue model and to the DFN Fracture Set (four sets modeled separately and than joint together and the P32 density) in Move 2011 and Petrel 2009; (c) results obtained using No flow Method applied to the Analogue Model and to the DFNs Fracture Set (four sets modeled separately and than joint together) and Fracture Network (all fracture together) and P32 density in Petrel 2009. See text for discussion.

influenced by the fractures network, with the largest flow values at the outlets of the fractures. In the third simulation, we add the open fractures to the block. These fractures have a thickness of 10 mm and a permeability of 10−5 m2 (=10,132,500 D). The fractures are now the preferential routes for the flow of CO2 , with the largest gas leakage values occurring at the intersection of two or more fractures (Fig. 13b). The fourth simulation measures the influence of the fractures classified as Small Faults, with a thickness of 15 mm and a permeability on the order of 10−5 m2 . Once again, the CO2 flow is controlled almost exclusively by the fracture network. The fifth and final model includes all the fractures that have been identified, mapped and classified in the field. As we observe in the results (Fig. 13e), the CO2 release distribution is more heterogeneous that in all the other models. The gas flow depends on the pressure gradient that it experiences and the changes in permeability. Hence, the velocity field will change spatially and with time as the fluid moves from one end of the block to the other.

One of the more complex aspects for the assessment of carbon capture and storage (CCS) sites is the evaluation of the reservoir permeability, controlled by the petrophysical parameters of the fractured rock, and the fractures/faults location, length, depth, state of stress. Research is usually focused on the reconstruction of the distribution of these parameters, based on high resolution seismic lines analysis and interpretation, integrated with field work studies (Cavanagh and Ringrose, 2011; Ringrose et al., 2011). These methodologies can reconstruct the discrete network and its parameters at regional scale. Equally challenging is the estimation of the permeability for a fractured volume of rock, especially in view of evaluating leakage potential. Fracture modeling modules use fracture data to describe large volumes of fractured rocks, in order to evaluate capability volume. Permeability is treated as related to aperture, the best known parameter. Other aspects, such as connectivity and geometry of the fracture network, are generally are generally poorly constrained. In our work, we tried to reproduce the natural fractured volume, being careful to preserve the geometry and the physical parameters of each fracture observed in the field. Statistical parameters were only used to evaluate the reliability of our approach and to compare our numerical modeling of CO2 leakage to soil gas surveys results. The simulated leakage of CO2 shows a number of similarities with the gas leakage patterns measured in the field. Gas leakage at the surface is concentrated in points aligned along a main discontinuity. These points at the surface correspond to gas vents and are usually aligned along the trace of the major faults and/or along highly fractured fault zones. These spotty anomalies (isolated points with high concentration values) are frequently observed and documented in the literature (Ciotoli et al., 2007; Zhou et al., 2010). Several gas vents are identified in the surveys of Figs. 2 and 14a, showing leakage at the regional (the area of the entire caldera of Latera) and local scale (the surface of the studied quarry). The map of gas emissions from the numerical experiments represents the smaller scale of the outcrop but the qualitative distribution of the CO2 leakage is very similar to that observed in the field, where leakage along the fractures is characterized by an irregular distribution of points of higher and lower emission (Fig. 14b). This irregular distribution can be observed in all the results of our numerical experiments (Fig. 13). Spotty anomalies are aligned along main discontinuities and at the intersection (or dip change) between planes belonging to the same fracture group. Higher leakage points are controlled by the intersections of high permeability or steeply dipping fractures planes (Figs. 13 and 14b), which means that not only the permeability but also the connectivity of a fracture system plays a key role on fluid flow and leakage control. In our Analogue Fracture Model (see Fig. 11, Ki direction), the upscaling analysis indicates the maximum permeability in the Ki direction, parallel to the intersection of the main conjugate fracture system and normal to the direction of the imposed flow (from bottom to the top). The flow is deviated repeatedly from vertical to horizontal and vice versa, as the gas finds a new discontinuity as preferential path to the surface (red arrows in the model of Fig. 13). Moreover the different sizes of the red arrows suggest a variation of flux intensity along the gas path, which, in nature, can imply pressure drops and mineral precipitation. This behavior can be compared to that observed at the outcrop scale at the quarry site at Latera. Here is not rare to observe sulfur deposits along the fracture surfaces, along a vertical path. The same traces, which indicate the occurrence of a pressure drop and deposition along fractures, suddenly disappear, suggesting a deviation of the main gas flow in a direction normal to the outcrop wall (Fig. 7c).

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

197

Fig. 13. Numerical modeling of gas flow in fractured media. Arrows indicate the flow direction; the arrow size indicates the flow velocity (m/s). Plots show the gas leakage at the free surface (upper boundary of integration volume): (a) open; (b) shear fractures; (c) closed; (d) small faults; (e) all fractures together, (f) fracture network imported from Petrel 2009. (For interpretation of the references to color in text, the reader is referred to the web version of the article.)

198

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

Fig. 14. Comparison between the soil gas survey and numerical model (e) from Fig. 12. (a) Gas release measured in the quarry (after Annunziatellis et al., 2008, modified). (b) Gas release distribution at the free surface (upper boundary of the integration volume) of the Analogue Model. The area covered by (b) corresponds to the black square in (a) indicating the location of the two outcrop walls where the field work analysis has been performed. Note that it comprises the main fault (white line and white ellipses in a) and its damage zone, corresponding to an area of leakage. The spotty anomalies characterizing gas emissions along this fault in (b) reproduce at smaller scale the trend of emissions along the main faults in the quarry, evidenced by white lines and white ellipses.

8. Conclusions In this work, we developed 3D Analogue Models of fractured rock volume by integrating field work and 3D fractures digitization. The approach to the three dimensional construction of the fracture network is deterministic (Data Driven Fracture Model) (van Dijk, 2002), and is based on values measured directly in the field. The applied methodology requires a minimal probabilistic component as well, because the reconstruction of the fracture network in a CAD, and the upscaling of the permeability, requires assumptions and simplification. The approach proposed is limited by the availability of a good outcrop. Time in the field is also a critical factor since the scale of observation range from centimeters to meters. Moreover, the fracture network reconstruction using a CAD (in our case AutoCAD 2010 by Autodesk) is time consuming. The approach adopted in measuring and preparing the field data can significantly influence the result. Field work is the critical component and could be improved in several ways: for example, employing measuring techniques that can help in preserving the geometrical relationship observed in the outcrop. Field analysis has the potential to classify fractures based on their relative permeability. This information can later be used to validate and verify fracture data collected and interpreted at larger

scale (∼km). Compared to other methods employed to image fractures at a similar scale, such as CT, X-ray or nuclear magnetic resonance (NMR), data are here completely acquired and verified manually. This allows us to capture a level of detail that cannot be reached by other methods, especially regarding the qualitative classification of fractures. In fact, the comparison between the Analogue and the stochastic Discrete Fracture Network models (Figs. 11 and 12) shows that the permeability distribution from stochastic DFN is closer to that of the Analogue Model when the different fracture sets can be treated and modeled separately. It should also be emphasized that although we stress the major role played by field work and manual data analysis in our work, we also made use of different software and modules to create the DFNs and run the numerical experiments. Some of these are specifically designed to evaluate the petrophysical feature of a reservoir, based mainly of the statistical approach. Moreover some of them are limited and we used them only to develop some steps of the entire workflow, to avoid unrealistic results. From a geometrical point of view, the Analogue Model shows the occurrence of anisotropies better than the statistical models. This is expected even in the upscaling results: here, in fact, DFN models from Move and Petrel estimate permeability values of about two orders of magnitude lower than the Analogue Model, and the permeability distribution does not reproduce the anisotropies seen in the field, especially when the aperture values are not considered (Fig. 12). Taking into account this latter aspect, it is clear that the statistical approach strongly depends on aperture values and their distribution within a fracture set, and give less importance to fracture geometrical relationships. The four different permeability values inferred from field data provide a quantitative indication of the control exerted on the gas migration by each group of fractures, in terms of both fluid volume and fluid pattern. Numerical experiments show that the fracture dip variation, as well as the orientation with respect to the flow direction, can also control the gas flow. In the models developed in COMSOL, fluid flow initially follows the direction of maximum permeability (horizontal in the model, corresponding to the direction of intersection of the conjugate fracture system), but the direction of flow can easily be deviated at the intersection with more permeable fractures (Fig. 13e); in some cases, as suggested by the different dimensions of the red arrows that represent the gas path, it also influences flux intensity and can be associated with fluid pressure variations. The result (see Fig. 13e) is a spotty distribution of fluid emission at the surface, a feature recognized in soil gas field surveys (Fig. 14). The fractured medium model used for the numerical experiments derived from the Analogue Model has a complex and irregular geometry, based exclusively on field data. The numerical simulations shown in Fig. 13 allow us to identify three major parameters controlling the migration pathway of the fluid: (1) the permeability of each discontinuity (from the field analysis); (2) the attitude of the planes (much more complex than in commercial software); and (3) the relative position of fractures within the model volume itself. The possibility to import the geometries measured in the field (even though still approximate) has highlighted the role of geometries as a controlling factor of connectivity, and as a consequence, of permeability distribution. To conclude, the Analogue Model may be time consuming to build but preserves distinctive characteristics of the natural fracture networks that are usually smoothed by the statistical analysis and the upscaling process. The possibility to preserve the geometry observed in the field is critical in the case of studies dedicated to the storage of CO2 . The identification of possible leakage structures, as well as of potential reactivated structures, is a critical constraint during both injection and the long term monitoring of the reservoir (Ringrose et al., 2011).

S. Bigi et al. / International Journal of Greenhouse Gas Control 18 (2013) 183–199

Acknowledgments Free academic licenses for Move 2011 and Petrel 2009 were generously made available by Midland Valley and Schlumberger; comments by two anonymous reviewers greatly helped to improve our manuscript. We thanks to all the people working at the Fluid geochemistry Lab of Earth Science Dept. of Sapienza University of Rome, and particularly to Stan Beaubien, for the their support, the useful discussions in the field and the revision of the manuscript. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:10.1016/j.ijggc.2013.07.011. References Annunziatellis, A., Beaubien, S.E., Bigi, S., Ciotoli, G., Coltella, M., Lombardi, S., 2008. Gas migration along fault systems and through the vadose zone in the Latera caldera (Central Italy): implications for CO2 geological storage. International Journal of Greenhouse Gas Control, http://dx.doi.org/10.1016/ j.ijggc.2008.02.003. Aydin, A., 2000. Fractures, faults, and hydrocarbon entrapment, migration, and flow. Marine and Petroleum Geology 17, 797–814. Baecher, G.B., 1983. Statistical analysis of rock mass fracturing. Journal of the International Association for Mathematical Geology 15 (2), 329–348. Bandis, S.C., Barton, N.R., Christianson, M., 1985. Application of a new numerical model of joint behavior to rock mechanics problems. In: Proc. International Symposium on Fundamentals of Rock Joints, Bjorkliden, pp. 345–355. Barton, C.A., Zoback, M.D., 1992. Self-similar distribution and properties of macroscopic fractures at depth in crystalline rock in the Cajon Pass Scientific Drill Hole. Journal of Geophysical Research, B 97 (4), 5181–5200. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier Publishing Company Inc., New York, USA. Brunton, G.D., 1988. Density and compressibility of Wyoming Bentonite particles. Clays and Clay Minerals 36 (1), 94–95. Caine, J.S., Evans, J.P., Forster, C.B., 1996. Fault zone architecture and permeability structure. Geology 24, 1025–1028. Caine, J.S., Forster, C.B., 1999. Fault zone architecture and fluid flow: insight from field data and numerical modeling. In: Haneberg, W.C., Mozley, P.S., Moore, J.C., Goodwin, L.B. (Eds.), Faults and Subsurface Fluid Flow, Geophysical Monograph 113. AGU, pp. 101–128. Cavanagh, A., Ringrose, P., 2011. Simulation of CO2 distribution at the In Salah storage site using high-resolution field-scale models. Energy Procedia 4, 3730–3737. Ciotoli, G., Lombardi, S., Annunziatellis, A., 2007. Geostatistical analysis of soil gas data in a high seismic intermontane basin: Fucino plain, Central Italy. Journal of Geophysical Research, 112-B05407. COMSOL. 4.2, 2011. Model Documentation: Discrete Fracture, www.comsol.com/ showroom/documentation/model/691/ Dershowitz, W.S., Fidelibus, C., 1999. Derivation of equivalent pipe network analogues for three-dimensional discrete fracture networks by the boundary element method. Water Resources Research 35 (9), 2685–2691. Gilmour, H.M.P., Billaux, D., Long, J.C.S., 1986. Models for calculating fluid flow in randomly generated three-dimensional networks of disk-shaped fractures: theory and design of FMG3D, DISCEL and DIMES. Lawrence Berkeley Laboratory, Earth Sciences Division, Berkeley, CA, LBL-19515, 144 pp. Konzuk, J.S., Kueper, B.H., 2004. Evaluation of cubic law based models describing single-phase flow through a rough-walled fracture. Water Resources Research 40 (2), W024021–W02402117. La Bella, R., 2006. La modellazione tridimensionale di reticoli di faglie e fratture negli studi di reservoir. Un nuovo approccio applicato ad un analogo naturale. Tesi Sperimentale in Geologia Strutturale, A. A. 2005/2006. Università degli Studi di Roma “La Sapienza”, pp. 86 pp. La Pointe, P.R., Hudson, J.A., 1985. Characterization and interpretation of rock mass joint patterns. Geological Society of America Special Paper 199, 37 p. Lanaro, F., 2000. A random field model for surface roughness and aperture of rock fractures. International Journal of Rock Mechanics and Mining Sciences 37 (December (8)), 1195–1210. Lapcevic, P.A., Novakowski, K.S., Sudicky, E.A., 1999. The interpretation of a tracer experiment conducted in a single fracture under conditions of natural groundwater flow. Water Resources Research 35 (8), 2301–2312.

199

Long, J.C.S., Aydin, A., Brown, S.R., Einstein, H.H., Hestir, K., Hsieh, P.A., Myer, L.R., Nolte, K.G., Norton, D.L., Olsson, O.L., Paillet, F.L., Smith, J.L., Thomson, L., 1996. Rock Fracture and Fluid Flow: Contemporary Understanding and Applications. National Academy Press, Washington, DC, pp. 551 pp. Maerz, N.H., Franklin, J.A., Bennett, C.P., 1990. Joint roughness measurement using shadow profilometry. International Journal of Rock Mechanics and Mining Sciences 27 (5), 329–343. Marrett, R., Ortega, O.J., Kelsey, C.M., 1999. Extent of power law scaling for natural fractures in rocks. Geology 27, 799–802. Masini, M., Poblet, J., Bulnes, M., 2010. Cross section restoration: a tool to simulate deformation. Application to a fault-propagation fold from the Cantabrian fold and thrust belt, NW Iberian Peninsula. Journal of Structural Geology 32, 172–183, http://dx.doi.org/10.1016/j.jsg.2009.11.002. Mazzoli, S., Szaniawski, R., Mittiga, F., Ascione, A., Capalbo, A., 2012. Tectonic evolution of Pliocene-Pleistocene wedge-top basins of the southern Apennines: new constraints from magnetic fabric analysis. Canadian Journal of Earth Sciences 49 (3), 492–509, http://dx.doi.org/10.1139/e11-067. Neuman, S.P., 2005. Trends, prospects and challenges in quantifying flow and transport in fractured rocks. Journal of Hydrology 13, 124–147. Oda, M., 1985. Permeability tensor for discontinuous rock masses. Geotechnique 35, 483–497. Ortega, O.J., Marrett, R.A., Laubach, S.E., 2006. A scale-independent approach to fracture intensity and average spacing measurement. AAPG Bulletin 90 (2), 193–208. Palladino, D.M., Simei, S., 2005. The Latera Volcanic Complex (Vulsini, central Italy): eruptive activity and caldera evolution. Acta Vulcanologica, Special Issue 17 (1–2), 75–80. Pettinelli, E., Beaubien, S.E., Zaja, A., Menghini, A., Praticelli, N., Mattei, E., Di Matteo, A., Annunziatellis, A., Ciotoli, G., Lombardi, S., 2010. Characterization of a CO2 gas vent using various geophysical and geochemical methods. Geophysics 75, B137–B146, http://dx.doi.org/10.1190/1.3420735. Priest, S.D., Hudson, J.A., 1981. Estimation of discontinuity spacing and trace length using scan line surveys. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 18, 183–197. Putra, E., Fidra, Y., Schechter, D.S., 1999. Use of experimental and simulation results for estimating critical and optimum water injection rates in naturally fractured reservoirs. In: Proceedings – SPE Annual Technical Conference and Exhibition: ‘Reservoir Engineering’, Houston, TX, USA; October 1999; Volume 1, 1999, Code56642. Ringrose, P.S., Roberts, D.M., Gibson-Poole, C.M., Bond, C., Wightman, R., Taylor, M., Raikes, S., Iding, M., Ostmo, S., 2011. Characterization of the Krechba CO2 storage site: critical elements controlling injection performance. Energy Procedia 4, 4672–46789, http://dx.doi.org/10.1016/j.egypro.2011.02.428. Sibson, R.H., 1996. Structural permeability of fluid-driven fault-fracture meshes. Journal of Structural Geology 18, 1031–1042. Snow, D.T., 1969. Anisotropic permeability of fractured media. Water Resource Research 5 (6), 1273–1289. van Dijk, J.P., 1998. Analysis and modeling of fractured reservoirs. In: European Petroleum Conference 1, pp. 31–43. van Dijk, J.P., Bello, M., Toscano, C., Bersani, A., Nardon, S., 2000. Tectonic model and three-dimensional fracture network analysis of Monte Alpi (southern Apennines). Tectonophysics 324, 203–237. van Dijk, J.P., 2002. Data Driven Fracture Models. SFERA Inaugural Meeting 2002 , Pescara (Italy), pp. 41–45, Abstracts Volume, Paper 8. van Dijk, J.P., 2013. Lo stoccaggio della CO2 nel sottosuolo Il ruolo delle Geoscienze. Atti del 1◦ Congresso dell’Ordine dei Geologi di Basilicata. In: “Ricerca, Sviluppo ed Utilizzo delle Fonti Fossili: Il Ruolo del Geologo”, Potenza, 30 Novembre – 2 Dicembre 2012, pp. 29–76. Walsh, J.J., Watterson, J., Heath, A.E., Child, C., 1998. Representation and scaling of faults in fluid flow models. Petroleum Geosciences 4, 241–251. Witherspoon, P.A., Wang, J.S.Y., Iwai, K., Gale, J.E., 1980. Validity of Cubic Law for fluid flow in a deformable rock fracture. Water Resources Research 16 (6), 1016–1024, http://dx.doi.org/10.1029/WR016i006p01016. Wu, H., Pollard, D.D., 1995. An experimental study of the relationship between joint spacing and layer thickness. Journal of Structural Geology 17, 887–905. Zhang, L., Einstein, H.H., Dershowitz, W.S., 2002. Stereological relationship between trace length and size distribution of elliptical discontinuities. Geotechnique 52 (6), 419–433. Zhou, X., Du, J., Chen, Z., Cheng, J., Tang, Y., Yang, L., Xie, C., Cui, Y., Liu, L., Yi, L., Yang, P., Li, Y., 2010. Geochemistry of soil gas in the seismic fault zone produced by the Wenchuan Ms 8.0 earthquake, southwestern China. Geochemical Transactions 11 (5), http://www.geochemicaltransactions.com/content/11/1/5 Zinszner, B., Pellerin, F.M., 2007. A Geoscientist’s Guide to Petrophysics. IFP Publications Edition Technip, Paris. Zito, G., Mongelli, F., de Lorenzo, S., Doglioni, C., 2003. Heat flow and geodynamics in the Tyrrhenian Sea. Terra Nova 15 (6), 425–432.