# AAPG Bulletin

- Copyright ©2016. The American Association of Petroleum Geologists. All rights reserved.

## ABSTRACT

Tidal heterolithic sandstones are commonly characterized by millimeter- to centimeter-scale intercalations of mudstone and sandstone. Consequently, their effective flow properties are poorly predicted by (1) data that do not sample a representative volume or (2) models that fail to capture the complex three-dimensional architecture of sandstone and mudstone layers. We present a modeling approach in which surfaces are used to represent all geologic heterogeneities that control the spatial distribution of reservoir rock properties (surface-based modeling). The workflow uses template surfaces to represent heterogeneities classified by geometry instead of length scale. The topology of the template surfaces is described mathematically by a small number of geometric input parameters, and models are constructed stochastically. The methodology has been applied to generate generic, three-dimensional minimodels (9 m^{3} [∼318 ft^{3}] volume) of cross-bedded heterolithic sandstones representing trough and tabular cross bedding with differing proportions of sandstone and mudstone, using conditioning data from two outcrop analogs from a tide-dominated deltaic deposit. The minimodels capture the cross-stratified architectures observed in outcrop and are suitable for flow simulation, allowing computation of effective permeability values for use in larger-scale models. We show that mudstone drapes in cross-bedded heterolithic sandstones significantly reduce effective permeability and also impart permeability anisotropy in the horizontal as well as vertical flow directions. The workflow can be used with subsurface data, supplemented by outcrop analog observations, to generate effective permeability values to be derived for use in larger-scale reservoir models. The methodology could be applied to the characterization and modeling of heterogeneities in other types of sandstone reservoirs.

## INTRODUCTION

Heterolithic sandstones are commonly generated by tidal processes in shallow marine environments, such as deltaic and estuarine depositional systems. In these tidally influenced environments, the main current direction varies depending on the relative strength of tidal currents over daily to twice-daily cyclical time periods and the interaction of tidal currents with waves and river currents (e.g., Dalrymple and Choi, 2007). Sand is transported as bedload by strong currents to form ripples and dunes during periods of rising (flood) and falling (ebb) tide, and mudstone drapes are deposited during intervening slack-water periods. Depending on the flow regime, the mudstone drapes are more or less continuous over the sandy bedforms (Reineck and Wunderlich, 1968; Reineck and Singh, 1980; Nio and Yang, 1991). This results in interstratified, millimeter- to centimeter-thick sandstone and mudstone layers that are deposited over multiple tidal cycles and form the fine-scale heterogeneities that are characteristic of heterolithic tidal sandstone reservoirs. The distribution of mudstones and sandstones is delimited by a hierarchy of stratigraphic surfaces including (in order of increasing length scale) (1) lamina boundaries and reactivation surfaces that record incremental migration of bedforms, (2) the erosional bases of beds and bedsets, (3) boundaries between facies and facies associations, and (4) sequence stratigraphic surfaces. These four levels of stratigraphic surfaces define the multiscale architecture and connectivity of mudstone and sandstone layers, which, in turn, exerts a key control on the flow of gas, oil, and water during field production (Weber, 1986; Jackson et al., 2003, 2005; Nordahl et al., 2005, 2006; Ringrose et al., 2005; Nordahl and Ringrose, 2008).

The presence of these multiscale heterogeneities in heterolithic tidal sandstone reservoirs ensures that the characterization of effective reservoir properties such as permeability, relative permeability, and capillary pressure is a recurring problem (e.g., Martinius et al., 2005). Effective reservoir properties are typically derived from subsurface well data such as wireline logs and well tests, combined with laboratory measurements on cores and core plugs. Laboratory-derived reservoir properties are measured at a length scale that is small (of the order of centimeters for a typical core plug) compared to the dimensions of grid cells in reservoir simulation models (of the order of tens to hundreds of meters in plan view and tens of centimeters to meters in the vertical direction). In the case of tidal heterolithic sandstones, lateral and vertical variations in the continuity and connectivity of sandstone and mudstone laminae (e.g., meters to tens of meters) are not sampled by either subsurface well data or laboratory measurements. However, effective reservoir properties in heterolithic units are highly dependent on the volume sampled (Norris and Lewis, 1991; Jackson et al., 2003, 2005; Nordahl and Ringrose, 2008). Consequently, effective reservoir properties derived solely from subsurface and laboratory data in such heterolithic units are not representative of reservoir behavior; instead, models are required that capture the continuity and connectivity of sandstone and mudstone laminae at the appropriate length scale.

Two different methodologies have been used to create such models, which both use stratigraphic surfaces to reproduce multiscale heterogeneities. The first approach mimics depositional processes by generating and translating bedforms with a particular geometry according to user-defined inputs such as current velocity and sediment accumulation rate through time (e.g., Rubin, 1987; Wen et al., 1998; Rubin and Carter, 2005). Cross-stratification is defined by the preserved remnants of the bedform-bounding surfaces, whereas lithologies are distributed according to the local current velocities during deposition. This process-based methodology has been used to generate highly realistic models of near-wellbore regions (with dimensions of the order of 0.3 × 0.3 × 2 m [∼1 × 1 × 6.6 ft]) (Nordahl et al., 2005; Ringrose et al., 2005). However, process-based methodologies suffer from two problems. First, the models cannot be conditioned directly to data available from outcrop or subsurface measurements. Second, the required input parameters describing ancient depositional properties, such as variations in current velocity and sediment availability, are highly uncertain and have to be selected so as to produce a model that matches the preserved rock architecture observed in core or outcrop; this is a complex and nonunique inversion problem that is difficult to solve.

The second approach uses geometric and lithologic data from the subsurface in conjunction with outcrop analogs to directly condition reservoir models (e.g., White and Barton, 1999; Willis and White, 2000; White et al., 2004; Jackson et al., 2009; Sech et al., 2009). Jackson et al. (2005) generated three-dimensional (3-D) models of rock samples (with dimensions of order 0.5 × 0.5 × 0.3 m [∼1.6 × 1.6 × 1 ft]) from heterolithic tidal sandstones observed at outcrop using serial two-dimensional (2-D) sectioning, scanning, and surface reconstruction techniques. Their methodology yields models that are directly conditioned to observed geologic data, but its application relies on selection of an appropriate analog (or analogs) for the reservoir facies to be characterized. Furthermore, such a method is time-consuming, difficult to replicate, and leads to the creation of deterministic models that do not capture uncertainty in sand body proportions, geometry, and connectivity.

In this study, a surface-based modeling workflow is presented, which is then used to produce stochastic models of heterolithic, cross-bedded tidal sandstones conditioned to outcrop or subsurface data. A cross-bedding template surface is used to define and populate a rock volume. The 3-D morphology of the template surface is defined by purely geometric input parameters that, in the case documented herein, were defined using measurements from an outcrop analog (the Eocene Dir Abu Lifa Member, Western Desert, Egypt; Bown and Kraus, 1988; Legler et al., 2013). The models incorporate three of the four hierarchical levels of heterogeneity for heterolithic tidal sandstone reservoirs described above: (1) lamina boundaries and reactivation surfaces, (2) erosional bases of beds and bedsets, and (3) boundaries between facies and facies associations. The paper has four objectives. First, we present the new surface-based modeling workflow. Second, we identify the geometric input parameters required for the modeling process and extract a range of values for these parameters from statistical analysis of the outcrop analog data set. Third, we describe two generic models that reproduce (1) trough cross bedding dominated by muddy toesets and with a relatively low sandstone content (89%) and (2) tabular cross bedding dominated by sandy foresets and with a higher sandstone content (94%). Finally, we use flow simulation to calculate the effective permeability of the models to demonstrate the effectiveness of the surface-based modeling workflow and its application to build models suitable for flow simulation. In a companion paper (Massart et al., 2016), the surface-based methodology has been used to create a set of minimodels to investigate the range of effective permeability in heterolithic cross-bedded tidal sandstone facies.

## METHODOLOGY

### Model Construction Methodology

The stratigraphic surfaces that define sedimentary structures within tidal sandstone reservoirs can be categorized by their 3-D geometries, irrespective of length scale: (1) planar surfaces (parallel bedding; erosional or conformable facies contacts), (2) concave-upward surfaces (sigmoidal bedding or cross-bedding structures; channelized erosional contacts), or (3) wavy surfaces (wavy-bedding, lenticular-bedding, and flaser-bedding structures; irregular erosional contacts). The surface-based methodology uses these scale-independent stratigraphic surface geometries by modeling rock volumes within which surfaces share a common geometric template. This methodology comprises the following four steps (Figure 1).

The volume of rock to be modeled is subdivided into elemental volumes delimited by a basal and a top surface. In each elemental volume, the heterogeneities are associated with stratigraphic surfaces that have the same 3-D geometry. The elemental volumes have uniform shapes, but their dimensions can be varied. The model volume is filled with elemental volumes until an appropriate 3-D density is reached, in an approach analogous to object-based modeling (e.g., Haldorsen and Damsleth, 1990). Rules of superposition and erosion are applied to the elemental volumes to mimic their chronostratigraphic ordering. For example, if the elemental volumes represent erosionally based sediment bodies, then each elemental volume is eroded by the basal surfaces of younger elemental volumes.

Each elemental volume contains only one type of stratigraphic surface, the geometry of which is defined by a template surface. The 3-D geometry of the template surface is defined mathematically. Each elemental volume is then filled with numerous stratigraphic surfaces derived from the single template surface, following rules introduced by the user to define, for example, the vertical and horizontal surface spacing. The vertical and lateral extent of the surfaces within each elemental volume is controlled by the vertical and lateral extent of the elemental volume.

Once every elemental volume has been filled with template surfaces, a facies code is assigned to the geologic domains defined by the surfaces or to the surfaces themselves. The facies codes constrain the modeling of fine-scale petrophysical properties such as porosity and permeability.

The surface-based model is then gridded for flow simulation. The grid is constructed around the stratigraphic surfaces, to retain the geometries defined by the surfaces and minimize the number of active grid cells required for flow simulation (Jackson et al., 2005, 2009, 2013, 2014; Sech et al., 2009). The resulting models are geometrically accurate and computationally efficient, although the complex grid architectures may introduce numerical artifacts in conventional reservoir simulators (described in more detail in Massart et al., 2016; see also Graham et al., 2015).

### Application to Heterolithic, Cross-Bedded Tidal Sandstones

The four-step methodology described above is applied herein to the modeling of heterolithic, cross-bedded tidal sandstones (Figure 1). Cross-bedded sandstones are common in a wide range of depositional environments, including those influenced by tides (e.g., Harms et al., 1982; Rubin, 1987; Ashley, 1990). Cross beds result from the migration of dunes (or megaripples, sensu Allen, 1968; or sand waves, sensu Allen, 1980) in response to a unidirectional current. Dunes develop straight crests (2-D dunes) under low current velocities, and sinuous or discontinuous crests (3-D dunes) under higher current velocities (Dalrymple et al., 1978; Allen, 1980; Elliott and Gardiner, 1981; Middleton and Southard, 1984). Any dip section (parallel to the main current direction) gives the same geometry for tabular (or planar) cross beds resulting from the migration of 2-D dunes, whereas trough cross beds resulting from the migration of 3-D dunes have a more variable dip section geometry. Each migrating dune is preserved as a cross-bed set with an erosional base, whose geometry and extent reflect the morphology and trajectory of the scoured area in front of the migrating dune. In the case of 2-D dunes, the unidirectional current is dispersed along a large area downstream of the dune crest, such that an extensive planar erosion surface of low scour capacity is formed (Harms et al., 1982). In the case of 3-D dunes, the current is focused downstream of the migrating dune into scour pits, which migrate to produce a curved, concave-upward erosion surface (Dalrymple et al., 1978; Harms et al., 1982). Cross beds produced by dune migration are commonly stacked into larger sediment bodies of characteristic internal architecture. For example, the deposits of larger bedforms, such as bars, accumulate via the accretion of cross beds that record the migration of smaller, superposed bedforms, such as dunes and ripples, across the bar surface. Tidal bars migrate laterally into adjacent channels because of changes in tidal flow patterns or interactions with other processes (e.g., variations in wave climate or fluvial discharge). Consequently, tidal bar deposits can be composed entirely of stacked cross-bed sets, corresponding to the preserved remnants of repeated dune migration (Allen, 1980; Dalrymple, 1984; Ashley, 1990).

### Modeling of Elemental Volumes

A volume of 9 m^{3} (∼318 ft^{3}) (3 × 3 × 1 m [∼9.8 × 9.8 × 3.3 ft]) of cross-bedded sandstone is considered in this study; in a companion paper (Massart et al., 2016), we demonstrate that this volume comfortably exceeds the minimum volume (the representative elementary volume) required to calculate representative values of effective permeability in these dune-scale cross-bedded heterolithic units. At this length scale, the elemental volumes comprise tabular and trough cross-bed sets, representing the preserved parts of 2-D and 3-D dunes in a tidal bar succession, respectively. In each cross-bed set, the key heterogeneities captured are mudstone drapes along foreset–toeset surfaces, and each set corresponds to an elemental volume. The model volume of 9 m^{3} (∼918 ft^{3}) here samples approximately 6 cross-bed sets and 600 foreset–toeset surfaces, based on outcrop analog data presented in a later section.

Cross-bed set boundaries correspond to the preserved remnants of the erosional surface developed downcurrent of migrating 2-D or 3-D dunes (Figure 2A). Observations of modern tidal dunes show that this erosional surface has a curved, elliptical shape in the strike direction (orthogonal to the main paleocurrent direction; Figure 2B). As the dunes migrate, the resulting erosional surface is a downstream-amalgamated composite of the elliptical strike sections that record the successive positions of the deepest part of the scour pool in front of the dune (Figure 2B). To mimic the 3-D geometry of this composite erosional surface, the corresponding elemental volumes have been modeled here as ellipsoids (Figures 3, 4). The model volume is thus subdivided into ellipsoidal elemental volumes that correspond to cross-bed sets, with tops that are truncated by the basal surfaces of overlying elemental volumes (Figure 1B). The elemental volumes are modeled stochastically using the input parameters summarized in Table 1. For each parameter, the modeling algorithm can use a single value or a distribution characterized by a mean value and a standard deviation.

### Modeling Template Surfaces within Elemental Volumes

Each ellipsoidal elemental volume representing a cross-bed set contains multiple foreset–toeset template surfaces of uniform geometry. The spacing of foresets and toesets, and their associated mudstone drapes, is typically rhythmic, reflecting a hierarchy of periodic cycles in tidal current velocity (e.g., Nio and Yang, 1991). The shortest tidal cycle is semidiurnal (∼12-hr period) and is characterized by the alternation of flood and ebb current stages, separated by slack-water periods when the current velocity is zero. During slack-water periods, mud particles and clay aggregates (flocs) are deposited to form mudstone drapes over sandy bedforms (Allen, 1981; Dalrymple et al., 2003). In an idealized semidiurnal tidal cycle, both the ebb and flood tides are recorded by deposition of a sand lamina on the lee face (foreset) of a dune (Visser, 1980). Slack-water periods are recorded by mudstone drapes that separate the foreset–toeset sandstone laminae representing the ebb tidal and flood tidal currents. The tide is typically asymmetric, such that the ebb tidal or flood tidal currents either are of unequal velocity or are physically separated around the bar form (Visser, 1980). The dominant tide is represented by thicker foreset–toeset sandstone laminae, and the subordinate tide is represented by either thinner laminae or erosion (reactivation) surfaces.

An idealized, fully preserved semidiurnal tidal cycle is thus represented by two sandstone laminae and two mudstone drapes (paired mudstone drapes or mud couplet; Visser, 1980) that constitute one tidal bundle (Boersma, 1969). Longer tidal cycles, which are commonly preserved as rhythmic variations in the thickness of sandstone laminae and mudstone drapes within cross-bed sets, are diurnal (∼24-hr period) and spring-neap (∼14-day period) cycles. Superposition of the different tidal cycles, combined with other sediment transport processes, leads to preservation of sandy foresets and muddy toesets. A vertical profile through dune toeset deposits typically exhibits rhythmic alternation of millimeter- to centimeter-thick, wavy-bedded mudstone and sandstone laminae (Reineck and Singh, 1967). The transition between the foreset and toeset of each lamina in a cross-bed set is marked by a gradual reduced downcurrent curvature. The resulting foreset–toeset geometry may be referred to as shovel shaped (van den Berg et al., 2007). In a dip section, the shape of the foreset part is therefore approximated by a parabolic curve, and that of the toeset part is approximated by a straight line:(1)where *x* is the dip direction coordinate and *z* is the vertical coordinate, relative to the junction point *O* between the flat toeset part and the concave-upward foreset part (which is defined to be the origin, *x* = 0, *z* = 0). The whole toeset–foreset surface is then rotated by an angle *α*, which corresponds to the dip angle of the toeset. Equation 1 becomes(2)Notice that both equations have the same derivative *z*′(*x* = 0) = tan *α* at the junction point *O*, so that the curve is continuous from the toeset part to the foreset part of the surface.

In a strike section with coordinate *y*, the foreset and toeset geometry reflects the erosional scour at the base of the migrating dune, so that the resulting cross section in the strike direction corresponds to trough or tabular cross beds. Successive foreset–toeset surfaces are parallel to each other and parallel to the erosional base of the cross-bed set (i.e., elemental volume). Consequently, equation 2 is generalized for any (*x, y*) direction:(3)where *B*(*x, y*) describes the 3-D ellipsoidal shape of the basal surface of the cross-bed set (i.e., elemental volume):(4)The term ∑_{1}^{i}*T*_{Ti} corresponds to the cumulative toeset thickness, with *T*_{T} corresponding to the individual toeset thickness. Every cross section of one foreset–toeset surface in the strike direction is an ellipse parallel to the erosional base of the cross-bed set. In particular, for *x* = 0, the strike cross-section curve links all junction points *O* of any given foreset–toeset surface, creating a junction line *Oy*, simplifying equations 3 and 4 to yield(5)To populate the ellipsoidal elemental volumes with foreset–toeset template surfaces, the input parameters summarized in Table 1 are required (Figures 1C, 5). Toeset thicknesses *T*_{T} (Figure 5) are generally too small to be routinely measured directly from cores and outcrop analogs with high accuracy (<1 cm). Therefore, we calculate *T*_{T} indirectly from two other parameters: the dip angle of the toesets *α* and the angle of dune climb *δ* (Figure 5). The dip angle *α* corresponds to the angle of rotation applied to the parabolic curve representing the foreset–toeset template surface. We then give *T*_{T} as(6)

### Modeling of Mudstone Drapes along Foreset–Toeset Surfaces

If the succeeding flood tide or ebb tide is sufficiently strong, then mudstone drapes can be partially or entirely eroded, such that only one mudstone drape and one reactivation surface may be preserved during one flood-and-ebb tidal cycle (de Mowbray and Visser, 1984). Thus, the foreset–toeset surfaces modeled in the previous step may not be entirely covered by mudstone. The extent and continuity of mudstone drapes is defined using a function to describe the mudstone frequency in the dip direction along the stratigraphic surfaces, relative to a well-defined position on the surfaces. Mudstone drapes are modeled as elliptical patches of mudstone that are placed stochastically on each surface. Where they overlap, new patches erode older patches, so the patches coalesce to produce drapes with complex geometries. The length and aspect ratio of each elliptical patch is also modeled stochastically. Patches are placed on the stratigraphic surfaces until a user-specified proportion of their area is reached, following the methodology of Jackson and Muggeridge (2000). The mudstone frequency function denotes the probability that a patch will be placed at a certain location along each surface. The foreset part of each surface is modeled first; the mudstone drape coverage is then calculated at the transition between foreset and toeset parts (line *Δ* in Figure 5), and this calculated value is then used as the target mudstone drape coverage for the toeset parts, to ensure mudstone drape coverage continuity between the foreset and toeset parts. Consequently, toeset and foreset parts of each surface can have different mudstone drape coverage, allowing us to capture the muddy toesets typically observed in outcrop. The distribution of mudstone drapes along each surface is controlled by the chosen mudstone frequency function *f*, which is determined here from outcrop analog data. The following equation has been used to define *f*:(7)where *x*_{O} corresponds to the coordinate of the junction point *O* between foreset and toeset sections; *x*_{F} corresponds to the coordinate of the point *F* marking the preserved top of the foreset; and *M*, *N*, and *O* are constants that are chosen to fit data extracted from the outcrop analog. Such data could also be extracted from process-based models.

Mudstone drape thickness is user defined in the models and, at present, is assumed to be constant for each drape. Because mudstones are modeled as barriers to flow, their thickness has no effect on their flow properties; however, drape thickness does affect the total volume of the model that is occupied by mudstone. Here we have assumed a mud drape thickness of 3.5 mm (∼0.14 in.), which is a typical mean value encountered in heterolithic cross-bedded tidal sandstones (Terwindt, 1971; Nio and Yang, 1991; Martinius and van den Berg, 2011). Measurements of mudstone drape thickness could be taken from core data sets for application to a specific reservoir or from a suitable outcrop analog. The input parameters required for modeling mudstone drapes are summarized in Table 1 (Figure 1D).

### Outcrop Analog Data Analysis to Define Model Input Parameters

The input parameters required to construct the models of heterolithic, cross-bedded tidal sandstones were collected from an exceptionally well-exposed outcrop analog (see below), which enabled the 3-D geometry of the elemental volumes, template surfaces, and mudstone distribution to be evaluated quantitatively.

The studied outcrop analog forms part of the Eocene Dir Abu Lifa Member, located in the Western Desert of Egypt (Figure 6). The Dir Abu Lifa Member was deposited in a shallow-marine environment protected from wave energy, resulting in a predominance of tidal processes (Abdel-Fattah et al., 2010; Legler et al., 2013). The lower part of the Dir Abu Lifa Member consists largely of tidal bar and channel deposits that are stacked laterally and vertically (Legler et al., 2013). The lower parts of tidal bar deposits typically comprise heterolithic, cross-bedded sandstones.

The lower Dir Abu Lifa Member is exposed in a continuous escarpment over 20 km long, which is cut by multiple canyons that provide some 3-D control (e.g., Legler et al., 2013). The data sets used in this study are taken from two locations, labeled Gecko Nose and Butterfly Canyon in Figure 6. Gecko Nose is a small promontory that is defined by two cliff faces that trend approximately west–northwest-east–southeast and south–southwest-north–northeast, nearly perpendicular to each other (Figures 6, 7). The promontory exposes stacked trough and tabular cross-bedded sandstones, interpreted as the deposits of tidal bars in a channel belt (the yellow channel in the Gebel Sagha area of Legler et al., 2013). The west–northwest-east–southeast-oriented cliff face (N110–N290) is 17 m (∼56 ft) long, and the south–southwest-north–northeast-oriented cliff face (N030–N210) is 12 m (∼39 ft) long. Paleocurrent measurements from the cross beds are oriented toward N230, indicating that the west–northwest-east–southeast- and south–southwest-north–northeast-oriented cliff faces provide close to strike and dip sections, respectively. The slight deviation from the mean dip direction indicated by the paleocurrent data does not significantly affect the geometry of the modeled cross-bed sets.

Butterfly Canyon contains a larger promontory than Gecko Nose, defined by two cliff faces that trend approximately north–south and west–east. The Butterfly Canyon outcrop also exposes stacked trough and tabular cross-bedded sandstones deposited in bars occupying an isolated channel in a tidal flat environment (the Wadi Ghorab area of Legler et al., 2013). Paleocurrent measurements from the cross-beds are oriented toward N196, indicating that the west–east- and north–south-oriented cliff faces again provide close to strike and dip sections, respectively. Tidal bar deposits exposed at Butterfly Canyon are sandier than those at Gecko Nose, and the two deposits are considered to be end-members of the same heterolithic, cross-bedded tidal sandstone facies.

High-resolution photographs and precise sketches were collected from the cliff faces of both localities, to capture the dimensions and geometries of cross-bed sets. Photographs were collected using no-distortion lenses. Each cross-bed set in the Gecko Nose outcrop has been reconstructed from the high-resolution photographs and scaled using the sketches. The boundaries of the cross-bed sets and their constituent foreset–toeset surfaces have been traced on the reconstructed pictures, enabling quantitative, statistically representative data sets to be compiled for the various input parameters of the modeling methodology described above. All values are summarized in Table 1.

To define the dimensions of ellipsoidal elemental volumes (*L*_{E}, *W*_{E}, and *H*_{E}) we used data from the Gecko Nose outcrop. We determined *W*_{E} and *H*_{E} from 12 trough cross-bed set boundaries (identified in Figure 8A) from the strike-oriented face, using the method presented in Figure 4. The data set was limited to cross-bed sets with sufficient exposure to allow a best fit elliptical curve, with dimensions corresponding to *W*_{E} and *H*_{E}, to be fitted to their erosional basal surfaces (Figure 4C). We estimated *L*_{E} from cross-bed sets exposed on the dip-oriented face. The basal boundaries of all trough cross-bed sets in the dip-oriented face were continuous and nearly planar over the 12-m (∼39-ft) extent of the face, suggesting *L*_{E} >> *W*_{E}. No pinch-outs were observed. The elemental volume density *D* was determined from the total of 90 trough cross-bed sets observed at the Gecko Nose location within a volume of 12 × 17 × 3 m (∼39 × 56 × 10 ft), such that *D* is equal to 0.15 elemental volumes per cubic meter (4.16 × 10^{−3} per cubic feet). The dimensions of the preserved parts of the ellipsoidal elemental volumes (*L*_{A}, *W*_{A}, and *H*_{A}) were determined using data from both outcrops. We determined *W*_{A} and *H*_{A} from the strike-oriented face of Gecko Nose (Figure 8B). Values of both *W*_{A} and *H*_{A} define lognormal distributions (Figure 9). All of the cross-bed sets observed in the dip-oriented face of Gecko Nose are laterally continuous, in which case *L*_{A} > 12 m (∼39 ft). At Butterfly Canyon, *L*_{A} is observed in one cross-bed set to equal 25 m (∼82 ft), which is the value used thereafter.

To determine the degree of curvature *A*, the 90 foreset–toeset surfaces contained in three well-preserved trough cross-bed sets in the Gecko Nose outcrop (numbered 34, 50, and 52 in Figure 8A) have been extracted from photomontages. The three cross-bed sets show clear, dip-oriented cross sections of the foreset–toeset surfaces. The foreset–toeset surfaces are rotated in our analysis so that their toesets are horizontal. For each foreset–toeset surface, the junction point *O* is identified. All the foreset curves are then translated to the same origin, and a best-fit parabolic curve is fitted to the data (Figure 10). To determine the foreset thickness *F*_{T}, the sandstone laminae thicknesses between the 544 foreset–toeset surfaces contained in 12 studied cross-bed sets (identified in Figure 8A) have been measured after extraction of the surfaces from photomontages. A lognormal distribution of *F*_{T} values is observed (Figure 11). The dip angle of the toesets *α* has been measured on photopanoramas of the north–northeast-south–southwest-oriented (oblique dip-oriented) face of Gecko Nose. The angle of dune climb *δ* has been determined by generating a best-fit line through the foreset-to-toeset junction points *O* of laminae in each of the studied cross-bed sets.

To define the mudstone frequency function *f*, the positions of mudstone drapes along the same 90 foreset–toeset surfaces of the three cross-bed sets used to determine the parameter *A* (numbered 34, 50, and 52 in Figure 8A) have been extracted from photomontages. From this data set, a frequency distribution of mudstone drape presence relative to position along the foreset has been determined using equation 7 to define a best-fit curve (Figure 12).

## RESULTS

### Models Constructed from Outcrop Analog Data

The 3-D models of heterolithic, cross-bedded tidal sandstones are based on those observed at the Gecko Nose and the Butterfly Canyon localities. Generic models have been generated using input parameters derived from both localities (Figures 13, 14, Table 1). The models are stochastically generated using the data reported in the previous section, except for the elemental volumes, whose coordinates inside the model are extracted directly from photomontages of the two outcrop localities so that the traces of the elemental volumes in cross sections of the model accurately reproduced the cross-bed set boundaries of the outcrop sections. Both models are 9 m^{3} (∼318 ft^{3}) in volume (3 × 3 × 1 m [∼9.8 × 9.8 × 3.3 ft]) and contain four partially preserved ellipsoidal elemental volumes in the case of the Butterfly Canyon model and six partially preserved ellipsoidal elemental volumes in the case of the Gecko Nose model. The model volumes are approximately five orders of magnitude larger than the volume of a typical core plug (∼20 cm^{3} [∼1 in.^{3}]). Approximately 500 foreset–toeset surfaces are populated in in the Gecko Nose model, whereas only 170 of the same surfaces are present in the Butterfly Canyon model (parameter *N*_{CB} in Table 1). Note that the two outcrop localities both display examples of tabular and trough cross bedding. However, the Gecko Nose model shown here contains only trough cross beds, and the Butterfly Canyon model contains only tabular cross beds, reflecting the specific parts of the outcrops modeled. The mudstone drape coverage that was chosen for both models is equal to 25% along the foreset parts and 57% along the toeset parts of the foreset–toeset surfaces; the foreset drape coverage was extracted from outcrop data, whereas the toeset drape coverage is given by the fraction at the foreset–toeset junction resulting from the chosen mudstone frequency function (Figure 12), as outlined in the methodology.

A comparison between the outcrop cliff faces of Gecko Nose and Butterfly Canyon and the corresponding generic models is presented in Figure 15. The models honor the geometry of the cross-bed set boundaries, in both strike and dip directions, which validates the choice of having cross-bed sets represented as ellipsoidal elemental volumes (Figures 2–4). The input average foreset thickness *F*_{T} is respected as observed at the outcrop locations, with *F*_{T} being smaller for the Gecko Nose model than for the Butterfly Canyon model. The Gecko Nose model is relatively mudstone-rich (sandstone volume fraction *V*_{S}/*V*_{T} = 0.89) because it comprises trough cross-beds containing a relatively high proportion of toesets and thin foresets (foreset thickness *F*_{T} = 5.85 cm (∼2.30 in.) and foreset-to-toeset volume ratio *R*_{F/T} = 6.5:1; Table 1) caused by the high dune climb angle of *δ* = 5°. The Butterfly Canyon model is comparatively mudstone-poor (sandstone volume fraction *V*_{S}/*V*_{T} = 0.94) because it comprises tabular cross beds dominated by thicker foresets (foreset thickness *F*_{T} = 10.0 cm [∼3.94 in.] and foreset-to-toeset volume ratio *R*_{F/T} = 24:1; Table 1), with a low dune climb angle *δ* = 0°.

The distribution of mudstone drapes along the cross-bedding surfaces closely matches the distribution observed at outcrop. Both in the models and at outcrop, some mudstone drapes appear continuous along the entire cross-bedding surface, from the toeset part to the top of the foreset at the top boundary of the cross-bed set. In most cases, the mudstone drapes in the models are discontinuous over the entire length of the cross-bedding surfaces in dip-oriented cross sections, but the discontinuities are limited in the strike direction, which is again a close match to outcrop observations. Discontinuities of the mudstone drapes in the models are mostly located at the top of the foreset part of the cross-bedding surfaces, following the trend of the input mudstone drape frequency function defined from statistical analysis of outcrop data (Figure 12).

### Calculation of Effective Permeability

The method to calculate effective permeability is presented in the companion paper (Massart et al., 2016), and we report only the results here for the two models shown in Figures 13 and 14. We report the effective permeability as a normalized value, expressed as a fraction of the sandstone permeability:(8)The results reported in this way are independent of the value of sandstone permeability used in the models; moreover, the normalized effective permeability can be rescaled to any value of sandstone permeability obtained from core or minipermeameter measurements.

The effective permeability of the model volume has been extracted in three orthogonal directions: the horizontal effective permeability down depositional dip *k*_{d}, the horizontal effective permeability along depositional strike *k*_{s}, and the vertical effective permeability *k*_{v}. The mean value for horizontal permeability *k*_{h} is defined as the arithmetic average between *k*_{s} and *k*_{d} such as *k*_{h} = (*k*_{d} + *k*_{s})/2 (Jackson et al., 2003). The results for each model are summarized in Table 2. Despite the relatively low fraction of mudstone in the models, the presence of the mudstone drapes significantly reduces both horizontal and vertical permeabilities and introduces permeability anisotropy in the horizontal as well as the vertical directions : *k*_{s} ≠ *k*_{d} and *k*_{v} << *k*_{h}. For the trough cross-bedded model of Gecko Nose, the dip-oriented horizontal permeability (*k*_{d}^{n} = 47.5%) is only 65% of the strike-oriented horizontal permeability (*k*_{s}^{n} = 72.8%) because the flow must cross a larger number of mudstone-draped foresets when flowing down depositional dip as opposed to along depositional strike. The *k*_{v}/*k*_{h} ratio is reduced to only 0.040, reflecting that vertical flow must also cross numerous mudstone-draped foresets. The horizontal permeability anisotropy of the tabular cross-bedded Butterfly Canyon model is less pronounced than in the trough cross-bedded Gecko Nose model, reflecting the lower mudstone fraction and greater strike-oriented continuity of foreset–toeset sandstone laminae in the former model: the dip-oriented horizontal permeability (*k*_{d}^{n} = 70.0%) is 78% of the strike-oriented horizontal permeability (*k*_{s}^{n} = 90.0%). Despite the lower overall mudstone fraction, the *k*_{v}/*k*_{h} of the tabular cross-bedded Butterfly Canyon model is one order of magnitude smaller (at 0.003) than the value of the trough cross-bedded model. Mudstone drapes are approximately three times more numerous in the trough cross-bedded model (*N*_{CB} ≈ 500) than in the tabular cross-bedded model (*N*_{CB} ≈ 170); they are more densely spaced and laterally continuous in the toesets of the trough cross-bedded model because of the high values of the toeset dip angle *α*. Moreover, the sandstone volume fraction in the toesets of the tabular cross-bedded model is 0.26, which is less than half of the sandstone volume fraction in the toesets of the trough cross-bedded model, despite the common value of mudstone drape coverage of 57% in both models. The smaller value of *k*_{v}/*k*_{h} ratio in the tabular cross-bedded model arises from the closer spacing of mudstone-draped toesets, higher density of toeset surfaces, and consequent lower sandstone volume fraction in the toeset parts of cross-bed sets.

For the two models studied, the normalized effective permeability values can be rescaled to any measured sandstone permeability to yield estimates of effective permeability suitable for use in larger-scale reservoir models (Figure 16). For example, if the measured permeability of the sandstone was 500 md, then the trough-cross bedded (Gecko Nose) model yields permeability values of *k*_{d} = 238 md, *k*_{s} = 364 md, and *k*_{v} = 12 md, whereas the tabular cross-bedded (Butterfly Canyon) model yields permeability values of *k*_{d} = 350 md, *k*_{s} = 450 md, and *k*_{v} = 10 md. If the permeability of the sandstone was lower at 100 md, effective permeabilities are proportionately reduced as well, yielding *k*_{d} = 48 md, *k*_{s} = 73 md, and *k*_{v} = 3 md for the trough cross-bedded model and *k*_{d} = 70 md, *k*_{s} = 90 md, and *k*_{v} = 2 md for the tabular cross-bedded model. Effective permeability values from a broader range of model geometries and mudstone fractions are reported in the companion paper (Massart et al., 2016).

The modeling workflow reported here can be applied to create appropriate models for the calculation of effective permeability values depending on the geometric characteristics of the heterogeneity surfaces of any tidal cross-bedded heterolithic sandstone observed at outcrop location or in subsurface. The required input parameters and methods for measuring each of them are summarized in Figure 1. The main orientation of cross-bed sets (i.e., the paleocurrent) and its standard deviation can be deduced from dipmeter logs. The style of cross bedding is easily recognizable from the trace of the cross-bedding plane around the core or from borehole image logs. The tracing of the cross-bedding plane on well imagery can be considered for more precision because core observations are typically only possible on one-half of the core. The foreset thickness *F*_{T} can be measured on core from a representative number of occurrences, even if the typical width of a core (8–20 cm [∼3–8 in.]) prevents observation of a complete spring-neap cycle that displays cyclical variation of foreset thickness. The foreset-to-toeset ratio *R*_{F/T} can be appraised in a similar way from core observations but with similar limitations on the degree to which core data can represent variation in the parameter. For both parameters, an outcrop analog(s) can provide a more complete data set. Because dune climb angle *δ* has typically small values, the ratio *R*_{F/T} remains relatively uniform in the cross-bed set (Figures 3, 5). However, no lateral variation in *R*_{F/T} can be deduced from core observations. Finally, the toeset dip angle *α* can be observed in core if toeset areas are sampled. All input parameters can be otherwise derived from statistical analysis of appropriate outcrop analogs, in a similar way to the analysis presented herein using the Dir Abu Lifa Member as an outcrop analog, with the important proviso that the degree of analogy between subsurface and outcrop cases must be established with due care.

## CONCLUSIONS

This study presents a novel reservoir modeling methodology that accurately and efficiently reproduces the geometry and connectivity of sandstone and mudstone layers in heterolithic, cross-bedded tidal sandstones by stochastically modeling stratigraphic surfaces and associated heterogeneity. The model input parameters are geometric and can be derived from subsurface cores and/or outcrop analog observations. The application of the modeling methodology is demonstrated via the construction of models that represent heterogeneity in significantly larger volumes (9 m^{3} [∼318 ft^{3}]) than those sampled by core plugs (∼20 cm^{3} [∼1 in.^{3}]), using input parameters derived from analysis of an outcrop analog. Quantitative outcrop analog data are collated and used to constrain the geometry and spatial distribution of the small-scale heterogeneity surfaces (i.e., cross-bed set boundaries, cross-bedding foreset–toeset surfaces, and mudstone drapes). The resulting models are a close visual match to the outcrop data, such that the complex mudstone and sandstone connectivity of the heterolithic tidal deposits is accurately reproduced. The surface-based methodology is dependent not on length scale but on the geometric configuration and hierarchical arrangement of geologic surfaces. The methodology can, therefore, be applied to a much wider range of reservoir types in which the heterogeneity style can be characterized by the 3-D shape and distribution of geologic surfaces.

## ACKNOWLEDGMENTS

We gratefully acknowledge funding for this work and approval for publication from Norske Shell, Statoil, Total, and Petoro. We thank Roxar for providing academic licenses for RMS^{®} software, with which models were gridded and visualized, and Schlumberger for providing academic licenses for Petrel^{®} and ECLIPSE 100 software. The authors would like to thank Rodmar Ravnås and Marcus Sarginson (Norske Shell) for advice and guidance in the research project, as well as M.Sc. students Ayana Bréhéret and Justine Daboer (Imperial College London). We gratefully acknowledge C. Hern, M. Sweet, F. Whitehurst, and an anonymous reviewer for their constructive reviews and editorial comments.

## Footnotes

- Manuscript received December 8, 2014.
- Revised manuscript received October 12, 2015.
- Final acceptance July 10, 2015.
- Final acceptance February 1, 2016.

Benoît Y. G. Massart is part of the mature-area developments and increased oil recovery team at the Statoil research center in Bergen, Norway. He holds a Ph.D. in petroleum engineering from Imperial College London and an M.Sc. in petroleum geology from the École Nationale de Géologie de Nancy, France. His research interests are in reservoir modeling and quantifying the influence of geologic heterogeneity on effective rock properties and fluid flow behavior, facies modeling from geophysical data, and facies and petrophysical property update in history matching in ensemble methods.

Matthew D. Jackson is the Total Professor of Geological Fluid Mechanics in the Department of Earth Science and Engineering, Imperial College London. He holds a B.S. degree in physics from Imperial College and a Ph.D. in geological fluid mechanics from the University of Liverpool. His research interests include simulation of multiphase flow through porous media, representation of geologic heterogeneity in simulation models, and downhole monitoring and control in instrumented wells.

Gary J. Hampson is a reader of sedimentary geology in the Department of Earth Science and Engineering, Imperial College London. He holds a B.A. degree in natural sciences from the University of Cambridge and a Ph.D. in sedimentology and sequence stratigraphy from the University of Liverpool. His research interests lie in the understanding of depositional systems and their preserved stratigraphy and in applying this knowledge to reservoir characterization.

Howard D. Johnson holds the Shell Chair of Petroleum Geology in the Department of Earth Science and Engineering, Imperial College London. His main interests are in clastic sedimentology, sequence stratigraphy, reservoir characterization, and basin studies. Previously, he spent 15 years with Shell working in research, exploration and development geology, and petroleum engineering. He holds a B.Sc. degree in geology from the University of Liverpool and a D.Phil. in sedimentology from the University of Oxford.

Berit Legler is a clastic sedimentologist with Wintershall, Germany. She graduated in geology with a diplom degree and Ph.D. from the Technical University Mining Academy Freiberg, Germany. She worked as a sedimentologist with RWE Dea (2005–2008), postdoctoral researcher at Imperial College London (2009–2011), and lecturer at the University of Manchester (2011–2013). She currently works as sedimentologist in Wintershall’s reservoir geology technical services team.

Christopher A.-L. Jackson is currently the Statoil Professor of Basin Analysis in the Department of Earth Science and Engineering, Imperial College London. He obtained a B.Sc. and Ph.D. from the University of Manchester. His research interests are in the tectonostratigraphic evolution of rifts and the application of three-dimensional seismic data to understanding the formation and filling of sedimentary basins.