Abstract: 

Conceptual and quantitative models of fluvial stratigraphy typically argue that alluvial architecture is driven by a combination of allogenic processes (e.g., tectonic subsidence, sea level, sediment supply), and autogenic behaviors (e.g., avulsion). We test these models via analysis of a large (c. 200 m thick by 100 km depositional-strike extent) alluvial-to-coastal-plain succession that records shoreline progradation in response to decreasing tectonic subsidence (Late Cretaceous Blackhawk Formation, Wasatch Plateau, central Utah, U.S.A.). Oblique aerial photographs and architectural panels of six nearly linear, nearly vertical cliff face “windows” were used to characterize the dimensions (apparent widths, thicknesses) and spatial distributions of channelized fluvial sandbodies. Two statistical measures, raster-based lacunarity and point-based L-function, are used to investigate whether the distribution of channelized sandbodies show significant regularity, randomness, or clustering.

Over 490 channelized fluvial sandbodies are identified in the six cliff-face “windows.” From base to top of the Blackhawk Formation, sandbodies broadly increase in width and decrease in overall abundance per unit area. Clustering of sandbodies occurs relatively frequently in lower-coastal-plain strata (< 50 km from the coeval shoreline). Spatial regularity of sandbody spacing is consistently more apparent in upper-coastal-plain and alluvial strata (> 50 km from the coeval shoreline). A strong negative correlation between lacunarity and stratigraphic position is also observed, such that a wider range and greater variety of spatial gaps occurs between sandbodies in lower-coastal-plain strata relative to upper-coastal-plain and alluvial strata.

Comparison with numerical modeling studies of fluvial stratigraphic architecture implies the predominance of an avulsion-generated pattern of sandbody distribution that includes an element of compensational stacking on the upper coastal plain and alluvial plain, for a range of distances from the coeval shoreline (c. 40–130 km) and tectonic subsidence rates (c. 40–700 m/ Myr). On the lower coastal plain (< 50 km from the coeval shoreline) localized clustering of channelized sandbodies is interpreted to have formed by avulsion of deltaic distributary channels downstream of delta-apex avulsion nodes, possibly modulated by low-amplitude (< 30 m) sea-level fluctuations.

Introduction

Since Leeder (1978), Allen (1978), and Bridge and Leeder (1979) first recognized the relationship between avulsion, sedimentation rate, and the stacking density of fluvial channel-belt sandbodies, workers have attempted to utilize this information to develop improved models of alluvial stratigraphic architecture (Shanley and McCabe 1994; Hajek et al. 2010). In addition to strengthening theoretical understanding of how rivers operate in time and space, the analysis of spatial organization of channelized fluvial sandbodies is important because it enables better interpretation of parameters such as climate, tectonic subsidence, and sediment and water supply in the geologic record. Given that stratigraphic architecture impacts the connectivity of hydrocarbon reservoirs and groundwater aquifers, an understanding of the controlling mechanisms that determine sandbody organization is also required to predict reservoir and aquifer performance.

Conceptual frameworks derived from numerical modeling, modern analogues, and outcrop studies typically argue that alluvial architecture is controlled by allogenic mechanisms and/or by autogenic responses of the fluvial system (Allen 1978; Bridge and Leeder 1979; Wright and Marriott 1993; Shanley and McCabe 1994; Hajek and Wolinsky 2012). Allogenic forcing is due to extrabasinal processes such as climatic oscillation, changes in tectonic subsidence, or eustatic perturbations that occur over timescales of 105–107yr (Blum and Törnqvist 2000; Slingerland and Smith 2004; Muto et al. 2007; Stouthamer and Berendsen 2007; Hofmann et al. 2011). In comparison, autogenic behaviors that result from floodplain processes (e.g., channel avulsion and meander migration) operate over shorter timescales (103–105 yr), impose first-order controls on stratigraphy, and generate self-organization in the absence of, or independent from changes in allogenic forcing (Mackey and Bridge 1995; Blum and Törnqvist 2000; Muto et al. 2007; Straub et al. 2009; Hajek et al. 2010; Hofmann et al. 2011; Hajek and Wolinsky 2012). Both allogenic forcing and autogenic behavior can control avulsion frequency and sediment accumulation rate over a range of spatial and temporal scales (Allen 1978; Leeder 1978; Bridge and Leeder 1979; Mackey and Bridge 1995; Heller and Paola 1996; Mohrig et al. 2000; Hajek et al. 2010). An understanding of the relationship between allogenic forcing and autogenic behaviors is still poorly developed, because autogenic sediment transport processes can obscure or destroy allogenically derived environmental signals in alluvial architecture (“signal shredding” sensuJerolmack and Paola 2010). This effect creates difficulty in determining the temporal and spatial scales at which allogenic forcing and autogenic behaviors operated (Karssenberg and Bridge 2008; Wang et al. 2011).

Autogenic avulsion is considered to be a stochastic process that can generate spatial organization in stratigraphy, including two end-member patterns of sandbody spatial distribution: compensational stacking and clustered stacking (Straub et al. 2009; Hajek et al. 2010; Hofmann et al. 2011). Compensational stacking results from preferential infilling of previously unoccupied interchannel topographic lows on the floodplain by deposition (e.g., Straub et al. 2009; Wang et al. 2011). This spatial pattern is generated by channel belts that avoid the positions of the preceding channel belts and instead infill previously unoccupied positions on the floodplain (Straub et al. 2009; Wang et al. 2011). Clustered stacking patterns result from channel belts reoccupying a preferred position on the floodplain (e.g., Leeder 1978; Hajek et al. 2010). Channel-belt clustering can be generated by floodplain topography, with remnant channel scours and topography acting to redirect and attract newly avulsed channels towards pre-existing sites (Leeder 1978; Shanley and McCabe 1994; Jerolmack and Paola 2007; Hajek et al. 2010), or by spatial variations in sediment supply and/or the creation of accommodation space. Difficulty in interpreting these two end-member patterns occurs where deposits resulting from clustered stacking develop in topographically low positions on the floodplain, and can appear to be driven by compensational stacking (Straub et al. 2009).

In this paper, we use data from one large-scale continuous outcrop of alluvial-to-coastal-plain strata in which allogenic forcing is well constrained (Late Cretaceous Blackhawk Formation, Wasatch Plateau, central Utah, USA) to: (1) quantitatively characterize the geometry and spatial distribution of channelized fluvial sandbodies, and (2) determine the relative roles of allogenic forcing and autogenic behaviors in controlling stratigraphic architecture, via comparison with three widely cited numerical models of fluvial stratigraphic architecture (Bridge and Leeder 1979; Heller and Paola 1996; Karssenberg and Bridge 2008).

Geologic Context and Stratigraphic Framework

Deposition of the early to mid-Campanian Blackhawk Formation took place along the western margin of the epicontinental North American Western Interior Seaway (Fig. 1B) within a flexurally subsiding foredeep 180 km wide (Kauffman and Caldwell 1993; Krystinik and DeJarnett 1995; Liu et al. 2014). The Cretaceous seaway lay within a broad (1500 km wide), north–south-trending basin that formed by a combination of short-wavelength thrust-induced subsidence along its western margin and long-wavelength dynamic subsidence that operated across the entire seaway (Liu et al. 2011, 2014). During the Late Cretaceous, the Wasatch Plateau region was situated at a paleo-latitude of c. 45° N and experienced a subtropical to warm-temperate climate with seasonality, as indicated by the type and abundance of plants (Parker 1976; Kauffman and Caldwell 1993; Roberts and Kirschbaum 1995). The age and duration of the Blackhawk Formation is constrained by ammonite biozones calibrated against radiometric dates in correlative shallow-marine strata (Krystinik and DeJarnett 1995), and it represents a time interval of 3.5–4.0 Myr (Hampson 2010). Sandstone petrography indicates a probable sediment source in the inactive Canyon Range thrust-sheet culmination (Fig. 1A) in the Sevier Orogenic Belt, a fold-and-thrust belt situated c. 80 km west of the Wasatch Plateau outcrop belt (DeCelles and Coogan 2006).

The Blackhawk Formation constitutes part of the Mesaverde Group, which includes the overlying river-dominated Castlegate Sandstone and underlying shallow-marine Star Point Sandstone (Speiker and Reeside 1925; Dubiel et al. 2000). These three lithostratigraphic units together form an eastward-thinning siliciclastic wedge that passes basinward into offshore deposits of the Mancos Shale and records overall eastward progradation into the basin (Young 1995; Hampson 2010). Although coals are thick and abundant in the lower part of the Blackhawk Formation, most of the formation consists of channelized fluvial sandbodies and fine-grained, nonchannelized overbank deposits (Marley et al. 1979; Flores et al. 1984; Dubiel et al. 2000; Hampson et al. 2012; Flood and Hampson 2014).

From base to top of the Blackhawk Formation, channelized fluvial sandbodies are documented to become generally thicker, wider, and more abundant, and many have a multilateral, multistory internal character (Hampson et al. 2012). These trends have been attributed to a combination of three mechanisms: 1) decreasing tectonic subsidence rate (from c. 700 to 80 m/Myr) that forced overall progradation of the siliciclastic wedge, based on regional subsurface and outcrop mapping tied to dated ammonite biozones; 2) increasing distance from the coeval shoreline (0–100 km; Figs. 1, 2), based on regional mapping of paleoshoreline positions at outcrop; and 3) a potentially higher avulsion rate during deposition of the lower part of the formation (Hampson et al. 2012, 2013). Internal stratigraphic subdivisions of the Blackhawk Formation (“lower Blackhawk Formation,” “upper Blackhawk Formation 1,” “upper Blackhawk Formation 2,” and “upper Blackhawk Formation 3” in Fig. 2) are defined by a series of laterally continuous coal zones, which are associated with the up-dip pinchouts of shallow-marine parasequences in the underlying Star Point Sandstone (Flores et al. 1984; Dubiel et al. 2000; Hampson et al. 2012).

Based on the distribution of regionally extensive coal zones in the Blackhawk Formation, low-amplitude (< 30 m) sea-level fluctuations are interpreted to have controlled cyclic, high-resolution (< 0.1 Myr duration) stratigraphic patterns within lower coastal plain strata, situated < 50 km from the coeval shoreline (Dubiel et al. 2000; Hampson et al. 2013). Such low-amplitude sea-level variations are consistent with previous estimates of glacio-eustasy in the Late Cretaceous (Miller et al. 2003). Localized clustering of channelized fluvial sandbodies documented in upper-coastal-plain and alluvial strata (> 50 km from the paleoshoreline) cannot be readily linked to high-frequency allogenic controls because they do not occur at particular stratigraphic intervals or paleogeographic locations, so are instead interpreted to reflect autogenic behaviors driven by nonrandom channel avulsion (Hampson et al. 2012, 2013).

Dataset and Methodology

The dataset for this study was assembled from exposures situated in a large (c. 100 km wide, c. 250 m thick), NNE–SSW-oriented transect through the Wasatch Plateau outcrop belt, which is aligned subparallel to regional depositional strike (Flores et al. 1984; Dubiel et al. 2000; Hampson 2010) (Fig. 2). Due to the large scale of these exposures and their limited accessibility, the interpretations presented herein are based on the analysis of oblique aerial photographs, captured from nearly perpendicular positions to the cliff face using a light airplane as a camera platform (Hampson et al. 2012).

Channelized fluvial sandbodies were interpreted in 88 oblique aerial photographs, and the results were compiled to produce two-dimensional (2D) cliff-face panels for six “windows” in the Wasatch Plateau outcrop belt (labelled A–F in Figs. 1, 2). These six “windows” were selected because over 60% of sandstone exposure is visible, and the size of each window is sufficient (c. 10 km wide by 250 m thick; Table 1) to enable relatively large quantitative datasets to be extracted. Although channelized sandbodies are well exposed in the “windows,” surrounding fine-grained overbank and floodplain deposits tend to be covered by scree and vegetation. The internal architectures of the channelized fluvial sandbodies cannot be confidently interpreted in the dataset due to partial vegetation cover, variable photograph quality, parallax effects, and limited exposure of some sandbody margins. Therefore we do not categorize channelized fluvial sandbodies in our analysis.

Measurement of Sandbody Dimensions

The apparent width and thickness of each channelized fluvial sandbody was measured from the six digitized cliff-face panels (labelled A–F in Fig. 3). Measurements from the digitized 2D cliff-face panels and the photographs used to compile them are considered to be accurate to the nearest 20 m laterally and 2 m vertically, because the selected cliff faces are nearly vertical and nearly linear, and lack heavy vegetation or scree cover, which minimizes perspective and parallax effects (Hampson et al. 2012). This interpreted accuracy of measurement is supported by comparison with high-resolution (c. 0.1 m) spatial data acquired by an oblique-view, helicopter-borne light detection and ranging (LIDAR) system over two sections of the outcrop belt (Rittersbacher et al. 2014). Relative discrepancies in measurement between the photograph-based data and LIDAR data are greater for sandbody thicknesses and widths than for sandbody positions (fig. 10 in Rittersbacher et al. 2014), because the latter are based on large measurements (Hampson et al. 2012). Errors due to differences in interpretation of scree-covered parts of the cliff face, which are independent of the data-collection method, are larger than measurement discrepancies between the two dataset types (Hampson et al. 2012). Errors introduced during the measurement of sandbody dimensions from the digitized cliff-face panels are < 4 m horizontally and < 8 cm vertically, and reflect image resolution and measurement repeatability. The measurement errors outlined above may cause very small, nonsystematic errors in lacunarity and L-function values, which are dependent on the width, thickness, and position of each channelized sandbody.

The measured, apparent sandbody widths are larger than true sandbody widths, except where sandbody axes are perpendicular to the cliff-face panels (e.g., Lorenz et al. 1985). Apparent sandbody width (Wa) is dependent on the true width of sandbodies (Wt) and the angle of the outcrop face relative to the orientation of the sandbody axis (θ):

 
formula

In the absence of paleocurrent data or detailed reconstructions of internal sandbody architecture (e.g., Hampson et al. 2013), it is not possible to accurately convert apparent widths to true widths for individual sandbodies. Sandbody orientation(s) relative to the cliff-face panel that contains them affects the calculated value of lacunarity, as discussed later. The orientations of the cliff-face panels are given in Table 1. Regional paleogeographic reconstructions indicate that the Star Point and Blackhawk shorelines were nearly linear and oriented approximately north–south (Fig. 1A; Hampson 2010; Hampson et al. 2011, 2012). Channelized fluvial sandbodies were likely oriented west–east on average over the whole outcrop belt, positioned perpendicular to the regional paleoshoreline trend, but with significant local variability.

Spatial Statistical Methods for the Analysis of Sandbody Distribution

A number of techniques have been applied to interpret the spatial distribution of points and objects in the biological, ecological, and geological sciences, including the quadrat method (Greig-Smith 1952), nearest neighbor distance method (Clark and Evans 1954), Ripley's K-function (Ripley 1977; Cressie 1991), Besag's L-function (Besag 1977), lacunarity (Plotnick et al. 1996), and compensation index (Straub et al. 2009; Wang et al. 2011). In this study, we use lacunarity and Besag's L-function (Besag 1977) to quantitatively characterize the spatial distribution of channelized fluvial sandbodies. The quadrat and nearest neighbor distance methods are inappropriate for this analysis, because they rely on locational information (i.e., quadrat location; Greig-Smith 1952) and permit comparison of results only for areas of similar size (Clark and Evans 1954), respectively. Compensation index quantifies the spatial variation of sediment-body distribution between a chronologically ordered series of depositional horizons (Straub et al. 2009), such as paleosols that are situated adjacent to channelized fluvial sandbodies. Only a small number of such depositional surfaces can be confidently mapped in the study dataset, which makes the application of compensation index problematic.

Lacunarity is a raster-based approach based on the multiscale analysis of spatial dispersion, and characterizes the distribution of gaps or spaces within a pattern as a function of scale (Allain and Cloitre 1991; Plotnick et al. 1996). Besag's L-function, a standardized version of Ripley's K-function, is a point-based method that determines the spatial scales at which points tend to be more or less dispersed than expected by chance (Besag 1977; Besag and Diggle 1977; Ripley 1977; Cressie 1991). Both methods have been used previously on geologic datasets, and can be applied readily to one-, two-, or three-dimensional datasets to give straightforward and interpretable graphical results (Plotnick et al. 1996; Plotnick 1999; Rankey 2002; Hajek et al. 2010; Roy et al. 2010; Zhao et al. 2011). The application of the two methods in this study is described below (Fig. 4).

Lacunarity.—

Lacunarity was selected for this study because it enables the length scales of randomness, clustering, or hierarchical structure to be determined in all of the 2D cliff-face panels (cf. Plotnick et al. 1993). In order to measure lacunarity in different paleogeographic locations and at different stratigraphic levels, the six cliff-face panels (A–F in Fig. 3) and their stratigraphic subdivisions (A1, B1-2, C1-3, D1-3, E1-3, F1-4 in Fig. 3) were analyzed. A binary, black-and-white image of each 2D cliff-face panel or one of its constituent stratigraphic subdivisions was first generated (Fig. 4A). A number of algorithms for computing lacunarity (e.g., gliding-box algorithm, standard box-counting algorithm, and differential box-counting algorithm) have been developed by various authors. The gliding-box algorithm (Allain and Cloitre 1991; Plotnick et al. 1996) was used in this study because it generates clear and interpretable results and is conceptually straightforward, and the large sample size can lead to better statistical results (Dale 2000; Saa et al. 2007). A box of a given length is placed at the top left of the binary image of the cliff-face panel or its stratigraphic subdivision, and the number of pixels representing sandstone (in black) within the box are subsequently scanned and counted. The box is then moved one column along to the right and the process is repeated until the entire image is scanned. In this study, 12 grid box sizes were applied to each panel and its subdivisions. The minimum and maximum box sizes were 2% and 45% of the total image area, respectively. The box size was chosen not to exceed 50% of the image area, because this may introduce point statistical errors (Karperien 2007).

One of the major uses of lacunarity is to compare spatial heterogeneity at different length scales (e.g., different grid box sizes), but little variation in lacunarity exists for the different grid box sizes used in this study. Instead, spatial heterogeneity is averaged across the length scales of all 12 grid box sizes to generate a mean value of lacunarity for each panel or its stratigraphic subdivision. Mean lacunarity was calculated using the following equation:

 
formula

This equation calculates the mean value for lacunarity (λ) for each panel based on the variation in pixels per box summarized over all grid orientations (Rasband 1997–2014; Karperien 2007). (F) refers to the total number of pixels that are considered as “sandbody” pixels (black in Fig. 4A) in the scanned part of the image per box count. Since lacunarity is a scale-dependent measure of heterogeneity (Allain and Cloitre 1991), a low value of lacunarity (minimum  =  0) indicates a low range of gap sizes that are homogeneous and translationally invariant (Fig. 4C). In comparison, patterns with high lacunarity values (maximum  =  1) are heterogeneous, display translational variance, and contain a wide range of gap sizes.

The values of lacunarity measured in this study are expected to vary with the three-dimensional (3D) orientation of channelized fluvial sandbodies relative to the 2D cliff-face panel that they intersect. This expected variation has two components. First, variation in the 3D orientations of the sandbodies relative to each other will result in a range of apparent sandbody widths, which cannot be deconvolved from variation in true sandbody widths, in a 2D cross section. Variation in relative sandbody orientation will also affect the horizontal spacing of the sandbodies in a 2D cross section. This first component of variation is inherent to the population of channelized fluvial sandbodies within each cliff-face panel, and is what we are seeking to measure using lacunarity. Second, the orientation of the 2D cliff-face panel introduces a systematic variation in apparent sandbody width and horizontal spacing that applies to all sandbodies intersecting the panel. This second, systematic component of variation is expected to be more pronounced for cliff-face panels that are oriented oblique to regional depositional strike (e.g., panels C, E, F oriented at N050, at 50° relative to the north–south-trending regional paleoshoreline; Table 1) than for panels oriented along regional depositional strike (e.g., panels A, B, oriented at N000, parallel to the north–south-trending regional paleoshoreline; Table 1). Although this second component of variation is expected to affect the value of lacunarity measured for each cliff–face panel, any trend in lacunarity between the stratigraphic subdivisions of the Blackhawk Formation within each panel (A1, B1-2, C1-3, D1-3, E1-3, and F1-4 in Fig. 3) should not be affected.

Ripley's K-function and Besag's L-function.—

The most common second-order spatial statistical technique for analyzing point patterns is Ripley's K-function, which determines how point pattern distributions change over different length scales in a dataset (Ripley 1977). In order to analyze sandbody centroid distributions using the L-function, a standardized version of Ripley's K-function, a series of panels showing the centroids of channelized fluvial sandbodies were generated (Figs. 4B, 5) for each of the six cliff-face panels (A–F in Fig. 3) and their stratigraphic subdivisions (A1, B1-2, C1-3, D1-3, E1-3, F1-4 in Fig. 3).

The K-function compares the predicted number of points (e.g., sandbody centroids) within a distance (h) of each point in each 2D area (e.g., cliff-face panel) to the average rate of the point process (λ) as outlined in the following equation:

 
formula

where λ is the number of centroid points in the 2D area of radius h (N) divided by the area of the study region, and E(N(h)) is the expected number of points in the same region (Cressie 1991). Besag's L-function (Besag 1977) is used here so that the K-function can be compared with its expected value and against a benchmark of zero (Rosenberg and Anderson 2011):

 
formula

If the expected value of L(h) found at a certain distance is equal to the number of points estimated, taking into account the intensity of the point process, then the distribution pattern is close to zero and represents complete spatial randomness (forumla; Besag 1977. If forumla < 0, then points are clustered, and the L-function plots negatively below the Monte Carlo simulation envelope (Fig. 4D; Besag 1977). In contrast, if forumla > 0, then points are regularly spaced, and the L-function plots positively above the Monte Carlo simulation envelope (Fig. 4D; Besag 1977). The range of critical values that define the minimum and maximum limits of complete spatial randomness are computed by Monte Carlo simulation algorithms (Besag and Diggle 1977). Randomization tests were run on each cliff-face panel and its stratigraphic subdivisions based on 99 repetitions (cf. Hajek 2010). The randomization tests involve the construction of a new point pattern within the confines of the study plot. Each of the new points were analyzed and the upper and lower 95% confidence limits are plotted in combination with the values obtained from analysis of the original coordinates (Monte Carlo significance level, α  =  0.05). Edge effects are generated when centroid points lie in close proximity to the boundaries of the area of study (Goreaud and Pélissier 1999). Several methods have been devised to correct for edge effects, but here we use a method that rescales counts based on the degree to which the area of a circle of radius (h) overlaps the boundary of the study area (Rosenberg and Anderson 2011).

The cliff-face panels (A–F in Fig. 3) and their stratigraphic subdivisions (A1, B1-2, C1-3, D1-3, E1-3, F1-4 in Fig. 3) were vertically exaggerated by ×55.8, which corresponds to the ratio of mean sandbody width to mean sandbody thickness across all six cliff-face panels, in order to minimize the effects of anisotropy in sandbody dimensions on the results of our L-function analysis. The advantages of attempting to introduce statistical isotropy to the study dataset are that data points display a constant mean and a constant variance in all directions, and that covariance is dependent only on distance (Masihi et al. 2006). The maximum distance between points that is considered in our application of the L-function is 25% of the vertically exaggerated thickness of the cliff-face panel or its stratigraphic subdivision, in order to avoid distortion by edge effects (Rosenberg and Anderson 2011). As a result, only centroid distributions at length scales smaller than this distance can be investigated. The position of sandbody centroids in the plane of a 2D cliff-face panel are independent of the 3D orientation of the sandbodies, but their horizontal spacing may vary according to the orientation of the cliff-face panel. The total number of centroids in the stratigraphic subdivisions of each cliff-face panel ranges from 12 (F2; Table 1) to 72 (A1; Table 1). The L-function is sensitive to the number of points and the size of the study area, such that the results for some stratigraphic subdivisions and associated cliff-face panels may not be statistically robust.

Results and Analysis

494 channelized fluvial sandbodies have been identified in the studied cliff-face panels (A–F; Fig. 3, Table 1). The limited grain-size variability within each sandbody prevents consistent interpretation of internal sandbody architectures in detail with the available photographic data. As a result, a rigorous architectural analysis of the sandbodies based on the hierarchical arrangement of basal and internal erosion surfaces (Holbrook 2001) is not attempted here. However, individual sandbodies whose internal structure can be interpreted from the available photographic data exhibit single-story, multilateral, and multistory architectures (Fig. 6), consistent with more thorough analysis of selected parts of the Blackhawk Formation outcrop belt on the ground (Adams and Bhattacharya 2005; Hampson et al. 2013).

Sandbody Dimensions

The apparent width and true thickness of each sandbody were measured from cliff-face panels A–F (Fig. 7) and for the four stratigraphic subdivisions of the Blackhawk Formation within each panel (Fig. 8).

Description.—

The mean apparent width of channelized fluvial sandbodies over the entire dataset is 380 m (standard deviation of 430 m), and the mean sandbody thickness is 7 m (standard deviation of 3 m). Large variations of apparent sandbody widths and associated large values of standard deviation are not dependant on paleogeographic location (Fig. 7) or stratigraphic interval (Fig. 8).

Panels C, D, E, and F display larger mean apparent sandbody widths (mean apparent widths of 750, 360, 360, and 370 m respectively; Fig. 7C–F) than panels A and B (mean apparent widths: 250 and 190 m; Fig. 7A, B). Sandbodies also tend to be thinner in panels A and B (mean thicknesses of 4 m; Fig. 7A, B) than in panels C–F (mean thicknesses of 7–10 m; Fig. 7C–F).

The apparent widths of channelized fluvial sandbodies increases from the “lower Blackhawk Formation” interval (mean apparent width of 350 m, Fig. 8A), to the “upper Blackhawk Formation 2” interval (mean apparent width of 420 m, Fig. 8C). The “upper Blackhawk Formation 3” interval contains sandbodies of smaller apparent width (mean apparent width of 390 m), but this value is based on a small dataset (n  =  13, Fig. 8D). Sandbodies situated in each stratigraphic interval (Fig. 8A–D) display larger apparent widths where they intersect the outcrop belt in cliff faces oriented at N050 and N030 (gray and white bars in Fig. 8), compared to cliff faces oriented at N010 (black bars in Fig. 8). There is no apparent stratigraphic trend in sandbody thickness (mean thicknesses of 6–8 m; Fig. 8).

Interpretation.—

The larger mean apparent sandbody widths in panels C, D, E, and F (Fig. 7C–F) than in panels A and B (Fig. 7A, B) are attributed to the more oblique orientation of these former panels (N050, N030; Table 1) relative to the interpreted north–south-trending regional paleoshoreline (Figs. 1A, 2). These results support the hypothesis that channelized sandbodies form a population with a mean west–east azimuth, normal to the overall paleoshoreline trend, over the whole outcrop belt. Estimates of mean true sandbody widths are generated by correcting for the orientation of the cliff-face panels relative to the assumed west–east mean sandbody azimuth; for panels A–F, these estimates of mean true sandbody width are 240 m, 190 m, 480 m, 350 m, 230 m, and 240 m. The absence of pronounced paleogeographic variation in estimated true sandbody width between the six cliff-face panels implies the presence of multiple feeder trunk rivers up-dip (west) of the outcrop belt with relatively uniform dimensions and lateral migration behaviors.

The upward increase in the mean apparent width of channelized sandbodies in the Blackhawk Formation (Fig. 8) is consistent with previous studies (Marley et al. 1979; Hampson et al. 2012). This trend suggests an increase in lateral channel migration and/or widening of channel belts with decreasing tectonic subsidence rate and increasing distance from the coeval shoreline (cf. Allen 1978; Shanley and McCabe 1994; Hampson et al. 2012).

Sandbody Distributions

Lacunarity and L-function were calculated to characterize sandbody distributions in cliff-face panels A–F (Figs. 9A, 10), and for the four stratigraphic subdivisions of the Blackhawk Formation in each panel (Figs. 9B–D, 11, 12, and Appendix Fig. 16).

Description of Lacunarity Results.—

Values of lacunarity for individual cliff-face panels range between 0.24 (Panel E) and 0.61 (Panel D) (filled circles in Fig. 9A). There is no correlation between panel location and lacunarity value. However, there is a weak negative correlation between the number of sandbodies per unit area and the south-to-north location of a panel (R2  =  0.31; red crosses in Fig. 9A) and a moderate positive correlation between the proportion of sandstone (i.e., net-to-gross ratio) and south-to-north panel location (R2  =  0.61; gray open circles in Fig. 9A). Lacunarity does not exhibit any relationship to apparent width of channelized fluvial sandbodies in the panels (Fig. 9A).

Panels A–F were subsequently separated into stratigraphic subdivisions to allow temporal and spatial trends to be interpreted (panel subdivisions A1, B1-2, C1-3, D1-3, E1-3, and F1-4; Fig. 3). Values of lacunarity show an overall decrease from lower to higher stratigraphic positions (R2  =  0.60; black filled circles in Fig. 9B). The overall upward decrease in lacunarity is associated with an increase in net-to-gross ratio (R2  =  0.68; gray unfilled circles in Fig. 9B), a decrease in the number of sandbodies per unit area (R2  =  0.51, red crosses in Fig. 9C), and an increase in the apparent width of channelized fluvial sandbodies (R2  =  0.50, blue squares in Fig. 9D).

Description of L-Function Results.—

Cliff-face panels A, C, D, E, and F display clustering of sandbody centroids over length scales that lie between c. 0.5 and 15 times the mean sandbody dimensions (mean apparent sandbody width and thickness of 380 m and 6.8 m, respectively) (Fig. 10A, C–F). Cliff-face panel B displays a random distribution of sandbody centroids over length scales of up to six times the mean sandbody dimensions (Fig. 10B).

In the “lower Blackhawk Formation” and “upper Blackhawk Formation 1” intervals, sandbody-centroid clustering occurs over length scales that lie between c. 0.5 and 7 times the mean sandbody dimensions (Figs. 11A–D, F, 12). Sandbody centroids are regularly spaced at a length scale of c. four times the mean sandbody dimensions in the “lower Blackhawk Formation” interval of panel E (Fig. 11E). In the “upper Blackhawk Formation 2” and “upper Blackhawk Formation 3” intervals, limited clustering of sandbody centroids occurs over length scales that lie between c. 0.5 and 2 times the mean sandbody dimensions (Fig. 12B and Appendix Fig. 16).

Interpretation.—

Lacunarity and L-function results suggest that localized sandbody clustering is more prominent (Fig. 12B) and stratigraphic architecture displays greater heterogeneity (Fig. 9B–D) towards the base of the Blackhawk Formation. Towards the top of the Blackhawk Formation, lacunarity and L-function results show a decrease in sandbody clustering and a corresponding increase in regularity of sandbody spacing (Fig. 12B), and that stratigraphic architecture is more homogeneous (Fig. 9B–D). The spatial organization of channelized fluvial sandbodies throughout the Blackhawk Formation in general exhibits more clustering than regularity, although random distributions are most abundant (Fig. 12B).

The absence of a consistent trend in lacunarity or sandbody-centroid distribution from cliff-face panels A to F (Figs. 9A, 10, 12A) suggests that there was no large-scale paleogeographic control on the distribution of channelized fluvial sandbodies. Trends are more apparent when the panels are broken down into stratigraphic subdivisions, implying a stratigraphic control on patterns of sandbody distribution, as reflected in lacunarity, net-to-gross ratio (Fig. 9A), mean apparent sandbody width (Fig. 9D), and patterns of sandbody-centroid spacing (Fig. 12B).

The overall upward decrease in values of lacunarity from the base to the top of the Blackhawk Formation is accompanied by increases in both net-to-gross ratio (Fig. 9A) and apparent sandbody width (Fig. 9D). Consequently, the variance of gaps between sandbodies decreases and stratigraphic architecture becomes more homogeneous. The upward increases in homogeneity and regularity of stratigraphic architecture at length scales smaller than the cliff-face panels likely result from autogenic, avulsion-related patterns of sandbody stacking that were controlled by the large-scale decrease in tectonic subsidence rate and increase in distance from the coeval shoreline.

Discussion: Comparison With Numerical Modeling Experiments

In order to aid interpretation of the spatial patterns of distribution of channelized fluvial sandbodies in the Blackhawk Formation dataset, we used lacunarity and L-function to characterize the distributions of similar sandbodies in three widely cited numerical modeling studies (Bridge and Leeder 1979; Heller and Paola 1996; Karssenberg and Bridge 2008) that studied the effects of channel-belt size relative to floodplain width, avulsion frequency and associated local sedimentation rate, and base-level fluctuations on alluvial stratigraphy. These modeling studies provide a benchmark for interpreting outcrop data, because simulated stratigraphic output from the models can be linked directly to mechanisms, boundary conditions, and input parameters. The models are generic, and simulate mechanisms that are inferred to have operated during deposition of the Blackhawk Formation, albeit with boundary conditions and input parameters that may be different. The modeling results presented below are ordered by increasing complexity of the numerical model, and are used to infer the plausibility of different mechanisms responsible for sandbody-distribution patterns in the Blackhawk Formation.

Simple 2D Model of Autogenic, Avulsion-Generated Architectures (Bridge and Leeder 1979)

The numerical model of Bridge and Leeder (1979) simulates avulsion during aggradation of fluvial strata using the assumptions that floodplain sedimentation accumulation rate decreases away from an active fluvial channel, and that newly avulsed channels occupy the topographically lowest part of the floodplain. Avulsion frequency is assumed to be spatially and temporally constant. These assumptions, in combination with differential compaction of channelized sandbodies and floodplain shales, result in architectures that are characterized by compensational stacking (cf. Straub et al. 2009; Wang et al. 2011). Figure 13A–C shows stratigraphic cross sections from three model runs, for a uniform number of channelized fluvial sandbodies of different dimensions and a floodplain of uniform width (modified from fig. 2 of Bridge and Leeder 1979).

L-function analysis of each cross section indicates regular spacing of sandbody centroids (Fig. 13D–F). The spacing of sandbody centroids, normalized to mean sandbody dimensions, decreases as sandbody dimensions increase; this result arises from the uniform number of sandbodies in each cross section (Figs. 12C, 13D–F). The cross sections all display low to moderate values of lacunarity (0.17–0.36; purple bars in Fig. 12C), consistent with a low degree of heterogeneity in sandbody spacing. The net-to-gross ratio for the cross sections illustrated in Figure 13A, B, and C are 0.33, 0.51, and 0.88, respectively. These results indicate that a combination of low to moderate values of lacunarity and regularly spaced sandbody centroids is characteristic of stratigraphic architectures generated by compensational avulsions. Regular sandbody spacing due to compensational stacking occurs over relatively short length scales (< 1–2 times the mean sandbody dimensions, for different model runs; Fig. 13), due to the relatively large number of sandbodies per unit area and relatively high net-to-gross ratio. In the case with the highest net-to-gross ratio, sandbody clustering also occurs, but over larger length scales (> 2 times the mean sandbody dimensions, Fig. 13C).

Regular spacing of sandbody centroids is weakly but consistently developed in the “upper Blackhawk Formation 2” and “upper Blackhawk Formation 3” intervals (red and yellow bars in Fig. 12B). Comparison with the models of Bridge and Leeder (1979) implies that this regular centroid spacing may have resulted from avulsion-generated compensational stacking of channelized fluvial sandbodies, to infill differential topography on the upper coastal plain and alluvial plain (cf. Straub et al. 2009; Wang et al. 2011). This inference is supported by the low values of lacunarity calculated for these stratigraphic intervals (red and yellow bars in Fig. 12B) and model cross sections (purple bars in Fig. 12C). Regular spacing of sandbody centroids is less common than clustering of sandbody centroids in the “lower Blackhawk Formation” and the “upper Blackhawk Formation 1” intervals (green and blue bars in Fig. 12B), implying that compensational stacking was less pronounced during deposition on the lower coastal plain, at the length scales of the cliff-face panel subdivisions (A1, B1-2, C1-3, D1-3, E1-3, F1-4 in Fig. 3). This contrast in sandbody distributions suggests that the boundary conditions during deposition of the “lower Blackhawk Formation” and “upper Blackhawk Formation 1” intervals (0–120 km from coeval shoreline, and sediment accumulation rates of 80–200 m/Myr; Fig. 2) were less conducive to compensational stacking by avulsion than the boundary conditions during deposition of the “upper Blackhawk Formation 2” and “upper Blackhawk Formation 3” intervals (40–130 km from coeval shoreline, and sediment accumulation rates of 40–700 m/Myr; Fig. 2). The degree of compensational stacking may also have been influenced by changes in: (1) avulsion frequency, (2) the ratio of sandbody width relative to floodplain width, (3) the ratio between local sedimentation and subsidence rates, and (4) paleotopography induced by peat formation and compaction (cf. Straub et al. 2009; Hofmann et al. 2011).

Complex 2D Model of Autogenic, Avulsion-Generated Architectures (Heller and Paola 1996)

Using a numerical model of aggradational alluvial architecture, Heller and Paola (1996) investigated the influence of the relationship between avulsion frequency (Fa) and local sedimentation rate (r) on stratigraphic architecture. We consider three of the cases modeled by these authors, in which avulsion frequency increases at slower rates than sedimentation rate (e.g. Fig. 14B), increases linearly with sedimentation rate (e.g., Fig. 14A, equivalent to the model of Bridge and Leeder 1979), or increases faster than sedimentation rate (e.g., Fig. 14C). Sandbody stacking density varies according to the relationship between avulsion frequency and sedimentation rate (e.g., compare Fig. 14A, B, and C). Figure 14A–C shows stratigraphic cross sections from three model runs in which this relationship differs, for the upstream part of a back-tilted foreland basin (modified from figure 9 of Heller and Paola 1996), comparable in tectonic setting to the Western Interior Basin foredeep, in which the study area was located (Fig. 1). L-function analysis of the cross sections indicates regular spacing of sandbody centroids for avulsion frequency increasing linearly with sedimentation rate (Fa  =  r1; net-to-gross ratio: 0.29, Fig. 14D), slower than sedimentation rate (Fa  =  r0.2; net-to-gross ratio: 0.26, Fig. 14E), or faster than sedimentation rate (Fa  =  r1.5; net-to-gross ratio: 0.37, Fig. 14F). Regular spacing is more pronounced when avulsion frequency increases faster than sedimentation rate (Fig. 14F). The three cross sections have relatively low values of lacunarity (0.03–0.27; brown bars in Fig. 12C), consistent with homogeneous sandbody dimensions and regular spacing.

Analysis of Heller and Paola's (1996) model indicates that regular spacing of channelized fluvial sandbodies, most likely due to compensational stacking (cf. Bridge and Leeder 1979), is generated by avulsion for a range of relationships between avulsion frequency and local sedimentation rate, and that similar lacunarity and sandbody-centroid distribution patterns can be generated independent of the relationship between avulsion frequency and local sedimentation rate (brown bars in Fig. 12C). These model results therefore imply that regular spacing of sandbody centroids, which is pronounced in the “Upper Blackhawk Formation 2” and “upper Blackhawk Formation 3” intervals (red and yellow bars in Fig. 12B) but also occurs sparsely in the “lower Blackhawk Formation” and “upper Blackhawk Formation 1” intervals (green and blue bars in Fig. 12B), can be generated by avulsion for a range of relationships between avulsion frequency and local sedimentation rates. The models of Heller and Paola (1996) further imply that the relationship between avulsion frequency and local sedimentation rate is important in controlling sandbody thickness and the degree of vertical stacking indicated by internal sandbody architecture; thicker, more aggradational sandbodies are developed where avulsion frequency increases more slowly than local sedimentation rate (Fig. 14B). Sandbody thickness shows a similar mean and distribution in each of the stratigraphic subdivisions of the Blackhawk Formation (Fig. 8), which implies the absence of a systematic stratigraphic control on the relationship between avulsion frequency and local sedimentation rate. There is some paleogeographic variation in the mean and range of sandbody thickness, with thicker sandbodies occurring in cliff-face panels C–F than in panels A and B (Fig. 7), which may indicate that avulsion frequency increased more slowly than local sedimentation rate towards the north of the study area. The apparently weak dependence or independence of stratigraphic architecture in the Blackhawk Formation on the relationship between avulsion frequency and local sedimentation rate is supported by analysis of well-exposed fine-grained floodplain deposits, which indicates that avulsion style did not vary significantly with stratigraphic position or paleogeographic location, and associated changes in local and regional boundary conditions (Flood and Hampson 2014).

Complex 3D Model of Allogenic and Autogenic, Avulsion-Generated Architectures (Karssenberg and Bridge 2008)

A complex 3D numerical model of channel-belt network and floodplain deposition in a delta is described by Karssenberg and Bridge (2008). The model simulates channel bifurcation and avulsion, driven by topographic gradients that result from local sedimentation histories and by thresholds in river discharge. Various model runs have been used to assess allogenic and autogenic controls on delta-plain dynamics and stratigraphic architecture, using input parameters derived from the Holocene-to-modern Rhine–Meuse delta (Karssenberg and Bridge 2008). Figure 15 shows 3D perspective views and 2D cross sections from two model runs, for uniformly rising base level (Fig. 15A, B; modified from figure 4 of Karssenberg and Bridge 2008) and for a cycle of base-level fall and subsequent rise (Fig. 15D, E; modified from figure 9 of Karssenberg and Bridge 2008). Both model runs contain avulsion nodes, most prominently at the delta apex near the modeled inflow point, and a network of downstream-narrowing channel belts (Fig. 15A, D). The model run for a base-level fall-then-rise cycle features the development, filling, and overtopping of incised valleys. The frequency of channel avulsion was also higher during base-level rise and associated high aggradation rate (Karssenberg and Bridge 2008).

Both model runs feature pronounced clustering of sandbody centroids, over distances of c. 0.5–7.5 times mean sandbody dimensions, as indicated by L-function analysis of three cross sections in upstream-to-downstream locations that lie downstream of the avulsion node at the delta apex (Fig. 15C, F). There is little difference between the results of L-function analysis for the three cross sections in either model run (Fig. 15C, F). Values of lacunarity occupy a narrower range in the base-level-rise model run (0.19–0.21; maroon bars in Fig. 12C) than in the base-level fall-then-rise model run (0.15–0.31; orange bars in Fig. 12C). The cross sections for the base-level-rise model run have net-to-gross ratios of 0.33–0.37 (Fig. 15B), and corresponding net-to-gross ratios for the base-level fall-then-rise model run are 0.20–0.27 (Fig. 15E).

These results indicate that stratigraphic architecture in delta-plain settings is dominated by avulsion-related clusters of channelized sandbodies that form in response to steady sea-level rise (Fig. 15A–C) or to cyclical sea-level rises and falls (Fig. 15D–F). Stratigraphic architecture becomes more heterogeneous from upstream to downstream, probably as a result of downstream narrowing of sandbody dimensions and increasing avulsion frequency downstream of the avulsion node (Mackey and Bridge 1995). The wider range of values of lacunarity for the base-level fall-then-rise model run (Fig. 15D–F) than for the base-level rise model run (Fig. 15A–C) suggest that a downstream allogenic control increases heterogeneity in stratigraphic architecture, despite a similar degree of sandbody clustering in the two model runs (compare Fig. 15F and C). The more heterogeneous architecture of the base-level fall-then-rise model run probably results from a larger range of sandbody dimensions and more ordered vertical distribution of sandbodies, relative to the base-level-rise model run (compare Fig. 15E and orange bars in Fig. 12C, with Fig. 15B and maroon bars in Fig. 12C).

The localized clustering of sandbody centroids and generally high values of lacunarity that characterize the “lower Blackhawk Formation” and “upper Blackhawk Formation 1” stratigraphic intervals in the studied cliff-face “windows” of the Blackhawk Formation (green and blue bars in Fig. 12B) are comparable to the cross sections in the delta-plain models of Karssenberg and Bridge (2008) (maroon and orange bars in Figs. 12C). Sandbody clustering in these models arises as a result of channel bifurcation and avulsion on the delta plain, which is fully consistent with interpretation of the “lower Blackhawk Formation” and “upper Blackhawk Formation 1” intervals as coastal-plain deposits, developed 0–120 km (but generally < 50 km) inland of paleoshorelines that comprised a series of wave-dominated and wave-influenced deltas with adjacent strandplains and barrier islands (Flores et al. 1984; Dubiel et al. 2000; Hampson et al. 2011, 2012). Furthermore, the distribution of extensive coal zones that formed during rising relative sea level and interpreted incised valleys that developed during falling and lowered relative sea level indicate that low-amplitude (< 30 m), high-resolution (< 0.1 Myr) cycles in relative sea level influenced stratigraphic architecture in the “lower Blackhawk Formation” interval (Flores et al. 1984; Dubiel et al. 2000; Hampson et al. 2012).

Wider Implications

Although multiple conceptual models predict patterns of fluvial sandbody dimensions and spatial distributions (e.g., Wright and Marriott 1993; Shanley and McCabe 1994; Richards 1996; Weissmann et al. 2013), this study is one of only a handful that quantitatively characterizes stratigraphic patterns in appropriately large outcrop datasets or data compilations of ancient systems (e.g., Hajek et al. 2010; Colombera et al. 2012; Pranter et al. 2013). The results suggest a template for spatial patterns of distribution of fluvial sandbodies that largely reflects autogenic behaviors, and that may be generally applicable to major (i.e., 1–10 Myr duration) regressive clastic wedges, at least in back-tilted foreland basins (e.g., Heller and Paola 1996; Hajek et al. 2014). Lower-coastal-plain strata are characterized by sandbody clustering that resulted from the internal dynamics of delta-plain distributary channel systems, potentially in combination with externally forced relative sea-level changes (e.g., Karssenberg and Bridge 2008) (Fig. 15). Upper-coastal-plain and alluvial-plain strata exhibit less pronounced sandbody clustering and more pronounced regular spacing, the latter due to compensational stacking of channel-belt deposits that were not constrained by the position of long-lived avulsion nodes at delta apices in upstream locations (e.g., Bridge and Leeder 1979) (Fig. 13). These spatial patterns appear to occur over a wide range of estimated sediment accumulation rates, which approximate tectonic subsidence rates (c. 40–700 m/Myr; Fig. 2), suggesting that there is no strong relationship between avulsion frequency and sedimentation rate (e.g., Heller and Paola 1996) (Fig. 14). Mean sandbody width increases only slightly from base to top of the Blackhawk Formation (Fig. 8), implying that there was only modest variation in the dimensions and lateral migration of paleochannels, as recorded in the width and thickness of channel-belt sandbodies. This inference is consistent with detailed architectural studies of selected channelized fluvial sandbodies in the Blackhawk Formation and the Castlegate Sandstone, which document only modest changes in fluvial style (Adams and Bhattacharya 2005; Hampson et al. 2013). High-frequency, allogenic cycles in relative sea level appear to have only a minor impact on sandbody clustering (cf. Karssenberg and Bridge 2008) (Fig. 15), outside of incised-valley fills that form thick (up to 25 m), laterally extensive (1–6 km), multistory and multilateral sandbodies (e.g., figure 13 in Hampson et al. 2012). Evidence of high-frequency, allogenic upstream controls (e.g., climate cycles, or “pulses” of tectonic uplift) also appears to be absent, because such signals have a subtle expression and/or have been lost during sediment transport (“signal shredding” sensuJerolmack and Paola 2010).

This stratigraphic-architectural template of distributions of fluvial sandbodies has several implications for the connectivity of hydrocarbon reservoirs and groundwater aquifers. The connectivity of channelized fluvial sandbodies is generally considered to be controlled by: (1) the proportion of channelized sandbodies in the reservoir (cf. net-to-gross ratio); (2) sandbody width and thickness; (3) sandbody sinuosity and range of sandbody orientations; (4) organization of sandbody stacking; and (5) the sandstone content of crevasse-splay and other nonchannelized floodplain deposits (Jones et al. 1995; Larue and Hovadik 2006). The stratigraphic-architectural template principally addresses the organization of sandbody stacking. Clustering of channelized fluvial sandbodies in lower-coastal-plain strata is likely to increase sandbody connectivity within the clusters (e.g., branching network of deltaic distributary-channel sandbodies) (Larue and Hovadik 2006), although the clusters themselves may be poorly connected. Compensational stacking and the resultant regular spacing of channelized fluvial sandbodies in upper-coastal-plain and alluvial-plain strata is likely to reduce sandbody connectivity (Larue and Hovadik 2006), although this may be offset by increases in the proportion of sandbodies and in mean sandbody width in these strata. These autogenically generated sandbody distributions occur in a predictable arrangement in the Blackhawk Formation, and potentially in other major (i.e., 1–10 Myr duration) regressive clastic wedges.

Conclusions

Data from a large, well-exposed outcrop belt (Upper Cretaceous Blackhawk Formation, Wasatch Plateau, central Utah, U.S.A.) are used to quantitatively characterize and analyze the geometry and distribution of channelized fluvial sandbodies in coastal-plain strata. The abundance and apparent dimensions of channelized fluvial sandbodies broadly increase from base to top of the Blackhawk Formation, resulting in an upward-increasing proportion of fluvial sandbodies (cf. net-to-gross ratio). These trends suggest an increase in lateral channel migration and/or widening of channel belts from base to top of the formation. Variations in apparent sandbody width between the cliff-face panels correlate to panel orientations, implying that channelized sandbodies form a population with a mean west–east azimuth, approximately perpendicular to the regional paleoshoreline trend. Lacunarity decreases from base to top of the Blackhawk Formation, indicating that a wider range and variance of gaps occurs between sandbodies as the proportion of channelized sandbodies and apparent sandbody width increase. Localized clustering of sandbody centroids also occurs preferentially in the lower part of the Blackhawk Formation, which comprises lower-coastal-plain strata (< 50 km from the coeval shoreline). Sandbody centroids show a weak tendency for regular spacing in the upper part of the formation.

Comparison with generic numerical modeling studies suggests that spatial patterns of sandbody distribution can be attributed to: (1) avulsion of deltaic distributary channels in locations downstream of long-lived avulsion nodes, which may also have been modulated by high-frequency relative-sea-level cycles, in the lower part of the Blackhawk Formation; and (2) compensational stacking of sandbodies in the upper part of the Blackhawk Formation. Sediment accumulation rate varied significantly during deposition (c. 40–700 m/Myr), but any potential variation in its relationship to avulsion frequency had little influence on avulsion style or, via comparison with numerical modeling studies, on sandbody-distribution patterns. These results imply that autogenic behaviors were the dominant control on stratigraphic architecture and patterns of sandbody distribution.

This study demonstrates the value of collecting large outcrop datasets, which enable quantitative characterization of sandbody distributions and related stratigraphic architectures using spatial statistical methods. The results also provide validation of numerical modeling studies of avulsion mechanisms and controls using exposures of ancient fluvial strata.

The authors thank the Department of Earth Science and Engineering, Imperial College London for support of YSF via a Janet Watson PhD Scholarship, and Chevron Energy Technology Company for additional support. We thank Bryan Bracken, Brian Willis, Michael Pyrcz, and Elizabeth Hajek for insightful discussions, and Ellen Chamberlin, Sian Davies-Vollum, and Gary Weissmann for their constructive review and editorial comments. Image J and FracLac were used for lacunarity analysis and PASSaGe v2 to apply Besag's L-function. We acknowledge Wayne Rasband for developing ImageJ (Research Services Branch, National Institute of Mental Health, USA), Audrey Karperien for creating the FracLac plugin for ImageJ (Charles Sturt University, Australia), and Michael Rosenberg and Corey Anderson for developing PASSaGE 2 (Arizona State University, USA). We appreciate the authors making the software readily available for public use.

Appendix Fig. 16.—