Abstract

Determination of fracture orientation can be an important aspect of structural analysis in reservoir characterization. The availability of ground-based laser scanner systems opens up new possibilities for the determination of fracture surface orientation in rock outcrops. Scanners are available in low-sample-density, low-accuracy, and fast, high-sample-density, high-accuracy models. These automatic laser scanner systems produce enormous volumes or “clouds” of point data at an instrument-dependent accuracy and resolution, which can be at the millimeter level.

This huge volume of data calls for an automated and objective method of analysis. We have developed a surface classification algorithm based on a multipass partitioning of the point cloud. The method makes use of both spatial proximity and the orientation of an initial coarse-grained model of the point cloud. Unsupervised classification of surface sets is demonstrated herein using the new algorithm.

Both previously mentioned types of scanners have been used to map the Jackfork sandstone outcrop at Big Rock Quarry in Little Rock, Arkansas. We apply the surface classification algorithm to these data to extract fracture surface orientations from the point cloud. The effectiveness of these new technologies when applied to fracture analysis is clearly demonstrated in this example. It is also shown that the low-density, low-resolution type of scanner is adequate to define general geomorphology but is inadequate for fracture definition. The surface classification algorithm can be used to reliably extract fracture and bedding strike and dip angles from the three-dimensional point locations acquired using centimeter-accurate, high-density laser scanner systems.

INTRODUCTION

Studies of natural fractures in outcrops are important because fractures can reveal the strain history of the host rock and they are economically important as paths for hydrocarbon and water flow. Since production can be improved or impeded depending on the placement of wells relative to the fracture orientations, detection and modeling of fractures can be an important part of reservoir characterization. Although fractures can be studied using borehole and geophysical data, this paper will focus on outcrop-based methods.

Typical methodology involves a detailed survey of fractures, in outcrop, by visual observation. This technique, although precise, is labor and time intensive and is hampered by the limited amount of data obtained, especially from large exposures, to be used for statistical analysis. Three-dimensional (3-D) geological structures such as faults, folds, fractures, and deformed sedimentary strata can be mapped using the global positioning system (GPS). This data-capture strategy consists of walking along faults, fractures, exposed contacts between formations, as well as other significant stratigraphic markers (Maerten et al., 2001), which is time consuming and effort demanding. Although fracture systems occupy three-dimensional volumes, they are most often studied in one-dimensional transects or two-dimensional maps and cross sections.

The approach of this paper is to use 3-D point clouds from ground-based laser scanners to identify fracture orientations in outcrops. Ground-based laser scanners have proven very useful for acquisition of high-quality, high-resolution, three-dimensional terrain data from outcrops (Xu et al., 2000; Bellian et al., 2002, 2005; McCaffrey et al., 2005). Surface geometry of the outcrop is remotely captured in the form of dense 3-D point clouds. Laser data have been used in other studies for 3-D modeling and mapping techniques as well as for fracture and fault recognition (Ahlgren et al., 2002; Ahlgren and Holmlund, 2003). The enormous volume of data produced by a laser scanner system demands an automated and objective processing technique. We have developed an unsupervised pattern classification technique that can extract surface orientation information directly from the point cloud. The result is the equivalent of taking hundreds of conventional strike and dip measurements in even the most inaccessible parts of an outcrop. This new surface classification algorithm is described in detail and demonstrated on synthetic data and on the Jackfork outcrop.

The methodology described in this paper provides a direct and easily interpretable visualization of the location, orientation, and geometry of zones of fracturing. Since outcrop geometry helps in identification and interpretation of fracture zones, our 3-D analysis of terrain surfaces of the outcrop at Big Rock Quarry shows how fractures with different orientations are distributed within the rock volume. Direct visualization of the 3-D fracture planes provides complementary information to the traditional stereonets, rose diagrams, and histogram plots.

METHODOLOGY

Outcrop Data Collection and Processing

Terrestrial laser scanner and differential GPS systems are now robust and portable. The availability of long-range, 3-D laser imaging systems offers the possibility of rapid precise mapping to distances of up to 1 km from a set-up point. In this study, several different scanner systems were used and compared. In these systems, a laser beam is scanned across the target surfaces automatically. The traveltime of reflected laser light pulses is used to calculate the distance to reflection points on the outcrop from the scanner. Each scanned point has horizontal coordinates (x and y) and elevation (z) information. All scanner setup points are tied together and to a global datum such as WGS84 (World Geodetic System 1984) using differential, real-time kinematic GPS. The rapid, automatic acquisition of data points creates a highly detailed point cloud. The point cloud can be used to create digital models of the terrain.

In this study, a mid-range laser measurement system (LMS), the Riegl Z360i scanner, and a long-range laser profile measuring (LPM) system, the Riegl LPM-i800HA, provided by Real Earth Models, LLC, were used to sample the outcrop topography. Laser scanner accuracy depends on both angle and range accuracy, and resolution depends on the ability to discriminate between objects in a point cloud. The resolution is determined by the smallest angle increment between successive samples (user definable) and the laser beam divergence rate. The LMS-Z360i scanner has an accuracy of 6 mm over a distance of 200 m and can capture up to 12,000 data points per second, with a resolution of 5 mm. The LPM-i800HA scanner has an accuracy of 15 mm over a distance of 800 m and can capture up to 1000 data points per second with a resolution of 1 cm. These systems produced a finely resolved data set with more than 1000 points per square meter.

Coarse-scale (decimeter-resolution), low-density (less than 15 points per 10 m2) capture of the same outcrop was accomplished using a Measurement Devices Ltd. robotic laser rangefinding system with a typical accuracy of 10 cm at 600 m. This system has a very low acquisition rate of ∼250 points per second. The high- and low-resolution data sets will be compared.

For a better definition of the outcrop surface geometry, the cliff faces were scanned from different locations (different angles) with both types of scanners. The multiple low-density scans (Measurement Devices Ltd. scanner) were aligned to each other and registered in the same coordinate system using the method described by Olariu et al. (2005). The Polyworks™ software (Alba and Scaioni, 2007) was used on the high-resolution point clouds from the Riegl scanners. The Polyworks™ alignment technique is based on the geometrical shape in overlapping scan areas, and, therefore, it does not require the use of reflective targets in each scan to align the point clouds. Overlapping portions of adjacent high-resolution scans should contain at least three features in common. Once these features are extracted, a six-parameter transformation (three rotations and three translations) can be estimated, and all the points of a scan can be transformed into a common coordinate system. The best sets of nonredundant data points in the overlapping regions are identified in the aligned 3-D scans and merged into a unified data set.

Surface Classification Algorithm Development

The laser scanner produces a 3-D point cloud that samples the outcrop topography at high spatial resolution at many thousands of individual points. With airborne lidar data, the quasi-normal direction to the surface is known to be vertical, and 2-D processing techniques may be applied. In ground-based 3-D scanner data for both natural and manmade objects, there can be many highly variable facet-normal directions. An automatic and objective method is desirable for the classification of the point data into discrete surface facets with well-defined orientations. The topography normal to each facet can be processed using well-known 2-D techniques for smoothing, decimation, and modeling, which do not readily generalize to 3-D. In some cases, such as fracture surfaces, a further analysis of the facet orientations themselves is of interest.

Unsupervised classification methods recognize patterns without the requirement for training data and with a minimal amount of human interaction. Many pattern recognition techniques make use of some form of cluster analysis. There are hundreds of published clustering algorithms (Hartigan, 1975; Spath, 1980; Kaufman and Rousseeuw, 1990) based on a wide variety of geometric, statistical, and heuristic concepts. The surface classification algorithm demonstrated here makes use of two different clustering methods, applied in two (or three) passes over the data in an unsupervised classification scheme.

A partition of a set of points in n dimensions is a division into k subsets. If a model, such as the subset centroid n-vector, characterizes each subset, then the partition can be optimized with respect to some measure of the model misfit to the data points. Optimized partitions for different integers, k, can be compared to each other, leading to the notion of partitioning into a certain number of clusters. Data points can be tested for membership in any particular cluster.

The surface classification method employs two partitioning steps. The first step (pass 1) is simply a coarse-graining of the data into groups, which are compact in 3-D configuration space (see Fig. 1A). These groups should be roughly equal in size (Navg) and have at least a minimum (Nmin < Navg) and no more than a maximum (Nmax > Navg) number of points. The parameter Nmin is determined by the requirement to compute cluster orientations for pass 2; Nmax must be small enough to prevent loss of spatial resolution after pass 1. This partitioning method does not need to have a well-defined cluster model but must be very economical to compute in order to cope with large point cloud data sets (Ntot >> Navg). The number of pass 1 clusters is k1 = Ntot/Navg. The improved, many-pass leader (multileader) algorithm (Hartigan, 1975) has been found to be a satisfactory pass 1 partitioning method. The second partitioning step (pass 2) must have a well-defined cluster model, but it can be less efficient in terms of computer resources since it will operate on a data set reduced in size by a factor of 1/Navg. The second partition will be created in six-dimensional space taking into account both the 3-D location and the 3-D orientation of the k1 clusters produced in pass 1. The partition in pass 2 is optimized with respect to an objective measure of model fit to produce a good descriptive model using the fewest number of surfaces (k2), as shown in Figure 1B.

It is anticipated that the clusters produced in pass 1 will occupy roughly disk-shaped volumes in 3-D configuration space. Each cluster therefore defines a facet, and subsequent aggregation of the facets will define surfaces. The points should spread out in two directions orthogonal to the disk-normal vector and be relatively confined in the normal direction. The normal direction can be specified by a 3-D unit vector. The centroid (mean position) vector and 3 × 3 covariance matrix is computed for each cluster, and eigenvalues and eigenvectors are then found. The unit-normal vector will be the eigenvector associated with the smallest eigenvalue (eigenvariance). The other two eigenvalues should be roughly equal in magnitude. Since nine parameters must be estimated in this computation, Nmin should be some small multiple (maybe 3 or 4) of nine. For all of the clusters, the centroid vectors are normalized (standardized) to approximately unit magnitude (on average), and the unit orientation vectors are appended to produce k1 points in a six-dimensional space for use in the next phase of the classification.

The k-means algorithm is a partitioning method that produces a specified number, k, of clusters, where each cluster is represented by its centroid vector in n dimensions (Hartigan, 1975). An iterative procedure is used to move points from cluster to cluster until a scalar error measure (least sum squared error) is approximately minimized. There are a number of variations on the basic idea that substitute different norms for the Euclidean distance (and least squares) used in the original method (Kaufman and Rousseeuw, 1990). Note that if k is equal to N, the number of points, the error can be made zero, and each point is its own cluster. The scalar error measure can be modified by the introduction of a penalty for the use of too many parameters (increasing k). Several such measures are available (Sakamoto et al., 1986; Hurvich, 1997). Bozdogan and Sclove (1984) applied one of these criteria (Akaike's information criteria or AIC) to cluster analysis. Possibly the best measure defined to date is the minimum description length (MDL) statistic (Rissanen, 1978, 1986; Barron et al., 1998). Minimization of MDL produces a parsimonious model with the minimum number of parameters. Use of this objective criterion permits the surface classification algorithm to function in an unsupervised mode.

The specific steps in the surface classification algorithm are:

(1) For Ntot data points in three dimensions, define the cluster size measures Nmin < Navg < Nmax and k1 = int(Ntot/Navg).

(2) Pass 1: Partition into k1 clusters using the multileader algorithm. Split any clusters larger than Nmax by application of the multileader algorithm. Determine the cluster medoids (the central data point in each cluster). Combine any cluster smaller than Nmin with the cluster having the smallest intermedoid distance. Adjust k1 as necessary.

(3) For each of k1 clusters, compute the centroid vectors, covariance matrices, and their eigenvalues and eigenvectors. Determine the cluster normal direction (eigenvector associated with the smallest eigenvalue). Normalize the centroid vectors and combine with the unit-normal vectors in a k1 by 6 data matrix. A weight factor can be used to favor either the location or orientation elements of the data. Note that without the normalization, the location will be favored over the orientation due to the large relative size of the 3-D coordinates (|(x,y,z)| >> 1).

(4) Pass 2: For a range of possible values kmin < k2 < kmax apply the k-means algorithm to the new data matrix. Use the minimum description length (MDL) statistic to find an optimal value of k2. Assign each of the original Ntot points to the appropriate cluster identified in pass 2.

The particular cluster methods (multileader and k-means) were chosen here for their availability as well as their effectiveness. There is currently a lot of work being done on new cluster analysis methods for application to very large data sets in data mining (Agarwal and Procopiuc, 2002; Baeza-Yates et al., 2003; Guha et al., 2003). Substitution of some of these new methods could greatly expand the applicability of the surface classification algorithm in practice.

Surface Classification Algorithm Applied to Synthetic Data

The algorithm was first tested and verified on synthetic data: 1338 points on six unfolded faces of a cube with known orientations, each of which had ∼223 uniformly distributed points (Fig. 1). Four faces of the cube were vertical planes; two of them were oriented north-south, and the other two were oriented east-west. The other two faces were horizontal planes. Therefore, we expected the final partition to produce three clusters corresponding to the three known orientations.

The first partition produced uniform clusters of about the same size (Fig. 1A). The algorithm was tested for different cluster sizes, but it produced better results for smaller clusters. A larger number of clusters for each face (fewer points in each cluster) proved to work better, especially in describing the face geometry at the edges. The second partition took into consideration the orientation of the clusters produced in the first partition and grouped the clusters with similar orientation together (Fig. 1B). The optimum number of clusters in the final partition was determined by the L-shaped break value of MDL criterion (Fig. 1C). Three clusters, corresponding to the three orientations of the six faces of the test shape, were identified. Parallel faces were assigned to the same cluster since they had the same orientation. Clusters located at the edges between faces were assigned to the face for which the cluster had the most points. Each pair of faces was correctly identified, and also their orientations were correctly estimated. Equal-area stereographic projection of the orientation data (Fig. 1D), showing the distribution of poles to the planes defined in the pass 1 partition, indicates that there is a strong clustering about the mean for each orientation.

GEOLOGIC SETTING AND STRUCTURAL DEFORMATION OF THE JACKFORK GROUP

The Jackfork Group provides classic outcrop exposures of deep-water siliciclastics deposited in an elongate east-west–trending Pennsylvanian basin. The Jackfork approaches 2100 m in original stratigraphic thickness and represents deposition during a period of 5 m.y. of structural quiescence (Coleman, 2000). Fractures in the Pennsylvanian section of the Ouachita belt document both regional and local stress fields present during the Ouachita orogeny. Most regionally traceable fractures are subparallel to the stress direction (N20W–N25E) within the Ouachita fold belt, and outcrops on fold noses and near thrust faults display fractures that are normal to the stress direction, indicating a local structural influence (Whitaker and Engelder, 2006).

Big Rock Quarry is located in the southeastern part of the Ouachita Mountains along the north bank of the Arkansas River in North Little Rock, Arkansas (Fig. 2). The cliff faces of the quarry expose a 3-D view of the lower part of the upper Jackfork Group (Jordan et al., 1993). The exposure has variable orientation and is almost 1 km long. The upper Jackfork at Big Rock forms an ∼60-m-thick succession preserved in the vicinity of Big Rock syncline, which extends over 16 km along the southern margin of the Maumelle zone of the northeastern Ouachita Mountains. The syncline trends N65E at its eastern extent but curves to an approximately N115E trend to the west (Hale and Connelly, 2002, 2003).

The Jackfork channel complex at Big Rock Quarry consists of erosive-based channels that stack up in a nested pattern. The channel complex infill consists mostly of amalgamated sand-rich individual channels. These massive sandstones comprise a significant proportion of the total channel fill, and their overall thickness decreases in an easterly direction where more shale is present (Slatt and Stone, 2001).

The sandstones exposed in the western part of the outcrop, which was the focus of our study, are predominantly well-cemented, fractured, amalgamated, and thick to very thick (several meters), tabular to irregular-bedded, fine-grained, and massive quartz arenite turbidites. The brittle, sandstone beds are bounded by thin (centimeter), ductile shale layers.

FRACTURE MAPPING AT BIG ROCK QUARRY

In the quarry, steep, nearly vertical wall exposures of thick, continuous Jackfork sandstones allow a 3-D visualization of the exposed fracture surfaces to be made. Two distinct natural extension fracture systems were identified in outcrop: quartz-filled veins and joints. Joints are planar tensile opening-mode fractures with little or no displacement parallel to the fracture plane (Price, 1966). They are associated with stress release during erosion and uplift (Plumb, 1982; Suppe, 1985; Price, 1990; Engelder, 1993) and are opened by movement perpendicular to the fracture planes. Fractures propagate through the rock layer perpendicular to the local minimum stress; therefore, their strike will mark the trajectory of the maximum horizontal stress at a specific time.

3-D Analysis of Terrain Surfaces Applied to the Low-Resolution Outcrop Data

The orientations of the cliff faces exposed on the western part of the quarry (Fig. 3) were determined using the surface classification algorithm. The 3-D point cloud that maps the entire northwest wall of the quarry (length 250 m, height 50 m) is made up of ∼3000 points. When compared with the high-density data set (Fig. 4C), the same number of points would define a surface only 2 m2. Due to the low resolution of the terrain data, the fracture pattern could not be identified, but the algorithm worked well in describing the overall cliff-face orientations (Fig. 5). Equal-area stereographic projection of the orientation data (Fig. 5D), showing the distribution of poles to the planes defined in the pass 1 partition, indicates that the clusters are oriented NS or NE-SW, which reflects the trend of the outcrop as can be seen in the map view in Figure 3.

3-D Analysis of Terrain Surfaces Applied to the High-Resolution Outcrop Data

Due to computer memory limitations, application of the surface classification algorithm to the whole outcrop at once was not possible. All processing for this paper was performed on a relatively modest laptop computer with 1.86 GHz clock, 504 Mb of central memory, and 900 Mb of virtual memory. Several areas (Fig. 4D) along the quarry wall have been selected to illustrate the efficacy of the classification methodology. These smaller portions of the outcrop were processed separately, and the results were compared. We tried to assess whether or not there were any lateral or vertical variations of the fracture pattern, as well as to appreciate if there were any changes in fracture orientations at different scales.

Fracture orientation measurements indicate that three major sets persist across the outcrop, and all of them are nearly perpendicular to bedding. The three identified orientations are: northeast (N21E/77E), northwest (N43W/79W), and north-northwest (N20W/81W). The NNW- and NE-striking fractures are subparallel to the regional northward compression during the Ouachita orogeny. The other fracture set is thought to have an orientation controlled by the principal stress axes that were active during folding in the region. Not all orientations were detected at every location. The NNW-striking one has smoother, more planar surfaces compared to the others, which have more irregular surfaces.

DISCUSSION

Fracture Attitudes at Different Localities

A comparison of fracture surfaces at different locations along the quarry walls revealed little difference in orientation between the dominant fractures sets (Figs. 6, 7, 8, and 9). The spacing between the various locations ranged from tens to several hundreds of meters (Fig. 4D), but even with larger distances, the fracture pattern was consistent, which is important from a reservoir point of view (Narr, 1991). Since the behavior of fractures did not change with location, we use a part of the outcrop—area 1 (Fig. 6), where the fracture planes are well defined—as a detailed example to show how the algorithm works, and for all the other locations, only the trend for the final orientations is presented.

Area 1 is a region of the outcrop where the fracture planes are easy to visualize. The thickness of amalgamated units reaches 2 m, and individual sandstone beds are ∼50 cm thick. Distinct planar orientations are shown as color-coded point clouds, where each color represents a particular orientation (Fig. 6C). All three fracture orientations were observed at this section of the outcrop: NE- (red), NNW- (yellow), and NW-striking (green points) fractures. All of them are nearly perpendicular to bedding (black points). The bedding plane surfaces strike approximately N26E and dip ∼6° to the east (Fig. 6E). The vector mean of the NE fracture strikes N23E. The NNW fracture set strikes N16W, and the NW fractures strike N53W.

In addition to extracting the main fracture orientations, the surface classification algorithm also gives the distribution of the identified orientations, which is of interest for reservoir modeling and simulation (Lorenz et al., 2002). The distribution of poles to fracture planes on the equal-area stereographic projection (Fig. 6E) indicates that there is more variation in strike orientation within each cluster than there is in dip variation, and the bedding planes have the highest variation in strike. We used the circular standard error (95% confidence cone semi-angle about the mean in degrees) as a measure of the dispersion of orientation about the mean; lower dispersion values result in tighter clusters in orientation about the mean (Fisher, 1993).

The NNW-striking fractures are more consistent in orientation (a smaller standard error) when compared to the other fracture sets 01(Table 1). For the NNW- and NW-oriented fracture surfaces, the strike dispersion exceeds dip dispersion, and, therefore, they are more planar when compared to the ones that strike NE (Fig. 6).

The fracture orientations determined by the surface classification algorithm for all studied locations at this outcrop were consistent with those obtained from other studies of the Pennsylvanian Jackfork Group in the Ouachita Mountains. Outcrop fracture characterization of turbidite sandstones using detailed geologic mapping of fracture pattern (Whitaker and Engelder, 2006) and linear scan-lines techniques (Combs-Scott and Smart, 2002) provided similar orientations (NNW- and NNE-striking fractures, with a range from N20W to N25E) for fracture surfaces.

Fracture Attitudes and Bed Thickness

The acquisition of data for different bed thicknesses resulted in fracture sets that had different distributions of both size and spacing resulting from the varied thickness of the sandstone beds. Fractures in the thin-bedded sandstone layers are densely spaced and confined to distinct layers. The relatively close spacing of the fractures in thinner beds is consistent with the relatively small bed thickness (centimeter to decimeter), while the thicker sandstone beds show an increase in fracture spacing with increasing bed thickness. The fractured, amalgamated sandstones comprise more than one stratigraphic bed, and the fractures are not limited to one single layer. The boundaries of the fractured sandstone layers are thin beds (centimeters) of shale.

The higher density of the fracture planes and the smaller scales for the thin-bedded (centimeters to decimeters) sandstone layers in area 2 made it more difficult to identify specific orientations (Fig. 7). The dispersion values 01(Table 1) for this section of the outcrop are higher when compared with the other locations where the sandstone layers are thicker. It was especially difficult to extract the orientation of the bedding planes. Even so, all three fracture orientations were identified, and the trends are close to the ones extracted from the previous locations.

For the thicker sandstone beds (thickness of amalgamated units is ∼2 m) in area 3 (Fig. 8) and area 4 (Fig. 9), it was easier to identify the fracture planes and to extract their orientations because of the larger number of points in the cloud to define particular surfaces. There was less variation along strike and dip for the larger-scale features. The NNW-striking fractures are more consistent in orientation (a smaller standard error) when compared to the other fracture sets 01(Table 1). The NW-striking fracture set has the highest dispersion, as was noticed also for the previous areas.

Measurements of orientations of the fracture surfaces at various locations along the quarry walls showed little variation of fracture orientation, even between widely separated areas (Fig. 10). The NNW- and NE-oriented fracture surfaces are more planar when compared to the ones striking NW, which have a rougher surface, and, therefore, they have lower dispersion values 01(Table 1).

There is a wide spectrum of fracture sizes found in fracture systems, and, in consequence, there is a scale dependence of fracture intensity. Because of limited direct fracture observation in the subsurface, especially for large fractures, different techniques (Narr and Lerche, 1984; Narr, 1991) have been developed that use the relationship between bed thickness and fracture density observed in cores/boreholes to extrapolate fracture spacing and orientation from 1-D to 3-D in the surrounding rocks (Wu and Pollard, 2002). This way, microfractures measured in cores/boreholes may provide estimates of macrofracture orientations (Laubach, 1997; Ortega and Marrett, 2000; Ortega et al., 2006).

Our results validate the previous findings about fracture variations with bed thickness, and the orientation information obtained from the 3-D outcrop can be extrapolated in the subsurface.

CONCLUSIONS

The main goal of this paper was to develop an automated and objective method to extract outcrop fracture surface orientations from 3-D point locations acquired using centimeter-accurate, high-density laser scanner systems. Our surface classification algorithm was effective in extraction of different orientations of the fracture surfaces. The identified trends in orientations were analyzed with respect to their variation with location and scale across the outcrop. Orientations obtained in multiple stratigraphic units indicated that the outcrop contains three distinct bed-perpendicular fracture sets. The main fracture orientations observed at this outcrop are consistent with those obtained from other studies of the Pennsylvanian Jackfork Group in the Ouachita Mountains.

The method proved especially successful with this type of outcrop—the exposure at Big Rock Quarry has a blocky aspect due to preferential weathering on the fracture planes. The blocky pattern of the terrain surface of the outcrop reveals the geometry of fractures and is easily captured by laser scanners. Although surface roughness noise is somewhat greater than for the larger-scale surfaces, the smaller-scale fractures in centimeter-thick beds can still be captured. More detailed point clouds are required for a better definition of the small surface facets.

*Corresponding author: ferguson@utdallas.edu

Bradley Bankhead of Veritas Exploration Services originally suggested this study. Norsk Hydro ASA (John Thurmond and Ole Martinsen) and the University of Texas at Dallas Digital Geology Consortium supported this project.