Hydrocarbon recovery in clastic reservoirs depends essentially on how well we understand the precise architecture of sand bodies and intercalated shaly baffles and barriers. Various methods have been developed for enriching the fundamental data collection from outcrop analogs; these include Terrestrial Laser Scanning (ground-based lidar), digital photogrammetry, high-precision GPS survey, etc. The three-dimensional outcrop data sets collected using these methods are critical for understanding the link between seismic-scale and well-scale data in the subsurface.

This study illustrates a methodology for integrating three-dimensional outcrop data, interpreting that data, and integrating the resulting interpretations with data from traditional outcrop measurements. A reservoir model of a fluvial sequence from the Escanilla Formation in the Ainsa Basin of northern Spain was produced using this methodology. Three-dimensional aerial photographs and laser scanner outcrop capture techniques provide a robust and flexible data set that can spatially constrain the modeling of observed features. The three-dimensional outcrop reconstruction, coupled with sequence stratigraphy concepts, enables the morphology, size, and distribution of key architectural elements to be modeled in subsurface reservoirs. The reservoir model constructed from these data allows geologists and reservoir engineers to evaluate the critical differences between real and modeled heterogeneities and provides a mechanism for an improved understanding of modeling subsurface reservoirs.

Braided fluvial reservoirs form some of the world's giant oilfields and are found in many petroleum provinces, notably the North Sea, North and South America, Yemen, Libya, and Australia. Many of these fields have been reviewed by Dreyer et al. (1993) and Martin (1993). Most of these fields are difficult to appraise and develop, due to small-scale geological variability often encountered in fluvial braided systems. Martin (1993) asked, “How can we adequately represent braided fluvial reservoirs in computer models?” and partially answered this question by proposing the use of stochastic modeling techniques when the average well spacing is less than the scale of the modeled heterogeneity.

Fluvial depositional systems are known to be highly dependent on allocyclic phenomena, especially in foreland basins, where the eustatic control on sedimentation is generally absent. The resulting fluvial architecture displays significant variability that has to be understood and quantified before it can be represented in reservoir models. Consequently, outcrop studies are essential for recording this variability, always bearing in mind that no outcrop analogue can match a reservoir perfectly (Kjemperud et al., 2003). From this perspective, there have been several studies of fluvial braided channels, and depositional processes are now consequently more clearly understood. However, a significant challenge exists in capturing the full three-dimensional geometry of channels from outcrop analogues. In general, the best constraint on three-dimensional geometry is provided by outcrops with irregular, highly indented topography; however, in this kind of outcrop, it is usually problematic to record and interpret the spatial variability of architectural elements using traditional field methods. For example, Zeito (1965) used a photomosaic panorama of cliffs in Colorado to measure the extent of sandstone and shale bodies in various depositional settings. In the case of highly indented outcrops, however, this methodology lacks precision and is difficult to manage and apply in three dimensions. The recent emergence of high-resolution positioning equipment (differential GPS) and new digital survey technologies (e.g., satellite imagery, digital elevation models, and airborne and terrestrial lidar acquisition) now allows the analysis of outcrop analogues in three dimensions (Lebel et al., 2001; Pringle et al., 2001; McCaffrey et al., 2005). The combination of these technologies, together with conventional field work, allows both sedimentology and stratigraphy to be measured and interpreted rapidly and precisely.

The outcrop described in this paper is located in the Ainsa Basin (South-Central Pyrenees, Spain) and is part of the Escanilla Formation (Late Eocene). The Escanilla Formation is a particularly well-exposed example of braided channel deposits, demonstrating rapid vertical aggradation on a short geological time scale. The highly irregular topographic surface provides an ideal test area for the application of the technologies and methods described below.

The workflow developed uses three-dimensional aerial photographs to define the regional stratigraphic framework, followed by terrestrial lidar acquisition to focus on specific areas of interest and to quantify the evolution of architectural elements seen in outcrop. Subsequently, the combination of lidar point-cloud data and field measurements is processed to establish quantitative morphological variations, which can then be superposed onto the previous stratigraphic framework.

The Ainsa Basin (also referred to as the Buil Syncline in the southern center of the Basin) is located at the western oblique margin of the South-Central Unit (Muñoz, 1992) on the Gavarine thrust sheet (Séguret, 1970) (Fig. 1). The central and western part of the southern Pyrenean foreland basin is formed by several N-S trending folds observed both in the Ainsa Basin and westward along the Sierra Exteriores (Kjemperud et al., 2003). The internal folds in the Ainsa Basin are interpreted as representing growth folds (Dreyer et al., 1999). The Ainsa Basin is bounded to the east by the Mediano Anticline and to the west by the Boltaña Anticline. Poblet et al. (1998) suggested that the Mediano Anticline was still active during the deposition of the Escanilla Formation (latest Eocene). The Boltaña Anticline, situated in the western part of the Basin, is a regional scale asymmetric anticline located above the western oblique ramp of the Gavarnie thrust sheet (Holl and Anastasio, 1995). Muños et al. (1998) suggest that the Anticline is a fault-propagation fold above a blind thrust. The Anticline developed between 50–36 Ma (Mutti et al., 1988; Puigde-fàbregas et al., 1989; Anastasio and Holl, 2001); i.e., active folding was occurring during the deposition of the Escanilla Formation.

The Escanilla Formation was deposited between the late Lutetian-Bartonian and the late Priabonian time (Bentham and Burbank, 1996) between ca. 41–34 Ma (Gradstein and Ogg, 2004). The drainage area of the Escanilla Formation is mainly from the Pyrenean massif through large valleys (Vincent, 2001). The Formation is divided into two members, the Mondot and Olson members.

The Mondot member is a transitional unit between the underlying deltaic Sobrarbe Formation and the alluvial Olson member. The Olson member consists only of alluvial deposits, and is unconformably overlain by alluvial fan deposits of the Collegats Formation (Dreyer et al., 1993; citing Atkinson, 1983). During the deposition of the Escanilla Formation, the transport direction changed between the Mondot and the Olson members (Fig. 2). The Mondot member is dominated by paleocurrent directions from southeast to northwest, and the Olson Member is characterized by paleocurrents converging from the north, northwest, and northeast, toward the south. These changes reflect the tectonic activity during sedimentation in the Buil Syncline.

The maximum preserved thickness of the Escanilla Formation is ∼1,000 m within the Ainsa Basin. The Formation thins toward the flanks of the Buil Syncline (Bentham et al., 1992). It is probable that the Escanilla Formation was originally deposited on top of the Mediano and Boltaña anticlines but has since been eroded (Bentham and Burbank, 1996).

The Escanilla Formation was subdivided by Kjemperud et al. (2003) into three main units based on changes in alluvial geometry and architecture. These three units are further subdivided into seven unconformity bounded sequences (Fig. 2). This outcrop study focuses on Kjemperud et al. (2003) sequence 4, located in the Olson area and described as having large lateral and vertical variations.

Large-Scale Investigation—Stratigraphic Framework

Fluvial systems respond to a number of allo-cyclic controls on sedimentation, including eustasy, climate, source area tectonism, and basin subsidence. The relative significance of these controls varies between the source area and the shoreline, as suggested by Shanley and McCabe (1994). According to previous criteria and the Shanley and McCabe classification, the Escanilla Formation was formed in the upstream area of the fluvial system, which is beyond any control and has no direct link with marine sealevel variations. In this setting, tectonism can play a major role in influencing the type of channel system deposited (Cant, 1978; Miall, 1981). The impact of tectonics on the architecture of upstream fluvial systems is illustrated in studies of modern systems (Coleman, 1969; Alexander and Leeder, 1987) and outcrop analogs (Heward, 1978; Gloppen and Steel, 1981; Lawrence and Williams, 1987; Nichols, 1987; Jollet et al., 1990; Turner, 1992; Garcia-Gil, 1993). The upstream portion of the fluvial system is also affected by climate and bedrock geology (Schumm, 1977; Westcott, 1993). Such climatic transitions have previously been observed in the Escanilla Formation (Dreyer et al., 1993; citing Atkinson, 1983; Kjemperud et al., 2003). The response of fluvial systems to climatic fluctuations is complex (Miall, 1996), influencing sediment load through time, vegetation types, and the degree of rainfall (Knighton, 1984). In spite of this climatic change overprint, the Escanilla Formation appears to be dominated by tectonic cycle events. Subsidence, superimposed on a background of long-term climatic variation, can lead to the development of cyclicity in fluvial deposits (Catuneanu, 2003). As with marine-influenced fluvial settings (middle to downstream fluvial regions), this cyclicity consists of a succession of fluvial sequences separated by subaerial unconformities. Fluvial successions isolated from marine influence are therefore dependent on tectonic events and/or climatic changes. According to Catuneanu and Elango (2001), sediment accumulation takes place during stages of flexural subsidence, whereas the bounding surfaces are related to stages of iso-static uplift.

In the study area, the lower boundary of the Kjemperud et al. (2003) sequence 4 is put at the base of a conglomeratic body that is basin-wide. The sequence is deposited in an alluvial plain setting, mainly covered by overbank deposits and some aggrading fluvial channel deposits. Aerial photographs draped over a digital elevation model (DEM), combined with lateral tracing in the field, were used to map the organization of larger-scaled braided channels and to define a regional stratigraphic framework (Figs. 3 and 4). Derived timelines converge toward the anticlines, underlining the strong tectonic influence on Escanilla Formation deposition. From a stratigraphic point of view, these lateral amalgamated channels can be linked to the low accommodation system tract, sharply eroding underlying floodplain paleosols. These extended braided belts are capped by maximum water table rises (base level rises) indicated by carbonate lacustrine deposits or wide rhizoco-dium or microcodium developments (Ryer, 1981; Atkinson, 1986; Kosir, 2004).

Aerial photographs also allow vertical aggra-dations of braided channels to be mapped. These aggradations are interpreted to have developed during the high accommodation system tract. The accommodation space in this setting is high during deposition, gradually decreasing upwards. The vertical aggradation of braided channels forms characteristic gravel chimneys (Fig. 5) with low lateral migration during deposition. The significant accommodation space, together with a high sediment supply, results in the deposition of thick floodplain mudstones that block lateral movement of the river systems.

At this scale, the main heterogeneity encountered is created by floodplain relics sandwiched between individual channels. These relics can either be encountered in lateral amalgamated low accommodation channel belts or in vertical amalgamated high accommodation braided channels (Fig. 6). Locally, lateral amalgamated channels can present lateral accretion features with their associated heterogeneity, formed by silty to shaly channel plugs (Fig. 7).

Channel Related Heterogeneities

Description of Channel Fill

Individual channel-fill sequences are typically between 5 and 10 m in thickness (Fig. 8), and at their base are composed of gravel dominated facies that sharply erode underlying reddened floodplain paleosols. Clasts range from 2 to 15 cm in maximum diameter and are sub- to well-rounded. They are organized along trough cross-stratifications, consisting of alternating conglomeratic and coarse-grained sandstones. The main heterogeneity encountered in this infill stage is formed by interbar mudstones (Fig. 8), which were deposited in protected or slack-water areas, comparable to those described by Lynds and Hajek (2006) from the North Loup River in Nebraska.

These gravel units are overlain by coarse- to medium-grained sandstones with trough cross-bedding to planar cross-stratification, fining upward to fine-grained sandstone with ripple laminations, mostly observed in the lateral amalgamated packages. The sequence ends classically in vertical amalgamated packages with abandonment deposits composed of bioturbated silty shales. This mottled facies is the second type of heterogeneity typically encountered in this channel infill sequence.

Interpretation of Channel Fill

This channel infill sequence is the result of a major braided-stream deposition. The coarse members are associated with coarse gravel transverse bar or bank-attached bar deposits, formed during the earliest stages of channel history. They tend to isolate protected or slack-water areas favorable for fine-grained deposition, preserved during channel avulsion stages. The local lateral accretion packages observed within lateral amalgamated belts are interpreted as lateral bar migrations, when the system tends to be more sinuous. The overall fining-up trend indicates a phase of channel system aggradation (Bentham et al., 1993). The sand-dominated planar cross-sets and trough cross-stratifications preserved within individual channel infills are interpreted in terms of the migration of transitory dune bed forms (Singh and Bhardwaj, 1991). The final infill stages, mostly dominated by mottled fine- to very fine-grained shaly deposits express the abandonment of the channels, with the emergence of vegetation and infauna.

Description of Floodplain Deposits

Extensive sheets ranging from bioturbated siltstones to fine-grained sandstones, with decimeter- to meter-scale thickness, have developed laterally to vertical amalgamated braided channels. These pinch out laterally away from adjacent channels into the surrounding floodplain shales. Locally, massive to planar stratified, medium- to fine-grained sandstones, with decimeter- to meter-scale thickness, are directly connected to braided channels and are tens of meters to hundreds of meters wide. These massive sandstones exhibit prograding features perpendicular to the channel paleocurrent directions.

Interpretation of Floodplain Deposits

Laterally extensive bioturbated siltstones are interpreted as overbank deposits that developed during flooding events or early stages of channel avulsion (Smith et al., 1989). The lack of sedimentary structures and high bioturbation rates (numerous horizontal burrows observed) are the results of extensive post-depositional reworking by flora and fauna. Massive prograding sandstones, interpreted as crevasse-splay deposits, occurring in the outer bank of channel meanders, are generated by a rupture in channel margin deposits created by rapid variations in flow velocity (Allen, 1965; Kraus, 1987).

Terrestrial laser scanning (TLS), also known as ground-based Lidar, is an ideal method for imaging surface topography in maximum detail and at rapid acquisition speed. Lidar technology (light detection and ranging) is now commonly used in a wide range of industrial, engineering, and scientific disciplines. Many different types of laser scanning equipment are currently available to respond to the varying needs of specific market segments.

The primary purpose of our data acquisition was to collect large-scale, channel geometry attributes over an area of up to 1 km2. The main requirement was a long-range scanner, capable of measuring several hundred meters across deeply incised valleys. Additional requirements were high-quality, true-color scan data and a means of pointing the scanner through a range of angles, upwards and downwards, in steep-sided canyons. These colored point-clouds are extremely useful during interpretation.

The Riegl LMS-Z420i scanner was chosen (Fig. 9). This scanner has a nominal range of 1,000 m and acquisition speeds of up to 12,000 points per second. On the top of the scanner, a tightly calibrated, precision-mounted, digital SLR camera provides true-color information precisely mapped onto the lidar point data. In operation, the scanner can rotate 360° around its principal axis, with an 80° field of vision through the scan window (Fig. 9). To scan at steeper angles, the scanner is tilted on its tilt mount at the relevant rotation axis (labeled in Fig. 9).

For our purposes, millimeter resolution and precision are unnecessary, and this scanner's angular resolution of 1 cm at 200 m was more than adequate. An additional requirement was differential GPS to georeference all the lidar data to a real-world global coordinate system. This was essential for integrating other spatially referenced data, such as GPS-located paleocurrent measurements and sedimentary logs, with the laser-scan point-cloud.

In the Olson area, 19 separate scans were acquired from eight tripod stations. The final data in the combined lidar point-cloud are in excess of 58 million points. Once lidar acquisition was complete, post-processing of the data was needed prior to sedimento-logical interpretation. Processing principally involves (1) applying color data from the digital photos onto each point in the point-cloud; (2) combining all the scans from the tripod stations within a local coordinate system; and (3) applying the GPS coordinates of known points from the scans to map the local coordinate system to real-world global coordinates. RiScan Pro software, delivered by Riegl, not only drives the LMS-Z420i scanner, but also enables us to carry out the post-processing procedure (Animation 11).

A deterministic model was constructed by combining sedimentological field data acquired using traditional field methods with geological parameters extracted from the lidar data. The area covered by the model is 900 m2 and 115 m thick. Approximately 1500 lines, representing facies transitions, geobody borders, and sedimentological features, were picked from the lidar point-cloud (Animation 2). Twenty-four sedimentological logs and 150 paleocurrent measurements were integrated with the lidar interpretation, using GPS locations.

The combination of all these data, together with the highly indented character of the out-crop, allows the deterministic reconstruction of individual channel morphology, using interpolation between outcrop faces constrained with paleocurrent measurements. Thirty individual channels were modeled and filled with identi-fied facies associations, based on lidar and field observations. In the final model, each individual channel is inserted in depositional order and erodes underlying previously modeled deposits. The final model is composed of more then 4,200,000 cells of 7 × 7 × 0.4 m.

This deterministic model provides the basis for characterizing heterogeneity linked with channel morphology and organization—e.g., for floodplain relics in which the resultant morphology is dependent on evolution of the sinuosity and spatial arrangement of individual channels. Hetero geneities can be extracted from the model to quantify their shape and dimensions (Animation 3). Such deterministic models can be used to generate flow simulations or to compute static connectivity between sand-dominated geobodies.

First observations of the outcrop suggest that vertical stacked channels are isolated. However, in-depth observations using the three-dimensional model, combined with static connectivity results, show that this is not the case. On the contrary, the resulting model emphasizes that static connectivity exists between laterally amalgamated channels. This connectivity is created by gravel chimneys, acting as vertical connections between the extensive low accommodation braided belts.

In addition to the deterministic modeling made possible by lidar acquisition, a significant number of morphological measurements an also be extracted from lidar point-clouds. These measurements allow a database that is capable of integrating morphological measurements of intervening architectural elements to be constructed. Channel morphology, over-bank, crevasse-splay, and associated heterogeneities were measured directly within the point-cloud and combined with all associated paleocurrent directions.

This databasecanbeusedtoconstructclassical outcrop dimension plots (Gibling, 2006)—e.g., to compare width and thickness ratios between different architectural elements (Fig. 10). The resulting values can be used in reservoir models to estimate the lateral extent of sandy geobodies at distances away from well bores.

We can glean additional information from the rich detail and high level of spatial precision of the lidar-based data set—for example, the study of the three-dimensional coherence of geobody distribution as well as its associated heterogeneity. This coherence also allows the vertical morphological evolutions of both sand geobodies and low permeability heterogeneities to be observed.

To perform a vertical analysis of the evolution of an architectural element, we extracted vertical facies proportions from the model (Figs. 11 and 12). The curve on the left side of Figures 11 and 12 represents the percentage of the different facies present along each layer of the gridded model. The high accommodation (HAST) and low accommodation (LAST) system tract stratigraphy is superimposed with plots of morphological measurements from different geobodies.

Low Accommodation System Tract Morphological Evolution

Low accommodation system tracts are characterized by a moderate to high sinuosity of individual channels, stacked laterally to form extended braided belts. These deposits, despite a high net-sand ratio ranging from 50% to 70%, are associated with a maximum preservation of channel clay plugs and interbar mudstones, inversely correlated with gravel preservation rate. Individual channel morphology can either be related to a morphological heritage, induced by early low accommodation erosion, or to the fluvial system response to increasing accommodation space. The fluvial system may maintain the gradient by increasing channel sinuosity (Richards, 1996). This channel sinuosity, which subsequently induces lateral accretion packages, allows the preservation of abandonment deposits marked by silty to shaly channel plugs. The frequent avulsion occurring in braided channel belts, together with the decrease in flow velocity, leads to local channel stream abandonment and preservation of slack-water deposits such as interbar mudstones. At the base of low accommodation system tracts, the gravel deposit percentage within channels indicates a coarsening upward trend, which is comparable to that described by Catuneanu and Sweet (1999) and is attributed to the gradual progradation of the fluvial system into the basin.

High Accommodation System Tract Morphological Evolution

High accommodation system tracts are characterized by an upward thickening and narrowing evolution of channel morphology. Channel sinuosity increases upward, associated with a decrease of the preserved gravel proportion fill. This evolution indicates a decrease of available accommodation space, coupled with a decrease of depositional energy. Maximum preservation of gravel fill is observed in the lower part of high accommodation system tracts, decreasing upward. This effect is interpreted as being generated by the hydraulic jump in conjunction with maximum base level rise.

Data integration using three-dimensional aerial photographs and terrestrial laser scanning (ground-based lidar), combined with field observations and measurements, allows stratigraphic horizons to be delineated in three dimensions and enables more precise three-dimensional morphological information to be extracted from the outcrop. The three-dimensional rendering and manipulation of the derived digital outcrop model provide access to a range of viewing perspectives, and architectural elements can be matched quantitatively through hillsides and across valleys. This is particularly relevant for irregular, undulating, and indented outcrops. This three-dimensional spatial coherence provides counterintuitive results, such as the static connectivity created by vertical gravel chimneys.

The main advantages of using terrestrial laser scanning (TLS) in this study include rapid data acquisition, rich data sets (tens of millions of points), and accurate high-resolution data (typically <5 cm). Standard field methods and subsequent interpretations are greatly enhanced. The method outlined in this paper provides a complex characterization of observed heterogeneities and geobodies, in which morphological evolution can be linked to a stratigraphic framework. These data can help to construct deterministic models that accurately reflect the spatial relationships of stratigraphic surfaces and architectural elements to a degree hitherto impossible to achieve.

The resulting quantitative spatial evolution of architectural elements can be readily incorporated into reservoir models; this increases the accuracy of predictions concerning various scales of heterogeneity in the subsurface. These quantitative spatial evolutions can be used as conditioning trends for stochastic modeling of subsurface geobodies and associated heterogeneities. The introduction of such spatial evolution of architectural element morphology is therefore a significant step forward in reducing uncertainties linked with random subsurface stochastic modeling.

The authors would like to thank Phillip Clegg for differential GPS measurements and moral support; Sophie Prunier-Dufarge, Stephanie Larmagnat, and Geoffray Musial for their help with field measurements; Dave Jennette and an anonymous reviewer for their constructive comments and remarks; all those students who did the “donkey work” (carrying the heavy equipment, etc.).

Supplementary data