## Abstract

Data from a large-scale outcrop analogue (Upper Cretaceous Blackhawk Formation, Wasatch Plateau, central Utah, USA) were used to construct three-dimensional, object-based reservoir models of low to moderate net-to-gross (NTG) ratios (11–32%). Two descriptive spatial statistical measures, lacunarity and Ripley’s *K* function, were used to characterize sandbody distribution patterns in the different models. Lacunarity is sensitive to sandbody abundance and NTG ratio, while Ripley’s *K* function identifies clustered, random and regular spacing of sandbodies. The object-based modelling algorithm reproduces sandbody dimensions and abundances, but patterns of sandbody distribution generated by river avulsion are poorly replicated because pseudo-well spacing provides only limited constraint on sandbody positions.

In common with previous studies, the connected sand fraction in the reservoir models increases with increasing NTG ratio and increasing range of sandbody orientations, but there is significant stochastic variation around both of these trends. In addition, low NTG reservoir models in which sandbodies exhibit strong clustering may also have a low connected sand fraction across the model volume because the sandbody clusters are widely spaced and, thus, tend to be isolated from each other. Consequently, connected sand fraction could be overestimated if avulsion-generated sandbody clusters are not identified and replicated in models of such reservoirs.

## Introduction

Sandbody connectivity is a key control on recovery from channelized fluvial reservoirs of low net-to-gross (NTG) ratio (*c*. <30%) (Jones *et al*. 1995; Larue & Hovadik 2006; Hovadik & Larue 2007). Factors that influence the connectivity of channelized fluvial sandbodies include: (1) the NTG ratio; (2) the width and thickness of the sandbodies; (3) plan-view geometry of the sandbodies, which is typically considered in terms of sinuosity and range of sandbody orientations; (4) the organization of sandbody stacking; and (5) the sandstone content of crevasse-splay and other non-channelized floodplain deposits (Allen 1978; Leeder 1978; Bridge & Leeder 1979; Allard & HERESIM Group 1993; Jones *et al*. 1995; North 1996; Ainsworth 2005; Larue & Hovadik 2006; Donselaar & Overeem 2008; Pranter & Sommer 2011; Pranter *et al*. 2014). Reservoir modelling studies have typically assumed a combination of deterministic stratigraphic layering (e.g. using sequence stratigraphic concepts for alluvial strata: Wright & Marriott 1993; Shanley & McCabe 1994) and stochastic placement of channelized sandbodies within stratigraphic layers (e.g. using object-based modelling algorithms conditioned to vertical sandstone-proportion curves: Clemetsen *et al*. 1989; Hirst *et al*. 1993; Georgsen *et al*. 1994; Holden *et al*. 1998; Deutsch & Tran 2002). In its simplest form, the assumption of randomly distributed channelized sandbodies allows a percolation threshold for sandbody connectivity to be identified in each stratigraphic layer based solely on its NTG ratio (King 1990; Larue & Friedmann 2005; Larue & Hovadik 2006; Hovadik & Larue 2007). This percolation threshold for sandbody connectivity in three dimensions occurs at a NTG ratio of 28% (King 1990).

More complex stratigraphic architectures have been recognized routinely in studies of Holocene–modern rivers and their floodplains (e.g. Smith *et al*. 1989; Farrell 2001; Stouthamer & Berendsen 2007), and in process-based numerical models of floodplain deposition (e.g. Allen 1978; Leeder 1978; Bridge & Leeder 1979; Mackey & Bridge 1995; Karssenberg *et al*. 2001; Jerolmack & Paola 2007; Lopez *et al*. 2008). The architectures in these latter studies result in large part from autogenic behaviours such as lateral channel migration, which constructs channel-belt sandbodies, and channel avulsion, which generates spatially distinct channel-belt sandbodies. Different avulsion styles may result in compensational stacking of poorly connected channelized sandbodies or sandbody stacking due to channel re-occupation (Mohrig *et al*. 2000; Slingerland & Smith 2004). Avulsion-related stratigraphic architectures have also been characterized in ancient fluvial successions (e.g. Kraus 1996; Kraus & Wells 1999; Mohrig *et al*. 2000; Jones & Hajek 2007; Straub *et al*. 2009; Flood & Hampson 2014), and recent work on large-scale outcrops of fluvial strata suggests that avulsion and similar autogenic behaviours may result in non-random stratigraphic architectures that are at least partly independent of allogenic controls on stratigraphy and palaeogeography (e.g. sea level, tectonics, climate) (Hajek *et al*. 2010; Wang *et al*. 2011; Flood & Hampson 2015).

In this paper, we use data from a play-scale outcrop analogue, the Upper Cretaceous Blackhawk Formation (Wasatch Plateau, central Utah, USA), to test whether industry-standard, object-based methods, constrained to outcrop-derived logs and statistics of object dimensions, can mimic avulsion-generated sandbody distributions in alluvial to coastal-plain reservoirs of low to moderate NTG ratio (11–32%). Such spatial distributions are not commonly considered in current reservoir modelling workflows. Previous work has established the stratigraphic framework and context of the outcrop analogue (Hampson *et al*. 2012), interpreted avulsion as the dominant control on stratigraphic architecture using internal sandbody architecture (Hampson *et al*. 2013) and character of fine-grained floodplain deposits (Flood & Hampson 2014), and quantitatively described avulsion-generated patterns of sandbody distribution using spatial statistical tools (Flood & Hampson 2015). In this paper, we use concepts and data from this previous work to construct and analyse a suite of three-dimensional (3D), object-based reservoir models for different stratigraphic intervals and palaeogeographical locations in the play-scale outcrop analogue. We present a new method for deriving plausible, approximate distributions of the orientations and true widths of channelized sandbodies from cliff-face data that describe only their apparent widths (after Lorenz *et al*. 1985), and illustrate the application of statistical tools to compare the spatial distribution of channelized sandbodies between outcrop analogues and reservoir models.

## Geological setting of outcrop analogue

Alluvial to coastal-plain strata of the Blackhawk Formation were deposited along the western margin of a wide intracratonic seaway, the Cretaceous North American Western Interior Seaway (Kauffman & Caldwell 1993) (inset map in Fig. 1), as part of a prograding siliciclastic wedge of 3.5–4.0 Myr duration (Krystinik & DeJarnett 1995). The strata comprise channelized fluvial sandbodies surrounded by shale-rich floodplain deposits, which are also interbedded with lagoonal shales and thick coals in the lower part of the formation (Marley *et al*. 1979; Flores *et al*. 1984; Dubiel *et al*. 2000; Hampson *et al*. 2012). Regional mapping of sequence stratigraphic relationships across large-scale exposures in the eastern Wasatch Plateau and contiguous Book Cliffs, over an area of approximately 7500 km^{2}, indicates that the Blackhawk alluvial–coastal plain lay behind a series of approximately north–south-trending, wave-dominated deltaic shorelines that built out progressively further east (Balsley 1980; Flores *et al*. 1984; Dubiel *et al*. 2000; Howell & Flint 2003; Hampson & Howell 2005; Hampson *et al*. 2011) (Fig. 1). Overall progradation of the siliciclastic wedge was driven by a reduction in tectonic subsidence rate (Taylor & Lovell 1995; Hampson 2010; Hampson *et al*. 2012). Sediment was supplied from the Sevier Orogen (inset map in Fig. 1), probably from the Canyon Range Culmination approximately 80 km west of the Wasatch Plateau outcrop belt (Fig. 1) (DeCelles & Coogan 2006). Deposition occurred under a seasonal and warm temperate to subtropical climate (Parker 1976) at a palaeolatitude of approximately 42°N (Kauffman & Caldwell 1993).

Across the eastern Wasatch Plateau as a whole, channelized fluvial sandbodies increase in size and proportion of the strata from the base (*c*. 10%) to the top (*c*. 30%) of the Blackhawk Formation owing to progressively decreasing tectonic subsidence rate and increasing distance from the coeval shoreline (Hampson *et al*. 2012), although there is localized variation in this overall trend (e.g. Rittersbacher *et al*. 2014*b*) (Fig. 2). The proportion of channelized fluvial sandbodies increases abruptly to >90% across the unconformable base of the overlying lower Castlegate Sandstone (Adams & Bhattacharya 2005). Locally, mapped channelized fluvial sandbodies trend SW–NE (fig. 13 in Hampson *et al*. 2012) to NW–SE (Hampson *et al*. 2013), and their internal architecture suggests deposition mainly as channel belts formed by the migration of single-thread rivers of low–moderate sinuosity and multiple-thread, wandering to braided rivers that exhibited little variation in fluvial style with stratigraphic position (Adams & Bhattacharya 2005; Hampson *et al*. 2013). The channel-belt sandbodies are wider and less sinuous than the river palaeochannels (Hampson *et al*. 2013). Sequence stratigraphic interpretations assign the Blackhawk Formation to a highstand systems tract that is truncated by a sequence boundary at the base of the lower Castlegate Sandstone (Taylor & Lovell 1995; Howell & Flint 2003) (Fig. 2). Flooding surfaces within the highstand systems tract are expressed as laterally extensive coal zones developed up to 50 km inland from the coeval shoreline in the lower part of the Blackhawk Formation (Flores *et al*. 1984; Dubiel *et al*. 2000; Hampson *et al*. 2012), but are absent from the middle–upper Blackhawk Formation in the eastern Wasatch Plateau outcrops (Fig. 2). Channelized fluvial sandbodies in the Blackhawk Formation display strikingly little vertical (stratigraphic) or lateral (palaeogeographical) organization, implying that their distribution is governed largely by autogenic behaviours such as avulsion (Hampson *et al*. 2012, 2013; Flood & Hampson 2015). Hampson *et al*. (2012) subdivided the Blackhawk Formation into four intervals, based on mapped and projected coal zones (‘lower’, ‘upper 1’, ‘upper 2’ and ‘upper 3’ Blackhawk Formation, see Fig. 2). Channelized sandbodies tend to be clustered in the ‘lower’ and ‘upper 1’ intervals of the Blackhawk Formation, most probably reflecting deposition of these strata downstream of avulsion nodes near the apices of deltaic distributary channel networks (Flood & Hampson 2015). In the ‘upper 2’ and ‘upper 3’ intervals, sandbodies exhibit random or weakly defined regular spacing, the latter probably resulting from the compensational stacking of channel belts (Flood & Hampson 2015).

## Data set and methodology

### Outcrop analogue data set

Along the eastern edge of the Wasatch Plateau, the Blackhawk Formation crops out in a nearly continuous, subvertical, SSW–NNE-orientated cliff face that exposes up to 300 m of alluvial to coastal-plain strata over a distance of approximately 100 km (Figs 1 & 2). The cliff face is orientated subparallel to coeval palaeo-shorelines (Fig. 1), and thus exposes a play-scale cross-section that is oblique to regional depositional strike. Sandbody dimensions and distributions have been measured over the outcrop belt using panels reconstructed from oblique aerial photographs (e.g. Figs 3 & 4). The simple photograph-based method used to construct the panels is described in Hampson *et al*. (2012). Comparison of the panels with LiDAR (Light Detection and Ranging) data collected over two small sections of the outcrop belt suggest that the photograph-based method is adequate to make measurements of sandbody position with reasonable accuracy (<10% error) in near-vertical, near-linear sections of the cliff face with light vegetation and scree cover (Hampson *et al*. 2012; Rittersbacher *et al*. 2014*a*). Large (>3 m thick, >60 m wide) channelized sandbodies are consistently identified by both methods. Sandbody dimensions measured using the photograph-based method are less accurate than those measured using LiDAR data, but the associated errors are smaller than discrepancies introduced by subsequent interpretation of parts of the cliff faces with light scree or vegetation cover (Hampson *et al*. 2012). The panels have been ground-truthed using a conventional measured section and sandbody mapping from Link Canyon (Fig. 4b), which allows safe access (Hampson *et al*. 2012, 2013).

Data from five ‘windows’ from three cliff-face panels are used in this study (Fig. 4: after Flood & Hampson 2015): (A1) ‘lower’ Blackhawk Formation interval in panel A; (C1) ‘lower’ Blackhawk Formation interval in panel C; (C2) ‘upper’ Blackhawk Formation interval 1 in panel C; (C3) ‘upper’ Blackhawk Formation interval 2 in panel C; and (F4) ‘upper’ Blackhawk Formation interval 3 in panel F. These five ‘windows’ were selected because they sample all stratigraphic intervals exposed in the outcrop belt (Fig. 4), are collectively representative of the range of proportions of channelized sandbodies and sandbody distributions measured in the outcrop belt (Table 1), and contain the most densely ground-truthed part of the outcrop belt (Link Canyon in panel C, Fig. 4b) (Hampson *et al*. 2012, fig. 6, 2013). Each ‘window’ is rectangular in area, which has allowed straightforward application of spatial statistical tools that characterize sandbody distributions (Flood & Hampson 2015).

### Sandbody geometries, dimensions and orientations

Flood & Hampson (2015) measured the maximum thicknesses and apparent widths of channelized sandbodies in six panels from well-exposed, near-vertical, near-linear sections of the cliff face (panels A–F of Hampson *et al*. 2012), including sandbodies within the five ‘windows’ used in this study (Fig. 4). The internal architectures of many channelized sandbodies are difficult to resolve in photographs because they exhibit little internal variation in lithology or grain size. However, architectural analysis of sandbodies that are accessible from the ground, and which are considered to be representative, implies that channel belts are the predominant sandbody type (Adams & Bhattacharya 2005; Hampson *et al*. 2013). Each belt consists of multiple, laterally stacked bars and palaeochannel segments (i.e. channel storeys) that cut down from the same approximate stratigraphic level to form a multilateral body (*sensu*Potter 1967). Minorities of the sandbodies studied in detail from the ground are single-storey bodies, consisting of single-bar macroform and adjacent palaeochannel fill, or multistorey bodies (*sensu*Potter 1967) that comprise vertically stacked channel belts. The inclusion of single-storey and multistorey sandbodies in our analysis, in addition to multilateral, channel-belt sandbodies, may introduce bias into our results. Specifically, clustering of channel storeys within channel belts, and clustering of channel belts within multistorey sandbodies are not recognized. The implications of these omissions from our results, and for the application of our results to analysis of subsurface data, are discussed later.

Plan-view geometries of the channelized sandbodies are poorly constrained by the cliff-face panels. Channel-belt sandbodies are invariably wider and less sinuous than the palaeochannels that deposited them (e.g. Donselaar & Overeem 2008). Architectural analysis of accessible sandbodies suggests that palaeochannels in the Blackhawk Formation had low–moderate sinuosity (Adams & Bhattacharya 2005; Hampson *et al*. 2013). We thus interpret channel-belt sandbodies to be nearly straight over along-axis distances of up to 10 times the sandbody width, such as those considered in our models, as a first-order approximation. Low values of sinuosity (<1.1), consistent with interpreted channel-belt geometries, are applied to all sandbodies in our models. However, single-storey sandbodies are likely to be more sinuous, while multistorey sandbodies may exhibit variable and complex geometries along their axes. Sandbody width and thickness will also vary along the axis of all types of channelized sandbody (e.g. Donselaar & Overeem 2008). Such variations in sandbody plan-view geometries are not explored in our simple modelling experiments, and the implications of these omissions are discussed later.

Channelized-sandbody orientations are also poorly constrained by the cliff-face panels, which are essentially 2D cross-sections orientated oblique to regional depositional strike. Ground-based mapping and architectural analysis of a small number of sandbodies in Link Canyon (Fig. 4b, e) suggest that they trend NW–SE (N150º) (Hampson *et al*. 2013), but the mean orientation of channelized sandbodies across the whole outcrop belt is probably west–east (N090º), assuming that they trend perpendicular to the coeval palaeo-shorelines (e.g. Fig. 1). These two inferred orientations are used in combination with distributions of observed apparent sandbody widths to develop plausible scenarios for distributions of true sandbody widths. We use a seven-step workflow to generate a single normal distribution of sandbody orientations in combination with a single triangular distribution of true sandbody widths that can account, as a first-order approximation, for the projection effects inherent to a distribution of apparent sandbody widths:

A modal value of sandbody orientation is inferred from available data. In this case, scenarios are developed for modal sandbody azimuths of N090º (all five ‘windows’; red data in Fig. 5) and N150º (‘windows’ C1, C2 and C3; blue data in Fig. 5).

The modal value of apparent sandbody width is selected from the frequency distributions measured at outcrop (black histograms in Fig. 5b, d, f, h, j).

The corresponding value of true width is calculated by accounting for the orientation of the cliff-face panel and inferred modal value of the sandbody orientation.

This value of true width is assumed to be a mode, and multiple distributions of apparent sandbody width in the cliff-face panel are calculated assuming a normal distribution of orientations, for a range of angular standard deviations, using the method of Lorenz

*et al*. (1985; their appendix A).The distribution of apparent widths that fits the long tail of the distribution of measured apparent sandbody widths is selected (red and blue curves in Fig. 5b, d, f, h, j).

The minimum value of measured apparent sandbody width is used to calculate a value of minimum true width, using the orientation of the cliff-face panel and the same inferred modal sandbody orientation as in step (3).

A triangular distribution of true sandbody widths is assumed, calibrated to the modal and mimimum values calculated in steps (3) and (6).

### Variations in NTG ratio

Channelized fluvial sandbodies consist principally of medium-grained, cross-bedded sandstones with subordinate fine-grained sandstones and mudclast conglomerates (Hampson *et al*. 2013; Flood & Hampson 2014). Floodplain and lagoonal deposits mainly comprise siltstones, mudstones and coals, but contain sheets and lenses of very-fine- to fine-grained sandstones that are interpreted as crevasse-splay deposits (Flood & Hampson 2014). Thus, there is a distinct grain-size contrast between the channelized fluvial sandbodies and the sandstones in floodplain and lagoonal deposits. In the models described below, we consider channelized fluvial sandbodies to be net reservoir rock, and floodplain and lagoonal deposits to be non-reservoir rock. NTG ratio therefore corresponds to the proportion of channelized sandbodies. This choice of net reservoir rock implies that crevasse-splay sandstones would not contribute to flow (cf. Pranter *et al*. 2014).

Horizontal variations in NTG ratio are shown for each ‘window’ in Figure 6, together with the position of photographic logs and the measured section at Link Canyon (Fig. 4a–c). ‘Windows’ A1 and C1 contain the most pronounced NTG variations. These ‘windows’ contain relatively high NTG regions that correspond to one or more clusters of channelized sandbodies at different stratigraphic levels, separated by low NTG regions of isolated channelized sandbodies (compare Fig. 4a, b with Fig. 6a, b). ‘Windows’ C2, C3 and F4 contain more subdued NTG variation (Fig. 6c–e ) for two reasons. First, these three ‘windows’ contain wider channelized sandbodies (black histograms in Fig. 5h, k, n) than ‘windows’ A1 and C1 (black histograms in Fig. 5b, e). Secondly, ‘windows’ C2, C3 and F4 are documented to contain randomly or regularly spaced sandbodies over most length scales, whereas ‘windows’ A1 and C1 contain pronounced clustering of sandbodies (Flood & Hampson 2015). Variation in NTG ratio occurs in areas of all ‘windows’ between photographic logs and the measured section (Fig. 6).

### Design and implementation of modelling experiments

An industry-standard object-based modelling algorithm was used to construct reservoir models that capture simple, first-order representations of the sandbody geometries and NTG ratios for the five studied outcrop ‘windows’ (Fig. 4a–c). All models are 3D and share the same areal dimensions of 6000 m (north–south) × 4500 m (east–west), corresponding to regional depositional strike and dip of the Blackhawk Formation palaeo-shorelines. This area is sufficient to contain the outline of an exposure ‘window’ (Fig. 4d–f), and corresponds to a typical moderately sized reservoir. All models used in the sensitivity tests, described below, were 100 m in height, whereas the thickness of models of specific exposure ‘windows’ varies according to the ‘window’ height (Table 1). The resulting model volumes are relatively small relative to the dimensions of the largest sandbodies, which affects measurement of connectivity, but this facilitates the use of conditioning data (e.g. pseudo-wells) directly from the exposure ‘windows’. The models use a regular orthogonal grid with cells of uniform dimensions: 15 m (north–south) × 15 m (east–west) × 0.5 m (thickness). These cells are less than half of the minimum width (50 m) and thickness (2.0 m) of the smallest channelized sandbodies, and are thus sufficient to capture them as continuous objects on the grid. The models contain up to 24 million cells.

Channelized sandbodies are represented by objects of nearly uniform plan-view shape that differ in their cross-sectional dimensions and azimuth (e.g. red and blue data in Fig. 5). Sandbody sinuosity is defined by triangular distributions of amplitudes (100–200–300 m) and wavelengths (1000–3000–5000 m) that result in low values of sinuosity (<1.1: Fig. 7). The resulting object shapes are appropriate for the channel-belt sandbodies that predominate in the outcrop analogue, but may poorly represent the minority of sandbodies that comprise single-channel storeys or vertically stacked channel belts. The sandbody objects have no internal heterogeneity, as a first-order approximation of the limited variation in lithology and grain size observed in the channelized sandbodies at outcrop (Hampson *et al*. 2013). NTG distributions are constrained in some models using pseudo-wells that correspond to photographic logs (e.g. Fig. 3) and the measured section at Link Canyon (Fig. 6). Pseudo-wells are thus located along the position of the outcrop belt, and tend to constrain sandbody positions in a north–south-trending strip through the centre of each model. Pseudo-well spacing varies from 200 to 1200 m (Figs 4d–f & 6), and is thus comparable to well spacing in a reservoir. However, as noted above, the pseudo-wells provide relatively sparse sampling of NTG variation (Fig. 6). All models are constructed using a standard implementation of the object-based modelling algorithm. Channelized-sandbody objects are drawn from the distributions of dimensions and orientations specific to a particular model (Table 2), and inserted randomly into the model volume from an initial seed. In models that are conditioned to pseudo-wells, sandbody objects that do not honour pseudo-well data are rejected. Ten stochastic realizations were generated of most models in each experiment described below. This number of realizations is sufficient to establish a range of outcomes (e.g. Seifert & Jensen 2000). More complete statistical distributions would require a larger number of stochastic realizations, but the time required for non-automated processing of stochastic-realization outputs for spatial statistical analysis was deemed prohibitive.

Three sensitivity tests were designed to investigate how parameters that are poorly constrained in object-based models (stochastic variation in sandbody insertion point) and that are uncertain in the outcrop analogue data (orientation of cliff-face panel relative to sandbody orientation, range of sandbody orientations) influence spatial statistical measures of sandbody distribution. The sensitivity tests all use the same sandbody-object dimensions and a fixed NTG ratio (11%), both derived from a single stratigraphic interval (‘lower’ Blackhawk Formation interval, characterized in ‘windows’ A1 and C1: Figs 4a, b & 5a, f), but do not use pseudo-well data for model conditioning (e.g. Fig. 7; Table 2).

In the first sensitivity test, 10 models were constructed in which the only varying parameter was the seed number, which sets the start of the random number sequence for the algorithms that distribute channelized-sandbody objects. Each model is thus a different stochastic realization generated from the same input data (Table 2). Sandbodies were assumed to trend west–east (N090). Cross-sections were extracted along the north–south centreline of each model (i.e. along depositional strike) for analysis of sandbody distributions.

One model realization from the first sensitivity test was used to investigate the sensitivity of spatial statistical measures of sandbody distribution to the orientation of the 2D cross-section, relative to sandbody orientation, in the second sensitivity test (Table 2). Ten cross-sections were extracted from the model, oriented at azimuths of 0°, 10°, 20°, 30°, 40°, 50°, 60°, 70°, 80° and 90° relative to the north–south centreline of the model (i.e. the cross-sections are orientated from depositional strike to depositional dip). All cross-sections pass through the centre of the model.

In the third sensitivity test, the range in sandbody orientations was increased in 20° increments while sandbody dimensions and geometry were held constant (Table 2). Using these parameters, four stochastic realizations of 10 models were constructed with sandbody-azimuth ranges of 0° (N090°), 20° (N080°–N100°), 40° (N070°–N110°), 60° (N060°–N120°), 80° (N050°–N130°), 100° (N040°–N140°), 120° (N030°–N150°), 140° (N020°–N160°), 160° (N010°–N170°) and 180° (N000°–N180°) (Table 2). Mean sandbody azimuth is west–east (N090°) in each model. Cross-sections through the north–south centreline of each model were extracted for analysis.

The sensitivity test results provide a benchmark for comparison with a series of models that were constructed using sandbody-object dimensions and pseudo-well data specific to each outcrop ‘window’ (Figs 5 & 6). Conditioning data and parameters used to define object dimensions, orientations and distributions in the various experiments are summarized in Table 2. One set of models was generated for ‘windows’ A1 and F4, with a single distribution of sandbody widths and orientations used for each ‘window’ (red data in Fig. 5b, c, n, o). Two sets of models were generated for ‘windows’ C1, C2 and C3, each set with a different distribution of sandbody widths and orientations derived from outcrop measurements of the ‘window’ (red and blue data in Fig. 5e, f, h–l). These various models are intended to test the ability of the object-based modelling algorithm and conditioning data to mimic the apparent sandbody dimensions and patterns of sandbody distribution observed in each ‘window’, and to assess their impact on connectivity. In this context, pseudo-wells were used to constrain the location of stochastically inserted sandbodies. Cross-sections were extracted from each model along a centreline with the same orientation as the relevant exposure ‘window’, for analysis of apparent sandbody dimensions and sandbody distributions.

### Measurement of sandbody distributions

A number of statistical techniques that allow characterization of the spatial distribution of points and objects have been applied in the biological, ecological and geological sciences. Several are inappropriate for our analysis because they rely on locational information (e.g. quadrat method: Greig-Smith 1952), only enable comparison of results for areas of similar size (e.g. nearest-neighbour distance method: Clark & Evans 1954) or require robust identification of a chronologically ordered series of depositional horizons such as palaeosols (e.g. compensation index: Straub *et al*. 2009). Herein, we use lacunarity and Ripley’s *K* function, which are described below, to characterize sandbody distributions in 2D cross-sections from our reservoir models (Fig. 8). The potential sensitivity of these values to the 3D variability inherent to channelized sandbody distribution and stacking is investigated in the sensitivity tests, which are then used to inform analysis of later results. Values of lacunarity and Ripley’s *K* function measured from the outcrop ‘windows’ by Flood & Hampson (2015) are also used as a reference for comparison with model results.

Lacunarity is a pixel-based method that was initially developed to calculate the fractal dimensions of gaps between solid objects, and describes patterns of spatial dispersion in one, two or three dimensions (Allain & Cloitre 1991). Lacunarity has been used successfully as a measure of spatial dispersion in landscape ecology (e.g. Plotnick *et al*. 1993) and in a handful of geological studies (Henebry & Kux 1995; Plotnick *et al*. 1996; Rankey 2002; Roy *et al*. 2010). In this study, lacunarity is calculated using a sliding-box algorithm that is applied to a binary image of a model cross-section in which channelized sandbodies in the ‘foreground’ (black in Fig. 8a) are distinguished from floodplain and lagoonal deposits in the ‘background’ (white in Fig. 8a). The algorithm uses boxes of different sizes to sample the binary image (Allain & Cloitre 1991; Plotnick *et al*. 1996). A box of a particular size is placed at the top left of the binary image, and the number of foreground pixels within the box are counted. The box then slides incrementally across the image, with the number of foreground pixels counted after each increment of slide, until the entire image has been scanned. Twelve grid-box sizes were applied to each model cross-section, and boxes range in size between 2 and 45% of the cross-section area. The maximum box size was chosen to be less than 50% of the cross-section area, to avoid point statistical errors (Karperien 1999–2013).

Lacunarity has been used to compare spatial heterogeneity at different length scales, by comparing the results for different box sizes (e.g. Roy et al. 2010; Plotnick *et al*. 1996). Since there is little variation in lacunarity for the different grid-box sizes used in this study, spatial heterogeneity is averaged across the length scales of all grid-box sizes and over all grid orientations to generate a single, mean value of lacunarity for each model cross-section (Rasband 1997–2014; Karperien 1999–2013) (e.g. on the vertical axis of Fig. 6c). Low values of lacunarity (minimum = 0) are obtained from spatial patterns with evenly distributed, translationally invariant gaps of similar size (e.g. regular patterns), whereas complex patterns that display translational variance and contain unevenly distributed gaps of heterogeneous size (e.g. clustered patterns) are characterized by high values of lacunarity (maximum = 1). Values of lacunarity are not dependent on interpretation of sandbody type or hierarchy as the method relies only on the assignment of pixels to ‘foreground’ or ‘background’ lithologies.

Ripley’s *K* function is a spatial point process method for detecting deviations from spatial homogeneity (Ripley 1976). Ripley’s *K* function has been widely used in geography and ecology, and has recently been applied to stratigraphic analysis of channelized fluvial sandbodies (Hajek *et al*. 2010). The method determines how point pattern distributions change over different length scales within a data set (Ripley 1977). Ripley’s *K* function, *K*(*h*), is obtained in a 2D plane, such as an outcrop cliff face or model cross-section, using circles of radius *h* with their centres at each point (e.g. Cressie 1993; Rosenberg & Anderson 2011). The average number of points inside of these circles is calculated and divided by the number of points per area to obtain *K*(*h*). The distribution and presence of points is then evaluated at different values of *h*. If the number of points found at a certain distance is equal to the number of points expected, taking into account the intensity of the point process, then the distribution pattern is random. If more points are found within a given distance than the number expected, then this indicates clustering. If fewer points are found, then points are distributed regularly. We use 99 Monte Carlo simulations of a completely spatially random point process to establish a probability distribution for the number of points expected for the studied range of *h* (Rosenberg & Anderson 2011). A variance-stabilized version of Ripley’s *K* function, Besag’s *L* function (Besag 1977), is used here so that the *K* function can be compared to its expected value and against a benchmark of 0 (Rosenberg & Anderson 2011) (e.g. Fig. 8b).

We use the centroids of channelized sandbodies in a model cross-section as the points of interest (cf. Hajek *et al*. 2010) (e.g. Fig. 8a). Model cross-sections are vertically exaggerated by ×55.8, as were cliff-face ‘windows’ in the analysis of Flood & Hampson (2015), which corresponds to the ratio of mean apparent sandbody width to mean maximum sandbody thickness over the outcrop belt (i.e. in the six cliff-face panels of Flood & Hampson 2015, fig. 3). This vertical exaggeration minimizes the effects of anisotropy in sandbody dimensions on the results of our analysis, such that the expected spacing of sandbody centroids displays a constant mean and constant variance in all directions. Length scale is expressed in multiples of mean apparent sandbody dimensions (labelled ‘×1’, ‘×2’, etc., on the horizontal axis of Fig. 8c), and vertical and horizontal spacings of sandbody centroids are scaled according to mean maximum sandbody thickness and mean apparent sandbody width, respectively. In order to avoid distortion by edge effects, the maximum distance between points that is considered in our application of the *L* function is 25% of the width or height of each model cross-section (Rosenberg & Anderson 2011). The identification of sandbody centroids is sensitive to interpretation of sandbody hierarchy, and results will be most robust if only sandbodies of a particular hierarchical level are included in the analysis. Most of the channelized sandbodies at outcrop are interpreted as channel belts. However, the inclusion of channel-storey sandbodies is likely to increase the number of sandbody centroids and decrease their spacing, while the inclusion of sandbodies comprising vertically amalgamated channel belts is likely to decrease the number and increase the spacing of sandbody centroids.

### Measurement of connectivity

In order to avoid results that are contingent on well placement, we use two forms of geobody (or sandbody) connectivity (e.g. Larue & Hovadik 2006): (1) the largest connected geobody; and (2) the fractional volume of sandstone that is connected between the northern and southern faces of the models (i.e. the connected sand fraction *sensu*King 1990). The largest connected geobody may not contribute to the connected sand fraction, where the geobody does not intersect either the northern or southern model faces. Conversely, the connected sand fraction may comprise more than one geobody, where such geobodies are spatially distinct and each intersects the two model faces. The use of a connected sand fraction enables our model results to be compared with the predictions of percolation theory (King 1990). In our implementation of either metric, sandstone grid cells must share a common face to be classified as connected.

## Results and analysis

### Sensitivity tests

#### Impact of stochastic variability in sandbody distribution

Cross-sections from the 10 models are characterized by a range of values of lacunarity, from 0.33 to 0.47 (Fig. 9). Ripley’s *K* function analysis indicates that sandbody centroids are distributed randomly in six cross-sections (shown as entirely grey bars, which lack clustered or regular spacing at any length scale, in Fig. 9), but are weakly clustered or regularly spaced over certain length scales in the remaining four model cross-sections (shown as the black portions of grey–black bars in Fig. 9). The NTG ratio in the models is held constant at 11%, but 15%−33% of the sand fraction in the models occurs within the largest connected geobody.

This sensitivity test demonstrates that a range of values of lacunarity and sandbody-centroid distribution patterns can be generated solely by stochastic variations for the same input parameters, in the absence of pseudo-well conditioning data. Lacunarity values for the sandbody distributions observed in ‘windows’ A1 and C1 occur within the relatively wide envelope of lacunarity values for stochastically generated model cross-sections (Fig. 9). Lacunarity values for the model cross-sections differ by up to 0.11 from those for the exposure ‘windows’ (Fig. 9). However, marked clustering of sandbody centroids, as observed in ‘windows’ A1 and C1, is only generated in two of the 10 model cross-sections (Fig. 9). Thus, as might be expected, it appears unlikely that patterns of sandbody distribution can be reproduced repeatedly by stochastic models in the absence of additional conditioning data, such as pseudo-wells.

#### Impact of cross-section orientation

The lacunarity value and pattern of sandbody distribution for sandbodies of true width is given by a cross-section with azimuth of 0° from the north–south centreline of the model (i.e. perpendicular to the sandbody orientation, and parallel to depositional strike: N000º) (Fig. 10). As the azimuth of the cross-section is increased, the apparent sandbody width increases and the number of sandbodies intersected by the cross-section decreases. For cross-sections with azimuths of 10°–50° from the north–south centreline (i.e. oblique to depositional strike; N010º–N050º), values of lacunarity are slightly and consistently increased (by up to 0.07: Fig. 10a), and the pattern of sandbody-centroid distribution observed along depositional strike is retained (in this case, randomly distributed centroids at all length scales: Fig. 10b). For cross-sections with azimuths of 60°–90° from the north–south centreline (i.e. oblique and parallel to depositional dip: N060º–N090º), values of lacunarity progressively decrease with increasing azimuth (by up to 1.6: Fig. 10a). Patterns of sandbody-centroid distribution become variable (e.g. regular spacing observed at some length scales for azimuth of N060º: Fig. 10b) or impossible to calculate using Ripley’s *K* function because too few sandbodies are intersected (for azimuths of N080º and N090º: Fig. 10b).

This sensitivity test demonstrates that lacunarity and Ripley’s *K* function are relatively insensitive to the orientation of the cross-section in which they are measured, relative to mean sandbody orientation, provided that the cross-section is orientated parallel or oblique to depositional strike. The cliff-face exposure ‘windows’ have azimuths of between N000º and N050º (Table 1), and are orientated at 40°–90° from the inferred mean orientations of sandbody axes (red and blue data in Fig. 5). Thus, the orientations of the ‘windows’ are likely to have relatively little influence on the measured values of lacunarity and sandbody-centroid distributions.

#### Impact of range of sandbody orientations

Lacunarity increases slightly as the range of sandbody orientations is increased, although there is much scatter around this trend (Fig. 11a). The trend is interpreted to reflect increases in the range of apparent sandbody widths, and in the number of intersections and amalgamations between sandbodies as the range of sandbody orientations is increased. Sandbody-centroid distributions across the four stochastic realizations of each sandbody azimuth range show little variation as the sandbody azimuth range is increased (Fig. 11b). The fraction of sand volume that occurs within the largest connected geobody increases as the range of sandbody orientations is increased, from 17 to 27% for models with 0° range of sandbody orientations to 32–92% for models with 180° range of sandbody orientations (Fig. 11c). This trend reflects the increased connectivity of sandbodies as their range of orientation increases, as noted in previous studies (e.g. Larue & Hovadik 2006). This sensitivity test demonstrates that lacunarity in 2D cross-sections depends to a degree on the range of 3D sandbody orientations, and implies that there may be a positive correlation between lacunarity and sandbody connectivity.

### Models of specific outcrop ‘windows’

On visual inspection, the model cross-sections appear qualitatively similar to the cliff-face exposure ‘windows’ in terms of apparent sandbody dimensions and occurrence of amalgamated sandbodies (Fig. 12). The mean sandbody dimensions and sandbody densities (i.e. the numbers of sandbodies per unit area) in the model cross-sections are also quantitatively similar to those measured in the exposure ‘windows’ (Fig. 13).

The model cross-sections are characterized by lacunarity values that bracket those of the outcrop ‘windows’ (i.e. grey–red and grey–blue bars representing the model cross-sections are positioned above and below the grey–black bars that represent the outcrop ‘windows’ on the vertical axes of plots in Fig. 14). Lacunarity values for the model cross-sections differ by up to 0.08 from those for the exposure ‘windows’ (Fig. 14). However, patterns of sandbody-centroid distribution observed in the ‘windows’ are only sporadically reproduced in model cross-sections (i.e. red and blue portions of the grey–red and grey–blue bars that represent clustering or regular spacing of sandbody centroids in the model cross-sections are generally a poor match for the corresponding black portions of grey–black bars for the outcrop ‘windows’ on the horizontal axes of plots in Fig. 14). Only 35–60% of the appropriate model cross-sections reproduce the clustering of sandbody centroids observed in ‘windows’ A1, C1, C2 and C3 at any length scale (i.e. have red and blue portions of the grey–red and grey–blue bars on the left-hand side of plots in Fig. 14a–d). Fewer model cross-sections (0–15%) reproduce regular centroid spacings, which is observed in ‘windows’ A1, C1, C3 and F4, at any length scale (i.e. have red and blue portions of the grey–red and grey–blue bars on the right-hand side of plots in Fig. 14a, b, d, e). There is also no apparent trend in the length scales or degree of clustering or regular spacing that occurs for the range of NTG ratios captured in the exposure-‘window’ models (i.e. by visual comparison of the five parts of Fig. 14, for ‘windows’ with NTG ratios of 11–32%). Clustering is more pronounced in some suites of models with the same NTG ratio (e.g. red data in Fig. 14b and blue data in Fig. 14c), but such clustering is not attributable to differences in the range of sandbody orientations (e.g. red and blue data were generated using the same range of sandbody orientations in Fig 14b, c: Table 2) or in the range of true sandbody widths (e.g. red data were generated using a smaller range of true widths than blue data in Fig. 14b, and blue data were generated using a larger range of true widths than red data in Fig. 14c: Table 2).

Values of connected sand fraction in the models generally increase as the NTG ratio becomes higher, although there is a large spread of values between the 10 realizations of each model (red and blue data in Fig. 15). These results are expected in low-NTG reservoirs, which lie close to the NTG threshold (30–35%) at which isolated sandbodies become well connected (King 1990; Larue & Hovadik 2006; Hovadik & Larue 2007). Connected sand fractions in the models are greater than those predicted by percolation theory (black and grey data in Fig. 15: after King 1990), which assumes uniform sandbody dimensions, stationary NTG ratio and random sandbody placement. These assumptions do not hold true in our models. In particular, the presence of large sandbodies relative to the model dimensions (e.g. maximum sandbody dimensions in Table 2) increases connectivity and the variance in connectivity because there is insufficient model volume to allow a random distribution of sandbodies to be generated (i.e. there is insufficient volume support *sensu*Larue & Hovadik 2006; Hovadik & Larue 2007).

## Discussion

In this section, we revisit the aim of the paper in light of our modelling results. We answer four questions about the extent to which avulsion-generated sandbody distributions, as measured in the Blackhawk Formation data set, are reproduced using object-based methods and conditioning data. We then relate these sandbody distributions and their model representations to sandbody connectivity, and consider the wider implications for reservoir characterization and modelling.

### How well do the models capture the outcrop data?

The similarity between mean sandbody dimensions in the model cross-sections and those in the exposure ‘windows’ (Fig. 13a, b) is to be expected given that the thicknesses and widths of sandbody objects were drawn from input distributions derived from the outcrop data. However, this similarity corroborates the approach used to generate plausible, first-order approximations of sandbody width and orientation distributions (red and blue data in Fig. 5). Discrepancies between sandbody dimensions in the model cross-sections and exposure ‘windows’ arise from simplifications in these distributions, which are used as model inputs. For example, the slight, but consistent, overestimation of mean sandbody thicknesses in the models (Fig. 13a) can be attributed to the use of a triangular distribution that introduces a greater proportion of thick sandbodies than the positively skewed thickness distributions measured at outcrop (Fig. 5a, d, g, j, m). The similarity between sandbody densities (i.e. the numbers of sandbodies per unit area) in the model cross-sections and those in the exposure ‘windows’ (Fig. 13c) is also expected, as sandbody density is a function of NTG ratio and sandbody dimensions, both of which are specified as model inputs. In summary, the object-based modelling algorithm reproduces sandbody dimensions and densities from the input data, as it is designed to do.

The spatial distributions of sandbodies in the model cross-sections are more poorly representative of those in the exposure ‘windows’ (Fig. 14). Values of lacunarity for stochastic model cross-sections constrained by pseudo-wells are closer to lacunarity values for the exposure ‘windows’ (difference ≤ 0.08: Fig. 14) than values of lacunarity for stochastic model cross-sections that are not conditioned on pseudo-wells (difference ≤ 0.11: Fig. 9). Thus, patterns of sandbody spatial dispersion are constrained to some extent by the use of pseudo-wells. The degree of mismatch between pseudo-well-constrained model cross-sections and the exposure ‘windows’ is comparable to that introduced by variations in cross-section orientation of up to 50° from depositional strike (difference ≤ 0.07 for cross-section orientations of N000°–N050° in Fig. 10a). Clustering and regular spacing of sandbodies are poorly captured by the models, even with pseudo-wells providing some constraint on sandbody positions, and regular spacing appears to be harder to reproduce than clustering (Fig. 14). The number and distribution of available pseudo-wells does not impose sufficient restrictions on sandbody positions to force the object-based modelling algorithm to generate either clustering or regular spacing of sandbodies. Random sandbody distributions tend to be generated instead.

Based on the results and analysis presented above, the variance between model cross-sections and with the exposure ‘windows’ could both be reduced in four ways.

More pseudo-wells, including those located away from the outcrop belt to guide sandbody orientation and plan-view geometry, could be used to provide greater constraint on the insertion points of sandbodies by the object-based modelling algorithm. As a conceptual upper limit, uniformly distributed wells with spacing equivalent to mean sandbody width would constrain the position of every sandbody such that patterns and locations of clustered or regularly spaced sandbodies could be consistently reproduced in all stochastic object-based model realizations, but at the expense of 494–1200 pseudo-wells for models of the various outcrop ‘windows’ (i.e. 18–44 wells per km

^{2}, equivalent to less than 20 acre spacing). Such densely spaced well data are unlikely to be available in practice, but an order-of-magnitude fewer wells (i.e. equivalent to 40–160 acre spacing) would significantly increase the proportion of models that reproduce clustered or regularly spaced sandbodies.The models described above assume that all channelized sandbodies in the outcrop data set represent channel belts, including channel-storey and multistorey channel-belt-complex sandbodies that are interpreted to be in a minority. Application of a hierarchical architectural scheme to distinguish these sandbody categories would enable a different distribution of sandbody dimensions and geometries to be applied to each category. Application of such schemes to 1D data (cores, well logs) is far from simple because there is uncertainty in the interpretation of different sandbody categories (e.g. Bridge & Tye 2000; North 1996). Notwithstanding this uncertainty, the use of a hierarchical architectural scheme would more tightly define sandbody-object dimensions and geometries. However, it would offer only limited improvement in prediction of sandbody distributions, where greater uncertainty resides in the application of the object-based modelling algorithm.

We can apply a different modelling algorithm(s) and approach, which incorporates geological knowledge of avulsion processes that control sandbody distributions in the studied outcrop data set. Such approaches are in their infancy, but use outputs from process-based forward stratigraphic models and process-mimicking forward geostatistical models that are conditioned to well and/or seismic data (e.g. Karssenberg

*et al*. 2001; Cojan*et al*. 2004; Pyrcz*et al*. 2009).Some of the variation in the measurements of the model cross-sections is likely to result from statistical (or ergodic) fluctuations arising from the limited spatial extent of the model, which is insufficient to sample a whole random distribution of sandbodies (Srivastava 1996). These fluctuations could be suppressed by applying a distribution transformation that ensures precise reproduction of input statistics (Pyrcz & Deutsch 2014, p. 338).

### Which parameters control lacunarity?

Sensitivity tests show that spatial heterogeneity, as measured by lacunarity, increases with an increasing range of sandbody orientations (Fig. 11a), which induces an increasing range of apparent sandbody widths and also increases amalgamation of sandbodies with different azimuths (Fig. 11b). These two effects cannot be distinguished from each other by our implementation of lacunarity as a single value averaged across multiple length scales, but can be differentiated using Ripley’s *K* function (as discussed below). Spatial heterogeneity also increases as the orientation of the cross-section in which it is measured changes from a dip orientation (i.e. subparallel to sandbody axes) to a strike orientation (i.e. perpendicular to sandbody axes) (Fig. 10a). This change in cross-section orientation is associated with a decrease in mean apparent sandbody width and an increase in the number of sandbodies cut by the cross-section.

In models based on the exposure ‘windows’, spatial heterogeneity appears to decrease as NTG ratio is increased (Fig. 16a) and as sandbody density (i.e. number of sandbodies per unit area) is decreased (Fig. 16b). However, the models are small in scale and thus lack the volume support to define a random distribution of sandbodies (cf. Larue & Hovadik 2006; Hovadik & Larue 2007). Similar trends are noted in the exposure ‘windows’ (Flood & Hampson 2015, fig. 9), which are also small in scale. The differences in spatial homogeneity due to varying NTG ratio and sandbody density would be minimized in models or ‘windows’ of infinitely large volume. These two scale-dependent effects mask some of the trends documented in the sensitivity tests, perhaps because the effects of multiple parameters are convolved in models based on the exposure ‘windows’. For example, the largest range of sandbody orientations is present in models of ‘window’ F4 (angular standard deviation of ±50°: Table 2), but these have the second lowest set of lacunarity values (Figs 14 & 16a). Sandbody clustering is common in many of the ‘window’-based models, but bears no clear relationship with lacunarity (Fig. 14).

### Which parameters control the spatial patterns captured by Ripley’s K function?

Clustered, regularly spaced and randomly distributed sandbody centroids, as indicated by Ripley’s *K* function, may be generated simply by stochastic variations in sandbody placement, although regular spacing is relatively rare (occurring in only 20% of model cross-sections: Fig. 9). Sensitivity tests suggest that these spatial patterns are independent of the orientation of the observed cross-section, provided that the cross-section is orientated with a greater strike than dip component (Fig. 10b). Increasing the range of sandbody orientations tends to accentuate sandbody clustering (i.e. the combined length of the black portions of the grey–black bars in Fig. 11b tends to increase as the range of sandbody orientations increases) because there is greater intersection and amalgamation of sandbodies with different azimuths.

Models based on the exposure ‘windows’ are dominated by randomly distributed and clustered sandbody centroids, although regularly spaced centroids are present in some cases (0–15% of model cross-sections in the different parts of Fig. 14). There is no apparent trend in the degree of clustering or regular spacing and NTG ratios (i.e. by visual comparison of the five parts of Fig. 14, for models with NTG ratios of 11–32%). We attribute these spatial patterns in the exposure-‘window’ models to three controls. First, as noted above, the spacing and distribution of pseudo-well data used to condition the models are not sufficient to tightly constrain sandbody distribution patterns, and the object-based modelling algorithm tends to generate randomly distributed sandbodies where such conditioning data are sparse. Secondly, sandbody clustering at small length scales (*c*. <×1 in Fig. 14) may be underestimated because we have not used a hierarchical architectural scheme. Since we have assumed that sandbodies in the model cross-sections represent channel belts, the stacking of channel storeys within channel-belt sandbodies and of channel belts within multistorey channel-belt-complex sandbodies is not accounted for in the distribution of sandbody centroids. The analysis of exposure ‘windows’ presented by Flood & Hampson (2015), which we use as a comparison for the model cross-sections, similarly lacks the effects of such small-scale clustering. Thirdly, the dominant control on sandbody distribution in the Blackhawk Formation outcrop belt is interpreted to be avulsion, and resulting spatial patterns of sandbody distribution are expected to be complex at the length scales of the exposure ‘windows’. The interpretation of a predominant avulsion control is based on the near-absence of regional palaeogeographical and stratigraphic trends (Hampson *et al*. 2012), the detailed stacking patterns of channel belts within multistorey channel-belt-complex sandbodies (Hampson *et al*. 2013), and the types and distribution of overbank depositional styles (Flood & Hampson 2014). Avulsion was governed by the internal dynamics of the Blackhawk river systems. Details of the resulting stratigraphic patterns are contingent on the history of local sedimentation rates and accumulation of palaeotopography, which are specific to a particular ‘window’. Such stratigraphic patterns, and resulting sandbody distributions, are likely to be inherently complex and spatially variable across different length scales (e.g. Smith *et al*. 1989; Mackey & Bridge 1995; Kraus & Wells 1999; Mohrig *et al*. 2000; Farrell 2001; Slingerland & Smith 2004; Jerolmack & Paola 2007; Stouthamer & Berendsen 2007; Straub *et al*. 2009; Hajek *et al*. 2010).

### Which parameters control sandbody connectivity?

Sensitivity tests demonstrate that geobody connectivity increases as the range of sandbody orientations is increased (Fig. 11c) because sandbodies orientated at different azimuths intersect each other more frequently than those with the same azimuth (e.g. Larue & Hovadik 2006). However, the tests also demonstrate that a considerable range in geobody connectivity arises from stochastic variations in the insertion points of sandbodies into the model volume, for low NTG ratios (11% in Fig. 11c).

Models based on each exposure ‘window’ exhibit a wide range in connected sand fraction between the northern and southern faces of the models due to stochastic variability (Fig. 15). However, the mean value of the connected sand fraction increases as the NTG ratio of the model volume is increased (Fig. 15), as found in many previous studies (e.g. Allen 1978; Leeder 1978; Bridge & Leeder 1979; King 1990; Larue & Hovadik 2006; Hovadik & Larue 2007; Pranter & Sommer 2011). The effect of a relatively high NTG ratio (21%) may be confounded with that of a wide range of sandbody orientations (angular standard deviation of ±50°: Table 2) in models of ‘window’ F4 (Fig. 15). Connected sand fractions in all of the models are greater than those predicted by percolation theory (black–grey data in Fig. 15: after King 1990) because the model dimensions are too small to provide sufficient volume support (*sensu*Larue & Hovadik 2006; Hovadik & Larue 2007). As a result, the large size and small number of sandbodies in the models violate the assumption of a stationary NTG ratio that underlies the predictions of percolation theory. Variations in NTG ratio in the pseudo-wells used to condition the models (Fig. 6) thus enhance the potential to generate local clusters of sandbodies that are marked by increased geobody connectivity.

The influence of variations in spatial heterogeneity, as measured by lacunarity, and sandbody distribution, as measured by Ripley’s *K* function, on connected sand fraction in the exposure-‘window’ models are summarized in Figure 17. Models dominated by clustered sandbodies that also exhibit relatively spatially heterogeneous patterns (lacunarity >0.35) tend to have low values of connected sand fraction. Clusters of sandbodies (i.e. connected geobodies) in such models may be widely spaced at low values of NTG ratio (e.g. Fig. 18a, b), such that they are susceptible to being isolated from each other. In contrast, connected sand fraction tends to be high in models dominated by regularly spacing of sandbodies and in relatively spatially homogenous models (lacunarity <0.30) (Fig. 17). Sandbodies in such models tend to be well connected and also widely distributed in the model volume, such that connected sand fraction is high even at low NTG ratios (e.g. Fig. 18c, d).

### Implications for reservoir characterization and modelling

The results presented in this paper demonstrate that avulsion-generated stratigraphic architectures are highly variable in low–moderate NTG fluvial strata, and do not simply correspond to random sandbody distributions (e.g. Fig. 14). Of particular significance for reservoir characterization and modelling are clustered sandbodies that are generated by avulsion: for example, in locations downstream of nodal avulsion points (e.g. Mackey & Bridge 1995). Sandbody connectivity within these clusters is high, but the clusters themselves tend to be preferentially disconnected from each other in low NTG reservoir models (e.g. Figs 17 & 18a, b). The locations of these clustered sandbodies are not predictable from their stratigraphic or palaeogeographical context, but instead arise, at least in part, from the localized crossing of internal thresholds within the fluvial depositional system (e.g. Mackey & Bridge 1995; Mohrig *et al*. 2000; Slingerland & Smith 2004; Jerolmack & Paola 2007; Stouthamer & Berendsen 2007). In the Blackhawk Formation, these clusters occur principally on the lower coastal plain, where they are likely to represent the products of deltaic distributary networks (Flood & Hampson 2015), but similar clusters are also documented in alluvial strata that lack a marine connection (Hajek *et al*. 2010).

Our results suggest that modelling studies of low–moderate NTG fluvial reservoirs should aim to include scenarios of avulsion-generated stratigraphic architectures if they are to assess the full range of potential sandbody connectivities that may exist within such reservoirs. This aim may be accomplished by using industry-standard stochastic modelling techniques that are conditioned to spatial statistical measures such as lacunarity and the Ripley *K* function for appropriate analogues, as attempted in this study using an object-based modelling algorithm. Such an approach is difficult because the resulting models cannot be directly conditioned on spatial statistical measurements. Instead, a high density of conditioning data (e.g. numerous closely spaced wells) is required to force spatial patterns into the models, and such dense data sets are unrealistic in many subsurface settings. We consider that an alternative approach based on process-based forward stratigraphic models and process-mimicking forward geostatistical models (e.g. Karssenberg *et al*. 2001; Cojan *et al*. 2004; Pyrcz *et al*. 2009) is more promising because such models can be formulated to explicitly include avulsion processes or their effects. The challenges in applying these forward models include validating their outputs against data-rich modern and ancient reservoir-analogues, and conditioning the models to subsurface well and/or seismic data.

The production character of low–moderate NTG fluvial reservoirs that contain avulsion-generated stratigraphic architectures is expected to be dominated by highly variable production rates and depletion profiles between closely spaced wells or groups of wells, which reflect: (1) the localized penetration of sandbody clusters of anomalously large connected volume; and (2) poor connectivity between sandbody clusters. Patterns of sandbody connectivity in such reservoirs may not readily correspond to interpreted sequence stratigraphic trends, as is the case for sandbody distributions in the Blackhawk Formation analogue. In the apparent absence of the simple quasi-deterministic trends implied in sequence stratigraphic models, it is likely that sandbody distributions and connectivities will be modelled using a probabilistic approach that uses a very wide range of input parameters. This approach may be justified by reference to a complex geological conceptual model that invokes multiple allogenic controls, rather than by reference to avulsion and similar autogenic behaviours. Many fluvial reservoirs fit this prognosis of complex production behaviour (e.g. Shanley 2004), such that improved characterization and modelling of avulsion-generated sandbody distributions is likely to have wide applications.

## Conclusion

Large-scale exposures of the Upper Cretaceous Blackhawk Formation in the Wasatch Plateau, central Utah, USA have been used to define reservoir-scale templates of stratigraphic architecture generated by river avulsion in alluvial to coastal-plain strata of low–moderate NTG ratio (11–32%). The apparent widths and thicknesses of sandbody populations in 2D cliff-face exposure ‘windows’ were used to reconstruct distributions of sandbody dimensions and orientations in three dimensions that are both parsimonious and consistent with regional palaeogeographical reconstructions. In combination with pseudo-wells extracted from the exposure ‘windows’, these distributions of sandbody dimensions and orientations were used to construct 3D, object-based, stochastic reservoir models using industry-standard software. Spatial patterns of sandbody distribution in the exposure ‘windows’ and corresponding model cross-sections are characterized using two descriptive spatial statistical measures: lacunarity and Ripley’s *K* function. Comparison of the exposure ‘windows’ and model cross-sections is used to assess the degree to which the object-based modelling algorithm and conditioning data enable avulsion-generated sandbody distributions in different stratigraphic intervals and palaeogeographic locations to be captured in the 3D models.

Distributions of sandbody abundance and apparent dimensions in the outcrop analogue are easily reproduced in the 3D reservoir models using the object-based modelling algorithm. Spatial patterns of sandbody distribution observed in the outcrop analogue are reproduced with only partial success in the 3D object-based reservoir models. Spatial heterogeneity in the distribution of sandstone and shale, as measured by lacunarity, is sensitive to NTG ratio and sandbody density (i.e. the number of sandbodies per unit area). Cross-sections of stochastic model realizations are characterized by lacunarity values that bracket those of the corresponding exposure ‘windows’. Ripley’s *K* function is an effective tool to identify clustered, random and regular spacing of sandbodies, although it requires interpretation of sandbody type for optimal implementation. Cross-sections of stochastic model realizations only sporadically reproduce patterns of clustered or regularly spaced sandbodies observed in the exposure ‘windows’ because pseudo-wells are too widely spaced to force the object-based modelling algorithm to consistently generate such patterns.

Connected sand fraction in the 3D reservoir models (i.e. the fraction of sandstone connected between opposite faces of the models) is controlled by four parameters, although there is also much stochastic variation for the studied range of NTG ratio. In line with previous studies, sandbody connectivity increases with: (1) increasing NTG ratio; (2) increasing range of sandbody orientations; and (3) increasing size of the largest sandbodies, which are large compared to the model dimensions. In addition: (4) pronounced sandbody clustering reduces the connected sand fraction in low NTG reservoir models in which the sandbody clusters are widely spaced, and thus poorly connected with each other in the model volume. Our findings indicate that modelling studies of low–moderate NTG fluvial reservoirs should include avulsion-generated sandbody distributions if they are to assess the full range of potential sandstone connectivities that may exist within such reservoirs.

## Acknowledgements and Funding

The authors thank Chevron Energy Technology Company for funding the collection of field data to support this work, and Schlumberger for providing access to the Petrel software package via an academic software donation. Lacunarity analysis was carried out using ImageJ and FracLac software, developed, respectively, by Wayne Rasband (Research Services Branch, National Institute of Mental Health, USA) and Audrey Karperien (Charles Sturt University, Australia). Analysis using Ripley’s *K* function was carried out using PASSaGE v2 software, developed by Michael Rosenberg and Corey Anderson (Arizona State University, USA). ImageJ, FracLac and PASSaGE v2 software are freely available, and we are extremely grateful to their authors. We are also indebted to David Larue, Michael Pyrcz and four anonymous reviewers for their constructive criticism and suggestions, which have greatly improved the paper, and to editor Sebastian Geiger.