Many key aquifers and oil reservoirs are in carbonate rocks. Understanding the flow behavior within this commonly complex pore space requires new perspectives and technology in order to improve carbonate aquifer and reservoir characterization. Dissolution of carbonates is related to flow; hence, quantifying the size of dissolution vugs on carbonate outcrops can help characterize controls on flow, namely matrix permeability and fracture connectivity. LIDAR (light detection and ranging) scans, combined with high-resolution photography, enable us to both measure vugs' areas and assess spatial relationships between vugs, beds, and fractures. We developed a method of obtaining and interpreting necessary vug, bed, and fracture data on the basis of these technologies. Application of this method on a Cretaceous Edwards Group outcrop in Texas (United States) revealed a significant correlation between the relative vug area of beds obtained remotely and air permeability measured in plugs extracted from these beds (R2 = 0.94, P = 0.001). The total area of vugs intersected by fractures was used to establish a probability density function of fracture lengths that can improve flow modeling of the reservoir. These findings show the potential of using LIDAR and photo images in reservoir characterization by data analysis of geological features, in addition to their use for accurate mapping.


Carbonate oil reservoirs hold significant quantities of the world's oil and gas reserves (Al-Anzi et al., 2003), and carbonate aquifers supply a great deal of fresh water globally (Ford and Williams, 1989). Nevertheless, modeling carbonate fluid reservoirs is complicated because of the spatial heterogeneity in their reservoir properties. These reservoirs are often discussed in terms of triple porosity and/or triple permeability (White, 2002; Camacho-Velázquez et al., 2005; White and White, 2005). Matrix, fractures, and dissolution-enlarged conduits are the three types of media in which fluids can flow within carbonate rock. The variability of permeability between and within each of these, and even within a single feature—i.e., a bed, a fracture, and a dissolution channel—can be very large (Nordqvist et al., 1996; Jennings et al., 2000). These heterogeneities add up to complications in modeling these reservoirs because of the discrete characteristics of fractures and conduits (networks) embedded in the matrix continuum, as well as the nature of fluid and contaminant exchange between them. Nevertheless, in this work we propose a conceptual model and demonstrate a remotely sensed data collection and analysis procedure that help characterize the aforementioned variability.

Conceptual Model for Fractured Carbonate Reservoirs

Carbonate reservoirs are viewed differently in different disciplines. In the petroleum geology and engineering communities, subcentimeter-scale reservoir properties of the carbonate matrix are of interest (e.g., Lucia, 1995; Lønøy, 2006), as well as reservoir-scale structured heterogeneity (e.g., Botton-Dumay et al., 2002). In hydrologic literature on carbonate aquifers, the matrix is often neglected, while development of hundred-meter- to kilometer-scale conduits by dissolution in fractures and fissures has drawn focus (e.g., Dreybrodt, 1996; Hanna and Rajaram, 1998; Kaufmann and Braun 2000; Gabrovsek et al., 2004). This difference reflects differences in the importance of matrix-conduit relations in the two porous-media fluid reservoir disciplines. Whereas oil production from a closed (not recharging) reservoir requires “squeezing out” the oil from the less permeable yet high-storage matrix, for sustainable production of water from a recharging carbonate aquifer, mobilization of the relatively stagnant water in the matrix is of lesser interest. We suggest a generalized conceptual view of carbonate-fluid reservoirs that can utilize remotely sensed observation of outcrops to improve reservoir characterization.

Figure 1 shows two carbonate outcrops of very different settings; one is from a Cretaceous limestone observed in a roadcut 45 km northwest of San Antonio, Texas, USA, and the other is from an Eocene chalk observed on an ephemeral stream bank 15 km south of Beer Sheva, Negev Desert, Israel (Kurtzman et al., 2005). However, a common characteristic relationship between the horizontal and continuous features, i.e., beds and bed partings, the steeply dipping fractures, and the dissolution features, hereafter termed vugs, can be observed in both outcrops. Vugs tend to develop at some intersections of steeply dipping fractures and horizontally continuous and relatively permeable features.

In Kurtzman et al. (2007a), it was shown theoretically that a positive linear relationship exists between dissolution intensity and flow rate. This positive relation was demonstrated by similarity between the spatial distribution of fluxes obtained by flow simulations and location of vugs on the physical outcrop (Kurtzman et al., 2007a). The autocatalytic process, in which enhanced dissolution caused by higher flow rates brings about a further increase in flow rates owing to enhanced permeability in the dissolution area, leads to concentration of flow in discrete features. Dissolution vugs at intersections of flowing fractures and permeable beds (or bed partings) are an outcrop manifestation of these discrete features (Fig. 1). The aforementioned process and observation of vugs and conduits at different scales are bases for the conceptual model of fractured carbonate reservoirs presented in the next paragraphs.

The conceptual model assumes that the reservoir consists of layered beds that are relatively continuous horizontally. These beds differ in their porosity and permeability because of their depositional circumstances, which lead to variation in grain size. The carbonate block that is now a closed reservoir or a recharging aquifer at one time underwent faulting and cracking. At some point an aggressive (with respect to calcium-carbonate dissolution) water source was connected by a fracture (or a fracture network) to a sink, enabling a continuous flow of aggressive fluid into the fracture, and fluid with higher concentrations of calcium out of it. The process of aperture and flow-rate enhancement in this fracture has been covered in depth in the hydrologic literature on carbonate aquifers (e.g., Dreybrodt, 1996; Hanna and Rajaram, 1998; Kaufmann and Braun, 2000). This growing preferential flow path attracts lateral flow from its surroundings. In the more permeable lateral features (permeable beds or bed partings), this attraction causes relatively high flow rates toward the fractures. Flow rates in fractures near these features therefore increase more than elsewhere, and hence dissolution occurs, leading to the observation in Figure 1.

On a larger scale, the most extensive faults (or fault zones) have higher chances of becoming preferential pathways due to their longer lateral extent (higher chances of connecting a source to a sink, directly or through a network connection) and higher probability for larger initial apertures. These features will attract lateral flow from even farther away, making it more probable that the entire fault zone will develop larger vugs at all scales.

Going up to the aquifer scale, many carbonate aquifers discharge their waters through one or a few high-flow-rate springs. Equivalent porous medium models that are often used to model these aquifers are calibrated using low hydraulic conductivities in the outcrops and increasing conductivities farther downstream toward the spring area (e.g., see Scanlon et al., 2003, for the Barton Springs Edwards aquifer, Texas). This is a manifestation of the process of concentrating flow to a narrow zone, increasing flux, increasing dissolution, and increasing drainage capability of the downstream zone, which is generally the same process that can be observed in vugs on outcrops similar to those shown in Figure 1.

On the subcentimeter scale (core-plug scale), the higher flow rate that develops within the permeable beds also causes small-scale dissolution vugs (submillimeter) to form preferential pathways within the matrix and, therefore, larger heterogeneity in plug permeabilities within a permeable bed. During evolution of the reservoir some flow may carry significant quantities of fine grains that can be deposited and cemented to the aforementioned submillimeter vugs, thus clogging and reducing core-plug-scale permeability below predissolution interparticle permeability. Dolomitizing and/or calcitizing fluid flows can alter these bed permeabilities as well. Such nondissolution processes can cause some zones within the permeable beds to have very low porosity and permeability.

The phenomenon of concentrating flow and enhancing heterogeneity caused by the dissolution of carbonate that influences permeability distribution at all scales enables us to draw conclusions both about average plug permeability and about structures on the scale of hundreds of meters, using observation of vugs at the 5 cm to 10 m scale. An efficient way to collect data on the three types of features (beds, fractures, and vugs) from reservoir-analog and/or aquifer outcrops is by using remotely sensed laser scans combined with high-resolution photos.

LIDAR, Photos, and Combined Images and Their Use in Earth Sciences

Scanning light detection and ranging (LIDAR) is a technology that has come to the forefront of surveying in the past 10 years. This laser-based measurement system allows rapid acquisition of detailed point data describing a terrain surface, from both aerial and terrestrial platforms (Bellian et al., 2005; Buckley et al., 2008). Use of this technology in earth sciences is growing rapidly.

The raw product of ground-based LIDAR is a point cloud of x, y, z coordinates and intensity of reflection from the point that was hit by the laser beam. Any number of scans, taken from the same and/or different points and having a significant overlap, are then merged so that a single coordinate system can be acquired (Bellian et al., 2005; Bonnaffe et al., 2007). Once the point data from each individual acquisition station are rectified and merged into a single spatially coherent point cloud, a triangulated surface is fit through the raw points to produce a triangular irregular network (TIN) model (Bonnaffe et al., 2007). Features like fracture traces or bed contacts may lack a strong enough topographic expression to be recognized directly on the point cloud or TIN. These features may be recognized on color photographs; therefore a high-resolution digital photograph is “draped” on the TIN surface through registration of selected points to form a photorealistic model (Xu et al., 2000; Alfarhan et al., 2008). The current work focuses on outcrop data analysis that is based on a photorealistic model of a carbonate outcrop and reservoir properties that can be estimated from it rather than the technical details involved in construction of the photorealistic model.

Schmugge et al. (2002), reviewing remote sensing in hydrology, acknowledged both large-scale aerial LIDAR scans for measuring channel and floodplain cross sections, as well as millimeter-resolution scans for micro-topography, influencing soil-water processes of interest in agriculture. An overview of LIDAR and its application to outcrop studies was provided in Bellian et al. (2005). Olariu et al. (2008) used laser scans to investigate fracture sets in a deep-water siliciclastic deposit. Labourdette and Jones (2007) studied elements of fluvial deposition sequences using a combination of LIDAR and aerial photography. Enge et al. (2007) presented a detailed workflow, starting with outcrop selection, data collection, construction of the virtual outcrop (similar to the previously mentioned photorealistic model of the outcrop), construction of a grid, and introducing offsets on faults according to the virtual outcrop data, attaching reservoir properties and running flow simulations. The current study also intends for outcrop-scanned data to be used in improving reservoir characterization, but rather than simulating flow through the outcrop, the scanned data are used to improve characterization of bed permeability and fracture length distribution—factors controlling flow and transport in the carbonate reservoir.


Site and Data Acquisition

A carbonate roadcut through the Cretaceous Edwards Group located near Lake Medina, ~45 km west of San Antonio, Texas (N29°33′.68, W98°54′.29), was selected for the study (Fig. 2A). The Edward's aquifer average annual recharge (A.D. 1934–2006) is ~900 × 106 m3. This aquifer serves approximately two million people in south-central Texas. The roadcut exposes the Dolomitic Member of the Kainer Formation in the lower section of the Edwards Group just above the confining Upper Glen Rose Formation (see stratigraphic chart in Fig. 2B). The Dolomitic Member consists of mudstone to grainstone and crystalline limestone with chert. In this specific outcrop, most beds consist of limestone, with fabrics ranging from wackestone to grain-dominated pack-stones. Some beds have been partly dolomitized and recalcitized. Beds are 30–150 cm thick and dip very slightly to the north (02°). Two major sets of fractures intersect these beds. Fractures from both sets strike between N45°E and N80°E but differ in their equivalent dips, which are 60°S and 60°N on the west wall and 70°S and 70°N on the east face (fractures of the two sets can be observed in Fig. 1A). A few faults of the aforementioned fracture orientations are observed in the roadcut; the largest vertical offset associated with these faults is 75 cm. Both east and west faces of the roadcut were laser scanned and photographed in order that a photorealistic model could be constructed.

LIDAR data were acquired using an Optech ILRIS-3D laser ranging and imaging system at average point spacing of 0.9 cm. A total of 33,006,544 points from 32 scans at 9 different survey positions was acquired on the west wall and 21,851,017 points from 26 scans at 8 different survey positions on the east wall of the road-cut. Alignment and processing of the data were performed using the Polyworks (Innovmetric Software Inc.) software suite. The merged laser point cloud was registered within Polyworks IMAlign. Optimized TIN surfaces of 7,870,341 and 7,094,058 triangles for the west and east faces, respectively, were generated from the original point clouds using Polyworks IMMerge software. This process removed excess triangles within specified tolerance on the basis of deviation from surrounding triangle normals. Photographs with a resolution of ~0.4 cm were taken by a Kodak DCS SLR/c 14-MP camera with a 35 mm lens. Photo registration transformation was outsourced to Xueming Xu at Real Earth Models, Inc. Registration that accounts for camera distortions (Xu et al., 2000) was established through 12 control points from each side of the roadcut to form photorealistic models of the 2 faces (Fig. 2B). Digitization and analysis of geological features were conducted using Poly-works IMInspect software.

Digitizing Bed Contacts, Fractures, and Vugs

Exposed fracture surfaces and distinct bed contacts (e.g., marly beds) can be directly mapped on the LIDAR point clouds, whereas less distinct bed contact and fracture traces that have no visible topographic imprint on the ~1 cm resolution point cloud are mapped on the basis of their distinctive color change on the photorealistic model. Figure 3 shows digitized bed contacts and fractures on the photorealistic model versus those digitized on a point cloud.

Digitizing vugs was more time consuming than digitizing fracture traces. The vug was generally first identified by dark shading on the photorealistic model, but its boundaries could be defined only by rotating the three-dimensional (3D) TIN model and finding where the rock surface bends in to form the vug (Fig. 4).

Estimating Vug Depths, Volumes, and Areas from LIDAR Data

Vug depth was defined as the largest distance between a point within the TIN representation of the vug and a plane that was best fit to the points of the vug's 3D closed polyline (Figs. 4 and 5). Vugs that had no measurable depth were discarded. The limited resolution of the LIDAR data and especially the decimated TIN combined with the not-always-orthogonal-to-the-face scanning resulted in an underestimation of topographical vug depth.

Figure 6 shows the relationship between vugs' depths measured in the field as the largest depth to a solid rock and those estimated in the procedure described in Figure 5. A correction factor of 2.4 (slope of the best fit line; Fig. 6) was applied to all vugs with remotely sensed depth <1 m. No correction was made for the few large caves in which LIDAR depth measured was >1 m, assuming that LIDAR penetration is good in these wide openings. Whether these vug depth estimates are useful for our purposes is discussed later.

Vug volume was calculated in Polyworks (surface-to-plane volume under the measure menu) by using the same plane that was used for depth. The factor of 2.4 was applied to these volumes, assuming that the average depth under estimate calculated for the largest depth of each vug is our best estimate for the underestimate for all points in the vug and therefore for the volume.

Because the vug area of a 3D closed polyline digitized in Polyworks is not well defined, vug areas of the 2D projection (on the Y [south-north] and Z [vertical]) of the closed polylines were used to form 2D vug polygons (Fig. 7). A correction factor of 1/cos(16.5°) = 1.043 was used to correct for the average trend of the outcrop relative to the north. The 2D projection was established by exporting the digitized data to text format and importing the Y and Z coordinates as XY data into ArcGIS (ESRI Inc.) software. Advantages of the latter software for further analysis are discussed later.

From Vug Geometry to Reservoir Properties

A linear relationship between a bed's plug permeability (geometric mean of one to five plugs per bed) and its relative vug area (RVA), defined in equation 1, was shown in Kurtzman et al. (2007a).  
where the index i is the bed number, Aj is the area of vugj, ni is number of vugs surveyed in bedi, and Bi is the area of bedi in which the survey was conducted. Vug area was estimated in that work (Kurtzman et al., 2007a) using horizontal and vertical manual tape measurements of vugs that can be approached from the road only on the east side of the same outcrop. In Kurtzman et al. (2007a), it was assumed that vugs that contain recent soil are flow-related vugs (FRV) rather than molds. On the basis of this assumption, they suggested a logistic regression model to help discriminate FRVs from molds for remotely sensed vugs in which the presence or absence of soil cannot be observed (equation 2).  
where S = [(vug area)0.5]/[vug depth], F—whether the vug was intersected by a fracture (1) or not intersected (−1) and V is the vug volume (cm3). Pr(FRV) is the probability that the vug is an FRV; in this work if Pr(FRV) > 0.5 (i.e., the logistic regression model estimates that the vug is an FRV with more than 50%), the vug is considered an FRV.

Using the RVA of the FRVs instead of all the surveyed vugs increased the correlation with plug permeability for manually measured vugs (Kurtzman et al., 2007a). The ratio of vug area to vug depth (S) was the most significant predictor in this logistic regression model, whereas vug volume had minor significance. In the current work we test how well these relationships hold up when outcrop and vug geometries are surveyed remotely using a LIDAR and a camera.

In Kurtzman et al. (2007a), the effect of fracture connection to a water source and sink on flow through the fracture and, hence, dissolution along it, was demonstrated through flow simulations. The longer a fracture extends laterally, the higher the chances are that it will be connected directly or through a fracture network to a flow passage. However, a fracture's lateral length (fracture length hereafter) cannot be measured on outcrops. Fracture-length distribution is also unknown, and the significance of estimations that are based on geometrical properties of fracture-trace data from outcrops is questionable (Kurtzman et al., 2007b). Fracture-length distribution is important when large numbers of fractures are modeled for a flow model of a fractured reservoir using a stochastic fracture generator (Chile's and de Marsily, 1993; Dershowitz et al., 1998). We therefore suggest basing the estimation of fracture-length distribution (for flow-modeling proposes) on the distribution of dissolution-vug area along fracture traces, rather than on trace length or other geometrical properties of traces.

Relating Vugs to Beds and Fractures

After the fracture traces, the bed contacts and vugs were digitized and vug geometry was estimated, as described earlier. The type of analysis suggested in the previous discussion (from vug geometry to reservoir properties) requires arranging the fracture, bed, and vug data in a tabular form to enable sorting, selection, and statistics. Another capability of the software that can help in the type of analysis suggested in the same discussion is the topological relationship between spatial features, which will enable spatial queries and selections (e.g., “Select all vugs whose center points [centroids] are within bed #2,” or “Sum the areas of all vugs that are intersected by fracture #44.”). Figure 8 shows these capabilities within the ArcGIS spatial database. Each spatial feature (vug, fracture trace, bed) is related to the database through feature identification (FID—green numbers for vugs, black for fractures; Fig. 8). Besides the FID and shape mandatory fields, we can add as many text and numerical fields as we wish. In this case, some fields were calculated using topological relations between the three layers of the spatial data (beds, fractures, and vugs). For example, for the selected vug (FID = 79, lower table), the field “Bed_center” is the answer to the spatial query “In what bed does the centroid of the polygon lie?” (for this specific vug, this is bed #2 [red numbers]). The field “Fracture” is 1 if the vug is intersected by a fracture. For the fracture selected in Figure 8 (FID = 44, upper table), a spatial-join operation resulted in a count of the vugs intersected by the fracture (the field “Count_”) and the sum of these vugs' properties (area, depth, and volume, from which only the sum of areas is necessary for our analysis).

Although ArcGIS is a platform that has these capabilities, its 3D capabilities are limited. Nevertheless, when dealing with outcrops in which the third dimension (perpendicular to the outcrop face) is relatively small, moving to 2D is beneficial in gaining database and topology capabilities. Results of the previous 3D analysis such as vug depth and volume can still be used, even though these values will have no geometrical representation. They will only be numbers in a table related to a 2D vug polygon.


Summary of the Features

Table 1 and Figure 9 summarize numerically and graphically the geological features that were digitized on the 3D photorealistic model and transferred to the geographic information system (GIS) spatial database.

Correlation between Bed Permeability and Its Scanned RVA

A distinction was made between bed-related vugs and larger vugs whose area intersects two or more beds, making it less easy to relate them to a specific bed. These larger vugs (caves) will be used later for characterizing the fractures to which they are related. Only the smaller, bed-related vugs were used for checking the correlation between beds' relative vug areas (RVA, equation 1) and air permeability measured in random ~2.54 cm (1 inch) plugs extracted from beds 1, 2, 3, 4, 5, and 7 on the east face (yellow and blue beds in Fig. 8 numbered from the bottom up). The largest bed-related vug is that with FID 79 that was selected in Figure 8, with an area of 0.61 m2, and the following analysis is based on vugs with areas equal to or smaller than this threshold. Some of the larger (and therefore important) bed-related vugs are not fully contained in the bed to which they are related, especially in permeable beds 2 and 6. Nevertheless, it is obvious that they are related to these permeable beds (e.g., the lower, grayish bed with vugs in Fig. 1A). Therefore, two criteria for assigning vugs to a bed were used in the spatial database, and the analysis derived from them was compared (Table 2; Fig. 10). A vug was assigned to the bed in which (1) its centroid point is contained, and (2) a point located at the vug center, horizontally, and at one-eighth, or 0.125 of the vug's height above the vug's lowest point, vertically (termed eighth hereafter).

The total number of vugs related to beds 1, 2, 3, 4, 5, and 7 from which plug permeabilities were available was 153 when vugs were assigned to beds with the eighth criterion and 157 if the centroid criterion was used. The number of vugs from each bed for which their area was summed in order to calculate beds' RVA is given in Table 2. Of the 153 vugs, 95 vugs follow the geometric-statistical criteria Pr(FRV) > 0.5 when using the logistic regression model (equation 2). Figure 10D shows the correlation obtained when only this subset of vugs was used.

The use of a point lower than the vug's centroid for automatically relating a vug to a bed resulted in slightly better correlation than that obtained using the centroid (Figs. 10A, 10B), because some large vugs associated with the more permeable beds (2 and 6) develop near the boundary of these beds and the less-permeable beds above.

Correlation between permeability and relative vug volume of a bed (RVV, defined as the sum of the bed's vug volumes, divided by the bed's area) is smaller than the correlation with RVA (cf. Figs. 10C and 10B). One explanation for this lower correlation is that the physical process that leads to the correlation is the linear relationship between dissolution intensity and flow (Kurtzman et al., 2007a) and the linear relationship of flow and permeability (Darcy's Law). Whereas vug area is dictated by the dissolution process, the remotely scanned vug volume is controlled largely by filling of the vug with recent soil, which is not part of the dissolution process. Another reason for the correlation of permeability with vug volumes being smaller than that of permeability and vug areas is due to the poor estimation of vugs' depth and volume obtained by LIDAR scans (Fig. 6). This is also the reason that implementation of the logistic regression model for discriminating flow-related vugs (equation 2) has not improved the correlation here when vug geometry was obtained remotely (cf. Figs. 10B and 10D), as it did in the case of manually surveyed vugs (Kurtzman et al., 2007a). Green points in Figures 10A, 10B, and 10D are the average permeability of bed 6 estimated by its RVA and the red line regression models. As expected, the estimated permeability of this grainier fabric bed is significantly higher than that found in muddier beds 1, 3, 5, and 7.

Estimation of Fracture-Length Distribution from Fracture and Vug Data

The spatial database of fractures built from remotely sensed data contains 126 fracture traces, in which the sum of areas of the vugs that they intersect is easily obtained using spatial selection and joining tools in ArcGIS (Fig. 8). Figure 11 shows a histogram of these fractures binned by the normalized sum of areas of the vugs that they intersect. Following the suggestion in the methods section, we maintain that for simulating fractures in the studied or analogous formation, to be inserted as high-permeability features in a flow model, fracture-length distribution can be estimated using this histogram. Vug areas in this histogram (Fig. 11) were normalized by the area of the largest vug recorded at the site, which is the shaft along the major fault on the east face (Figs. 2B and 9). The rationale behind this normalization of the vug area to the largest vug on site is based on the multiscale occurrence of dissolution and flow concentration to the dissolved area (described in the introduction). The estimated normalized histogram of fracture length can therefore be used at different scales. This can be done by assigning the fracture (fault) on which the largest dissolution vug develops a length on the order of the simulation domain, and other fracture lengths will be proportional to the relative area of the vugs formed along them.

A probability density function (PDF) was fitted to the histogram in Figure 11 (blue lines and text). This PDF is divided into three intervals: 0 < x < = 0.01, 0.01 < x < = 0.1, and 0.1 < x < = 1. The three integrals of f(x) over these three intervals result in 0.873, 0.119, and 0.008, respectively, which add up to 1, as required from a PDF. Using this PDF for relative fracture lengths, as suggested above, implies that fewer than 1% (uniformly distributed) of the fractures will have a length on the same order of magnitude as that of the simulation domain; ~12% (uniformly distributed) will have a length one to two orders of magnitude shorter than the length of the simulation domain; and ~87% (exponentially distributed) will be two orders of magnitude shorter or less.

Beyond the Single Fracture and Bed Scale

Usually we will not be able to model down to the resolution of a single fracture. However, the fracture and vug spatial database can help us conceptualize important features of higher reservoir scale (Fig. 12). There is evidence that the orientation of the major fault is the crucial orientation in this block of rock, and that some high-permeability stripes in this orientation can be indicated by long fracture traces and large vugs on both sides of the roadcut. Between these sections, the concentration of fractures and vugs is smaller and/or high only on one side of the roadcut (Fig. 12). This type of interpretation can help in well location, spacing, and other operational decisions.


Dissolution of carbonate rocks is a major source of heterogeneity in reservoir properties of carbonate fluid reservoirs at all scales. Dissolution intensity is related fundamentally to flow rate through the carbonate rock, which in turn is closely related to fracture connectivity and bed matrix permeability. Dissolution vugs observed on carbonate outcrops are a manifestation of dissolution intensity; therefore, their spatial relationships with fractures and beds can help us model the aforementioned reservoir properties.

Combining LIDAR scanning and high- resolution photography to a photorealistic model enables accurate digitization of vugs, beds, and fractures over large outcrops. Exporting the digitized data to a GIS makes it easy to assess vug-bed and vug-fracture spatial relations.

A significant correlation (R2 = 0.94, P = 0.001) between (1) relative vug area of beds obtained by analyzing a photorealistic model of a roadcut cutting through the lower section of the Edwards Group in south-central Texas, and (2) air permeability measured on 2.54 cm plugs extracted from the different beds was observed. Vug depth and volume estimations based on the remotely sensed data were not good enough. On the other hand, vug-area estimates were accurate and served as a good proxy for dissolution intensity. Statistical discrimination of flow-related vugs and molds did not improve the correlation with permeability for the remotely sensed data, as it did when vug-depth measurements were obtained physically on the outcrop (Kurtzman et al., 2007a).

Flow-through fractures depend on connectivity of the fracture to a high-pressure source and low-pressure sink. The linear relationship between dissolution intensity and flow rate, the fact that there is no good way to estimate the fracture's lateral length from its trace, and the influence of fracture length on its connectivity led us to suggest a method for evaluating fracture-length distribution on the basis of the area of the vugs that they intersect. A PDF of relative fracture length for simulating stochastic fractures was fitted to the fracture and vug data extracted from the photorealistic model of the Edwards outcrop. This PDF indicates that ~1% (uniformly distributed) of the fractures should have a length on the order of the length of the simulation domain; 12% of the fractures (uniformly distributed) should have lengths between one and two orders of magnitude smaller; and 87% (exponentially distributed) should have a length between zero and two orders of magnitude smaller than that of the simulation domain.

The need for improved characterization of carbonate reservoirs is rising as the demand for fluids extracted from these reservoirs increases. The current work shows how advanced remote sensing technology can be used, together with an understanding of the role of dissolution in these reservoirs, for improving reservoir characterization.

This research was performed by the Reservoir Characterization Research Laboratory (RCRL). We are grateful to all the RCRL industrial associates for their funding. Kurtzman thanks the Volcani Center for time given to write the paper. Support for this work was provided in part by the John A. and Katherine G. Jackson School of Geosciences and the Geology Foundation at The University of Texas at Austin. Publication authorized by the Director, Bureau of Economic Geology, The University of Texas at Austin.