Reconstructing subsurface sandbody connectivity from temporal evolution of surface networks

Characterization of groundwater aquifers and hydrocarbon reservoirs requires an understanding of the distribution and connectivity of subsurface sandbodies. In deltaic environments, distributary channel networks serve as the primary conduits for water and sediment. Once these networks are buried and translated into the subsurface, the coarse‐grained channel fills serve as primary conduits for subsurface fluids such as water, oil or gas. The temporal evolution of channels on the surface therefore plays a first‐order role in the 3D permeability and connectivity of subsurface networks. Land surface imagery is more broadly available than topographic or subsurface data, and time‐series imagery of river networks can hold useful information for constraining the shallow subsurface. However, these reconstructions require an understanding of the degree to which channel bathymetry and river kinematics affect connectivity of subsurface sandbodies. Here, we present a novel method for building synthetic cross sections using overhead images of an experimental delta. We use principal components analysis to extract river networks from surface imagery, then couple this with an inverse‐CDF method to estimate channel bathymetry, to generate a time‐series of synthetic delta topography. This synthetic topography is then transformed, accounting for deposition and subsidence, to produce synthetic stratigraphy that differentiates coarse‐grained channel fill from overbank and offshore deposition. We find that large‐scale subsurface architecture is relatively insensitive to details of channel bathymetry, but instead is primarily controlled by channel location and kinematics. We analyse the connectivity of sand bodies and the geometries of barriers to flow and find that periods of rapid sea‐level rise have more variability in sand body connectivity. We also find that barrier width decreases downstream during all sea‐level phases. Our method generates synthetic stratigraphy that closely resembles the large‐scale architecture and 2‐dimensional connectivity of the real stratigraphy built during the experiment it was based on. We anticipate that it will be broadly applicable to other experimental and field‐scale scenarios.


| INTRODUCTION
Constraining the three-dimensional architecture of groundwater aquifers or hydrocarbon reservoirs relies heavily on predicting the vertical and lateral connectivity of sandbodies in the subsurface. In deltaic environments, distributary channel networks serve as the primary conduits for water and sediment. Once these networks are buried and translated into the subsurface, the coarsegrained channel fills serve as likely conduits for subsurface fluids such as water, oil, or gas (Fogg, 1986;Jones & Scott, 1984;Livera, 1989;Sawyer et al., 2015). The 3D permeability and connectivity structure of subsurface networks is therefore largely driven by the temporal evolution of channels on the surface (Friend et al., 1979;Galloway & Sharp, 1998;Larue & Hovadik, 2006;Xu et al., 2021). This temporal evolution involves both horizontal migration of the channels and vertical translation of the channel network through net deposition. In effect, the channel network on the surface is translated into a 3D network in the subsurface via a complex, poorly understood set of transformations involving channel migration (both continuous and abrupt), sediment reworking and net deposition (e.g. Bridge & Leeder, 1979;Fisk et al., 1954;Friend et al., 1979;Leeder, 1978;Sahoo et al., 2020). This paper presents a first step in quantifying these transformations using a well-constrained experimental data set.
The approach to predicting the 3D connection structure in the subsurface starts with observations of channel networks in the modern world. Imagery is the most common and obvious source for information on network geometry. A single snapshot captures only one instant in time and provides no information about migration and channel kinematics. Time series imagery is sometimes available, but still usually captures only a limited range in time. Furthermore, to begin modelling stratigraphy, we need information on the third dimension -channel depth -which is much less available than imagery.
Thus, we identify two major problems in linking channel imagery to 3D subsurface 'plumbing', or connectivity: first, constraining channel depth from 2D imagery, and second, linking static channel networks to network kinematics. We present some initial steps to address both problems in this paper.

| Subsurface connectivity
The size and quality of a reservoir or aquifer depend on the grain size and sorting of the sediment as well as the geometry and connectivity of the various depositional elements (e.g. channel fills, accretion packages, crevasse splays) (e.g. Colombera & Mountney, 2021;Dreyer, 1990;Dreyer et al., 1990;Willis & Tang, 2010). At the most basic level, coarse-grained sedimentary rock such as sandstone or conglomerate is typically porous and can serve as highquality reservoir rock while fine-grained sedimentary rock forms barriers to subsurface fluid flow (Beard & Weyl, 1973). In groundwater systems, these sandbodies serve as aquifers and as conduits for contaminant transport. In deltas with widespread contamination issues, such as the Ganges-Brahmaputra-Meghna delta in Bangladesh (Edmunds et al., 2015;Fendorf et al., 2010), subsurface sedimentary structure has been shown to be a critical factor in sustainability of uncontaminated water (Khan et al., 2016;Mozumder et al., 2020). The focus of this study is sand connectivity in deltaic settings, and in particular, the ways in which the preservation and architecture of deltaic sandbodies in the sedimentary record reflect the temporal evolution of a delta topset (i.e. the subaerial part of a delta). In the topset, the deposition of coarse-grained sediment is primarily associated with distributary channels while floodplain deposition is dominated by mud and silt (Fisk, 1955;Oomkens, 1970). We thus simplify the complex problem of relating surface networks to subsurface architecture by focusing on channel stacking as a primary on. We anticipate that it will be broadly applicable to other experimental and field-scale scenarios.

K E Y W O R D S
aquifer, delta, physical modelling, reservoir, river, sandbody

Highlights
• We present a new method for generating synthetic stratigraphy using overhead images and limited topographic data. • Large-scale architecture of channel sand bodies is driven primarily by river kinematics. • Connectivity of subsurface channel sand bodies is most variable during rapid sea-level rise. • Fine-grained barriers to subsurface fluid flow are widest in proximal reaches of the delta in all sea-level phases.
control on permeability structure (Daws & Prosser, 1992;Doyle & Sweet, 1995). If we can better understand the linkages between channel kinematics and translation of channel bodies into the subsurface, then we can better predict size and connectivity of fluvio-deltaic reservoirs and aquifers.

| Channel depth
Estimating river bathymetry is necessary for hydrodynamic modelling, flood management and for predicting the dimensions of channel bodies that are ultimately preserved in the subsurface. Techniques used to estimate river channel depth from remotely sensed information include interpretation of airborne water reflectance data (Kinzel et al., 2013;Legleiter et al., 2004Legleiter et al., , 2009, operating on hyperspectral imagery with machine learning algorithms (Pan et al., 2015) and developing empirical relationships between channel geometry and catchment area (Leopold & Maddock, 1953). However, our ability to estimate river bathymetry from remotely sensed information remains limited by the availability of field measurements for calibration and the spatial and temporal resolution of satellite data (e.g. Moramarco et al., 2019;Tuozzolo et al., 2019). One common situation is to have limited bathymetry but more frequent overhead imagery because the latter is far easier to obtain. Is there a way to combine these two data sources to get a better idea of depth variability in space and time? We study this problem here, using abundant experimental imagery and high-quality, but much less frequent, topography to (1) provide rough estimates of local channel depth and (2) estimate the extent to which uncertainty in depth degrades estimated subsurface connectivity.

| Channel kinematics
Sandbody connectivity is strongly affected by channel kinematics, that is, the temporal evolution of a channel network. As rivers traverse their floodplains, they migrate laterally through erosion of channel banks and deposition of channel bars, and/or change course via avulsion. Channel mobility, taken here to include both migration (i.e. continuous lateral movement) and avulsion (i.e. abrupt shift in channel location), is influenced by a variety of factors such as sediment flux, discharge, sediment composition and sea-level change (e.g. Bufe et al., 2016Bufe et al., , 2019Constantine et al., 2014;Dunne et al., 2010;Mohrig et al., 2000;Wickert et al., 2013). River avulsions can take many forms, which may build sharp-based sand packages encased by floodplain sediment or instead may develop progradational sandbodies with more gradational transitions with their floodplain deposits (Chamberlin & Hajek, 2015;Hajek & Edmonds, 2014). Rates of lateral migration, frequency and styles of channel avulsion, channel reoccupation and rates of vertical aggradation all play significant roles in the connectivity structure in the subsurface (e.g. Allen, 1965;Bridge & Leeder, 1979;Leeder, 1978;Sheets et al., 2002). For example, high rates of lateral migration with limited vertical aggradation may result in a sandbody with high lateral connectivity but low vertical connectivity. Frequent avulsion may potentially lead to lower overall connectivity; however, the likelihood of a river to revisit previously occupied locations may result in multistorey sandbodies with high vertical connectivity (e.g. Chamberlin & Hajek, 2015;Friend et al., 1979;Gibling, 2006;Heller & Paola, 1996;Mohrig et al., 2000). The major challenge we face in linking surface kinematics to subsurface connectivity is that we generally either have the subsurface record, with limited knowledge of the surface evolution, or we have records of modern systems but limited information on the associated stratigraphic architecture (e.g. cores, radar, shallow seismic). Here, we use the high-resolution imagery and intermittent topographic scans of an experimental delta to build a simplified stratigraphic record that allows us to explore the relationships between surface evolution and subsurface channel sand connectivity. Our particular focus is to explore the potential benefits and drawbacks of combining historical records of channel position with crude constraints on depth to estimate stratigraphic architecture and connectivity.

| XES10 experimental data set
We use a data set from experiments described by  in the Experimental Earthscape (XES) facility at St. Anthony Falls Laboratory, University of Minnesota. The facility is designed to allow precise control over basin floor subsidence, basin water level, sediment supply and water supply (Paola et al., 2001). The data used here are from the XES10 experiment, which occupied a 5.72-m-long, 2.98-m-wide portion (108 subsidence cells) of the XES basin and featured back-tilted subsidence and eustatic sea-level cycles . Subsidence rates were highest 1.2 m downstream of the sedimentwater inlet, where they reached ~3 mm/h, and decreased linearly basinwards at a rate of ~0.6 mm/h per meter downstream. Water level in the basin ("sea level") underwent eight cycles: one slow cycle, one rapid cycle, and one slow cycle with six superimposed rapid cycles (Figure 1). The sediment mixture used in this experiment was 63% quartz sand (110 μm), 27% coal sand (bimodal at 460 μm and 190 μm) and 10% kaolinite. Water supply was held constant at Q w = 25.0 kg/min throughout the experiment. Sediment supply during the first 132 h of the experiment was held at Q s = 0.58 kg/min and was subsequently decreased to Q s = 0.50 kg/min for the remainder of the experiment. Cyan dye and titanium dioxide were added to the water to enhance visibility. Overhead photos were captured every 10 s and high-resolution topographic scans were collected approximately every 4 h during the 310-h experiment.

| Extraction of river networks
The water in the XES10 experiment was dyed a cyan colour, which helps to distinguish water from dry land. A simple way to extract the distributary river networks from overhead photographs is to select a threshold of the red colour band in RGB space because cyan is the complement of red in an RGB colour model (Figure 2d,e). However, a simple threshold of the red colour band leads to the incorrect classification of the dark coal sand (a proxy for fine-grained sediment) as water ( Figure 2). Instead, we use principal components analysis (PCA) to extract river networks from XES10 imagery (Figure 2a-c). PCA is a ubiquitous tool for extracting meaningful correlations in large, multi-dimensional data sets (e.g. Jackson, 1991;Jensen, 2005;Jolliffe, 2002). PCA has also been successfully applied to hyperspectral imagery and 3-band aerial photography for detection of landscape characteristics that are otherwise difficult to identify remotely, such as soil salinity (e.g. Csillag et al., 1993;Metternicht & Zinck, 2003), and mapping and monitoring large woody debris in river basins (e.g. Leckie et al., 2005;Marcus et al., 2003;Smikrud & Prakash, 2006). In the XES10 overhead imagery, the red, green and blue colour bands are highly correlated to one another (e.g. the red band value decreases, and the blue band value increases in river channels), and PCA allows us to transform the RGB data into new variables that represent the maximum variances in the data set. We extracted the PCA coefficients from a representative example image using the open-source python function sklearn.decomposition.PCA (Pedregosa et al., 2011) and used the coefficients to transform all RGB images into PCA images. We find that variation in the first principal component (PC1; Figure 2a) most strongly corresponds to changes in the dry land surface, the second principal component (PC2; Figure 2b) most strongly corresponds to changes in wet vs. dry land and the third principal component (PC3; Figure 2c) contains limited information regarding surface characteristics. Importantly, the variance in PC2 is high for water but not for coal, making it a more reliable tool for river extraction than a simple red-band threshold ( Figure 2). We therefore use a threshold PC2 variance value to extract river networks for all overhead images in the XES10 experiment. PC2 values range from 0 to 255 and the threshold used in this study to distinguish water from land was 90. This threshold was determined using a histogram of all PC2 values, which shows a division between wet and dry land ( Figure S1), and by visually examining the river edge in images and finding the PC2 pixel value that corresponded to the transition from wet to dry. A range of PC2 values span the transition from wet to dry land and the choice of threshold is somewhat subjective; however, by comparing the histograms with maps of PC2 values and multiple photos ( Figure S1), a PC2 threshold can be set that represents a reasonable channel bank value. Additionally, the channel banks typically represent a large jump in PC2 values ( Figure S1), so small changes in threshold values are unlikely to significantly change the widths or distribution of channels. An example of the PCA method used in this study is available in the data repository (see data availability statement). A complementary Matlab example is available by Chadwick et al. (2021).

| Spatial alignment of photoscan pairs
Our strategy was to combine the imagery with the much more limited topographic data set to constrain channel location and local depth for each image. Overhead photos F I G U R E 1 Ocean-level variations through the course of the XES10 experiment. The ocean level is recorded as mm below the basin rim. Our analysis focuses on phases one through seven (P1-P7) because subsequent phases were so short lived that they deposited relatively thin stratigraphic sections were captured every 10 s over the 310 h of the XES10 experimental runtime, corresponding to a total of 111,600 photos. A total of 124 of the overhead photos could be paired with contemporaneous topographic scans ( Figure 3) and were georeferenced with their scan using QGIS 3.4 (QGIS.org, 2021), allowing for comparisons between surface characteristics apparent in the photos and surface elevation from topographic scans. PC2 images were generated for all photos and imported into QGIS, creating a stack comprised of a topographic scan, a colour photo and a PC2 image (Figure 3). Cross-section lines were defined across a wide range of distributary channels, and information was compiled from all data sources, starting with channel depth, which was compared with channel width, pixel RGB values and PC2 variance ( Figure 3). Here, channel width is measured by water level in photos, and in this case is generally equivalent to bankfull width. Note that all listed parameters can be easily extracted from imagery for comparison with channel depth, which requires a topographic scan. We aim to define a relationship between channel depth and characteristics measurable from imagery in order to constrain channel depth in the many photos that do not have correlative topographic scans. Bathymetric profiles from channels spanning a range of sizes, downstream distances and base-level conditions were measured to encompass the variability in channel depth within the experiments. This compilation eliminates potential trends in channel depth related to proximity to the delta apex, base-level conditions, or subsidence patterns, but simplifies the amount of information needed to generate synthetic stratigraphy.

| Constraining channel depth
The channels in the XES10 data set have a mean depth of 6.57 mm and an average width-to-depth ratio of 27.5. In natural systems, width-to-depth ratios span a wide range of values, from 5 to 200 (Larkin & Sharp, 1992;Schumm, 1963). The width-to-depth ratios of the XES10 experiment falls within typical ranges reported for threads of braided gravel-bed rivers in natural systems (e.g. Metivier et al., 2016), which is expected because the XES10 experiment has bedload-dominated channels with low bank stability (Schumm, 1963).
Channel centrelines and widths can be calculated from photos relatively easily, but channel depth cannot be precisely estimated from photos. After aligning photos with topographic scans in QGIS, we measured depth profiles across 170 channels in the XES10 data set. A histogram and a cumulative distribution frequency (CDF) curve of channel depths were built from the measured channel bathymetry ( Figure 4). As mentioned above, bathymetric profiles from channels spanning the full spectrum of conditions were measured to encompass the variability in channel depth within the experiments. There may be trends between channel depth and downstream distance or base-level conditions; however, we expect associated changes in channel width or pixel intensity to reflect these F I G U R E 3 Topographic scan paired with colour photo. We measured channel bathymetry, red-band intensity and PC2 intensity for channels throughout the experiment to build a database of channel depth, width, red-band intensity and PC2 intensity for our data set. variations. Values of channel depth, pixel colour and PC2 variance were averaged across a channel cross section and compared to channel width to provide broad constraints on the relationships between channel depth and width, pixel colour or PC2 value. For example, average red-band values are commonly associated with shallow channels ( Figure S2). We focused on channel width, rather than PC2 value or pixel colour because it is the easiest to measure and incorporate into the workflow. Precise depth measurement from image data is generally not possible; instead, our goal is to make reasonable depth estimates, constrain the inevitable uncertainty in these estimates, and explore how uncertainty affects connectivity results. We considered two potential methods for estimating channel depth. There is a continuous relationship between channel width and depth (SI 2) and, in the first method, channel depth could be estimated using the linear regression between measured channel widths and depths. In the second method, the channels can be binned based on their relationship between channel width and depth, and histograms of channel depths can be built within those width categories. Channel depth can then be estimated using conditional probability, that is, by randomly drawing from the distribution of depths measured for a channel in a given width category. The trade-off between these two methods -width as a continuous predictor vs. conditional probability-is a balance between representing the residual (non-systematic) variance versus representing the trend in the data set. Linear regression prioritizes the trend in our data and explains nearly one-third of the variability in the data set (r 2 = .28; p < .001), whereas the conditional probability method prioritizes the residual variance in the data while reducing the trend to a set of points equal to the number of bins created. The data also show heteroscedasticity (the variability varies), which further supports the use of binned data. The more bins used in the conditional probability method, the more closely this method follows the linear regression. However, the number of bins is limited by the number of data points available. On balance, we choose to use conditional probability because we think it better captures the main variability in the data set, but either method is workable.
We categorize our channels as 'wide', 'medium', and 'narrow' and build histograms of channel depth for each of these categories (Figure 4). Wide channels are classified as those wider than 23 cm and narrow channels are those narrower than 12 cm. These width ranges were selected by looking at width histograms ( Figure 4a) and by searching for width categories that provided distinct depth distributions without significantly reducing the number of data points in each category (Figure 4b). We used conventional Monte Carlo methods to sample our depth distributions, which confirmed that the 'wide', 'medium' and 'narrow' distributions differ significantly enough from the bulk channel depth distribution to justify categorizing our channels in this way. We start with this basic channel depth estimate and evaluate the sensitivity of subsurface architecture to this assumption. We used two methods to test the sensitivity of the synthetic stratigraphy to channel depth estimate: a) we ran 100 iterations of the synthetic stratigraphy code at 16 different locations in the delta, to generate 100 cross sections with distinct depth profiles for all channels and b) we generated a modified histogram of channel depths ( Figure 5), with a similar skewness, but a broader range of channel depths to draw from and again ran 100 iterations of the synthetic stratigraphy code with this exaggerated depth distribution. This modified depth distribution does not reflect the channel depths measured in the experimental data set but was instead used to estimate how sensitive the stratigraphy might be when there is a high likelihood of selecting an unreasonable channel depth. The measured channel depths range from 0.61 to 23.19 mm, with a median and mean of 6.08 and 6.57 mm respectively, and a positively skewed distribution. The modified channel depths range from 1.9 to 35.8 mm, with a median and mean of 12.4 and 13.2 mm, respectively, and a positively skewed distribution. The modified depths and modified depth distribution (magenta dots). This is for stratigraphy generated during a slow sea-level fall (phase two) at the location X-X′. Stratigraphy for this location and phase is shown in (e). (d) Connectivity values for 100 iterations of synthetic stratigraphy using measured channel depth distribution (black dots) and modified depth distribution (magenta dots). This is for stratigraphy generated during a slow sea-level rise (phase three) at the location X-X′. Stratigraphy for this location and phase is shown in (f). (e) One iteration of synthetic stratigraphy generated at X-X′ during the slow sea-level fall (phase two). White represents channel fill and black represents overbank sediment. This synthetic stratigraphy was generated using the measured depth distribution, but as can be seen in (c), the stratigraphy is not sensitive to using measured vs. modified depth distributions. (f) One iteration of synthetic stratigraphy generated at X-X′ during the slow sea-level rise (phase three). This synthetic stratigraphy was generated using the measured depth distribution, but as can be seen in (d), the stratigraphy is not sensitive to using measured vs. modified depth distributions Connectivity were randomly generated using the scipy.stats.skewnorm function in Python 3.8 (Virtanen et al., 2020).

| Building a synthetic cross section
The generalized workflow for building stratigraphy from a succession of overhead photos involves: (a) defining a cross-section line, (b) identifying channel locations across that line for a single image, (c) estimating an approximate topographic profile for the given image by assigning depths to all channels, (d) subsiding the profile based on subsidence and sedimentation rates and (e) repeating the process for each successive image, vertically stacking individual profiles and allowing channels to 'erode' underlying stratigraphy where current elevation is below previous values. This process is described in more detail below. The synthetic stratigraphy code and the XES10 data used in this study are publicly available in the data repository (see data availability statement).
Defining channel boundaries from overhead photos is not straightforward, particularly in experiments that are prone to sheet flow. Channel banks are frequently overtopped, making dye colour somewhat misleading for defining channelized flow. Visual comparison suggests that PC2 images most clearly distinguish channelized flow from overland flow and/or dry land; therefore, we use this component for automatically extracting channel widths from imagery (Figures 2 and 3). A PC2 threshold can be defined to identify 'wet' pixels, such that a Boolean variable can be generated to indicate whether each pixel along the cross section lies within a channel or outside of a channel. Each channel is assigned a sinusoidal, concave-up profile with a maximum depth drawn from the histogram of depths compiled from the photo-scan pairs and a halfwavelength corresponding to the channel width.
Each topographic profile for a given image contains 'overbank' with elevation value 0 and channels with negative elevations (Figure 6b). This is an oversimplification of a topographic profile; however, we do not expect the large-scale connectivity structure to be sensitive to mmscale variations in topography adjacent to channel banks. A delta with more cohesive banks may build more floodplain relief than observed in the XES10 experiments, which could potentially influence subsurface connectivity. In the application of this method to a less mobile and more cohesive system, it may be reasonable to give a more complex floodplain elevation profile, such as a decrease in elevation away from the channel margins. Aggradation in the XES10 experiment is a combination of local subsidence (i.e. subsiding basin floor) and base-level driven accommodation. Rates of strike-averaged surface aggradation or degradation throughout the experiment are added to basin-floor subsidence to create net sedimentation rates through time at any location in the delta. Subsidence is applied to each profile based on these sedimentation rates (Figure 6b). These simplified profiles, once many thousands of them are stacked together, provide a first-order estimate of the subsurface distribution and connectivity of channelized sands.
Synthetic stratigraphy generated using time-series of topographic scans has been shown to provide good approximations of true stratigraphy (Paola et al., 2001;Straub et al., 2009;Yu et al., 2017). The novelty of our proposed method lies in the generation of synthetic stratigraphy primarily from sequential photos of experiments paired with much more limited topographic data. Following the XES10 experiment, the deposit was sliced in dip-oriented sections F I G U R E 6 Example of topographic profiles generated from imagery. (a) Principal Component 2 (PC2) image used to distinguish active channels from floodplain. (b) A topographic profile generated along the line X-X′. Floodplain is given no elevation, and channels are given a depth drawn from a measured distribution ( Figure 7) and a concave-up profile. This profile is subsided to ~875 mm below the surface based on aggradation rates calculated for the experiment. (c) The final stratigraphic surface after successive time steps generate channels that may erode into pre-existing surfaces. This is similar to many stratigraphic surfaces preserved in natural systems, which are often composite, timetransgressive surfaces (along the length of the delta from the apex to the shoreline) and photographed. Composite, strike-oriented images were created by stitching together thin slices of the photos. We compare our synthetic stratigraphy to the composite images to see whether the synthetic stratigraphy captures the broad patterns recorded in the real stratigraphy from the XES10 experiment (Figure 7). The synthetic stratigraphy tracks only subaerial, topset stratigraphy, so any subaqueous deposits such as foreset or bottomset strata appear as black bands in the synthetic stratigraphy (Figure 7b,c).

| Measuring subsurface connectivity
We define a sandbody as a single connected component, that is, all pixels within that sandbody are connected to one another using 4-way connectivity in two dimensions. We use Connectivity, C, (Equation 1) as a simple measure of subsurface architecture, and define it as the ratio of the largest sandbody to the sum of all sandbodies in the cross section (Hovadik & Larue, 2007 F I G U R E 7 Comparison between synthetic stratigraphy and real stratigraphy generated during the XES10 experiments. In synthetic stratigraphy, white represents channel fill and black represents overbank sediment. The XES10 experiment was cut in a dip-orientation (i.e. from the apex to the toe of the delta), so high-resolution photos of the stratigraphy are oriented perpendicular to our synthetic stratigraphy. However, composite strike-oriented photos were created by compiling slices of dip-oriented photos. Missing information (black voids) in the photos is due to sediment failure during cutting. (a) Base-map showing the location of sections shown in (b, c). (b, c) Real stratigraphy (left) and synthetic stratigraphy (right) for the locations indicated in (a). Annotated versions of both real and synthetic stratigraphy are included, which highlight the locations of foreset and bottomset strata (grey polygons). An elongate barrier composed of overbank sediment is captured by the synthetic stratigraphy and highlighted with the blue arrow. The synthetic stratigraphy generally captures the largescale structure well, including the locations of delta foresets (grey polygons) and the scattered nature of the fine-grained sediments (black in photos and in synthetic stratigraphy). Some preserved fine-grained barriers that are not captured in the synthetic stratigraphy likely represent mud-filled channels; our method assumes that all channels will be sand-filled, which is not always the case (c) B. C.

(a) (b)
where A B is the area of the largest sandbody and ∑ A i is the sum of all sandbody areas. The area of the largest sandbody (A B ) is used over the mean sandbody area because the size of the largest connected sand package represents the longest potential flow pathway, that is, the 'connected sand fraction', which is an important parameter in reservoir modelling (Allen, 1978;Hovadik & Larue, 2007;King, 1990). The XES10 delta was a bedload dominated delta with no floodplain cohesion, meaning that channels were extremely mobile and produced highly connected stratigraphy. Because connectivity is so high in this system, we consider the possibility that the distribution of fine-grained barriers to flow may be more critical to flow pathways. Here, we consider any floodplain deposit as a potential barrier to flow and, assuming a barrier of any thickness has the potential to block or divert subsurface flow pathways, we use the barrier width as a measure of subsurface connectivity structure. We calculate the maximum width of each barrier (i.e. distance between the farthest left and right pixels), and the maximum barrier width (W) is the width of the widest barrier in a given cross section.
Each cross-section panel spans a range of sea-level conditions and, therefore, distances from the shoreline at any given point, so we break up the synthetic stratigraphy into sections based on sea-level conditions (Figure 1). We focus on the following sea-level phases: P1 -constant sealevel, P2 -slow fall, P3 -slow rise, P4 -constant sea-level and P6 -rapid rise. These are the phases that preserved the thickest stratigraphic sections. The P5 phase of rapid sea-level fall was net degradational and did not preserve basin-wide stratigraphy. (1)

| large-scale architecture
Synthetic stratigraphy, created here by stacking simple topographic profiles, highlights large-scale stratigraphic architecture in which frequently channelized regions are easily distinguished from those where channels rarely reworked the delta surface (Figure 8). Regions with little to no channelization build strata that are dominated by overbank sediment and appear devoid of channels in cross section. These overbank-dominated strata have extremely limited 2D vertical or lateral connectivity due to a lack of coarse-grained conduits.
Our synthetic stratigraphy does well at representing the general distribution of fine-grained barriers and amalgamated sand bodies in the topset deposits (Figure 7). Slight discrepancies between the synthetic stratigraphy and the real stratigraphy exist because our method cannot account for mud-filled channels, which form when active channels are abandoned and filled in by fine-grained floodplain sediment. However, the broader floodplain deposits that are preserved because of floodplain aggradation without extensive reworking or channelization are well characterized by the synthetic stratigraphy (Figure 7; blue arrow). Any subaqueous deposits, that is, foresets and bottomsets, appear as black bands in the synthetic stratigraphy (highlighted by grey polygons in Figure 7). When compared to the composite imagery of the real experimental stratigraphy, these black bands align well with the foreset and bottomset stratigraphy (Figure 7b,c). Generally, our method of building stratigraphy using only surface information can accurately characterize large-scale subsurface architecture (Figure 7). By breaking the stratigraphy into separate sea-level phases prior to calculating connectivity metrics, we prevent foreset strata -which are simplified here to be entirely fine-grained deposits -from affecting measures of connectivity.
Overall, we find that large-scale subsurface architecture is so dominated by channel kinematics coupled to sedimentation that it is relatively insensitive to the uncertainty in estimates of channel depth (Figure 9). At basin scales, subsurface connectivity structure is controlled primarily by channel location, which in this case is measured via time-lapse photography (Figures 7 and 8). This is highlighted by the three iterations of synthetic stratigraphy shown in Figure 9. In this example, the three iterations are qualitatively very similar, and the large-scale architecture is relatively unaffected by the variations in channel depth estimates. The subtle differences among these three panels are primarily in the shape of the fine-grained barriers. However, these variations do not change the 2D connectivity structure. We explore these comparisons quantitatively in the following sections.

F I G U R E 9
Comparison between multiple iterations of the same 2-dimensional section. (a) Base-map showing the location of synthetic stratigraphy. (b-d) Three different iterations of synthetic stratigraphy generated in the same location. Each iteration re-draws the channel depth, so the only differences between these will result from slightly different depth assignments to channels. The large-scale structure of these sections is clearly not sensitive to uncertainties associated with estimating channel depth

| Sandbody connectivity
Overall, the connectivity in these stratigraphic sections is extremely high (Figure 10). The majority of sea-level phases and locations have sandbodies that are over 90% connected (C > 0.90). A comparison of connectivity values across 100 iterations of the same cross section (i.e. only the channel depth draw varies) shows very little variation, confirming that the depth estimates are not significantly affecting large-scale sandbody connectivity ( Figure  10). Furthermore, comparison between the 100 iterations using the measured depth CDF compared to the expanded depth CDF shows little difference between connectivity values ( Figure 5). The expanded depths do provide slightly higher connectivity measurements, 0.95 versus 0.98; however, these slightly higher values are not expected to result in a great difference in fluid flow. An error in the sediment feeder during the experiment caused an increase in sediment supply at the boundary between Phase 3 and Phase 4 (132 h into the experiment).
This change in sediment feed rate resulted in a water-tosediment-supply ratio of 43:1 during sea-level phases 1-3 and 50:1 during all subsequent phases. There is no discernible change in sand connectivity associated with the slight change in sediment supply ( Figure 10). Box plots of connectivity values show a subtle increase in connectivity downstream; however, again, this increase is fairly small overall ( Figure 10). These changes in connectivity are subtle because of the high sand content in this system, and we expect that trends might be similar, but more noticeable, in a more mud-prone system. Slightly lower connectivity values in the proximal sections are most likely due to lower channel mobility and fewer distributary channels near the delta apex, leading to a higher potential for overbank sediment preservation. Connectivity is more variable during phases of rapid sea-level rise (Phase 6 and Phase 9; Figure 10e,g). Not only does connectivity show variability between locations, but it also appears more sensitive to channel depth estimates; for a given location, the stratigraphy deposited during rapid sea-level rise has a broader F I G U R E 1 0 Connectivity boxplots for sand bodies in synthetic stratigraphy. Each box-and-whisker set represents 100 iterations of the synthetic stratigraphy for a given location and sea-level phase (see Figure 9 for an example of 3 iterations). Boxes extend from quartile one to quartile 3 of the data, line is at the median, and whiskers extend across the range. Circles plotting outside the whisker range are considered outliers. The very high connectivity values in all locations suggest that most sand in this experiment is essentially fully connected. The spread of the box-and-whiskers indicates the variations in connectivity due to the estimation of channel depth. Overall, there is very little variation due to channel depth estimates. Lower connectivity values and/or more variation in connectivity between iterations (spread of boxes and whiskers) appear during phases or rapid sea-level rise (e & g) and near the zones of maximum subsidence (~100 cm downstream; a, e-g). High subsidence rates and high rates of sea-level rise are associated with high sediment aggradation, which may result in less sand connectivity  (Figure 10), which we interpret as a reflection of greater aggradation rates in the delta plain (i.e. topset). Because we focused our analysis on strike-oriented cross sections, and because river networks often have high downstream connectivity, the 3-dimensional connectivity is likely to be even higher. Early fluvial models linked greater riverplain aggradation rates to reduced channel connectivity (Allen, 1978;Bridge & Leeder, 1979;Leeder, 1978) and promoted the importance of fluvial scour in generating connections to underlying sandbodies (Salter, 1993). The balance between the conduits for fluid flow (channels) and barriers to flow (overbank fines) is set in large part by the ability of channels to erode and rework overbank deposits and provide a connection to previously deposited sandbodies (Salter, 1993). During high rates of aggradation, floodplain deposits may become thick relative to channel scour depths, resulting in a lower likelihood of vertical connection and therefore more sensitivity to variations in channel depth ( Figure 11). These findings agree with recent numerical modelling results that find lower sandbody connectivity under rapid rates of sea-level rise .

| Barriers to fluid flow
In sandy and highly connected systems like the XES10 experiment, sandbody connectivity metrics are all very high and may provide relatively little insight into subsurface structure. In these cases, characterizing the distribution and dimensions of the fine-grained overbank deposits, that is, barriers to flow, becomes an important constraint on subsurface fluid flow. We suggest that in sand-dominated stratigraphy, narrow barriers are unlikely to change subsurface flow patterns, but laterally extensive barriers may hinder fluid flow or at least greatly increase the length and tortuosity of flow pathways. For each synthetic stratigraphic section, we calculated the maximum barrier width, W. Each cross section spans the entire width of the delta, so we divided W by the width of the cross sections to give a normalized width W* (Figure 12). Most barriers span <30% of the crosssection width (W* < 0.3), and the period of slow sea-level rise has slightly wider barriers than the constant sealevel phases or the slow sea-level fall (Figure 12). Barrier width decreases downstream, which is likely a reflection of channel bifurcation leading to a greater number of distributary channels downstream. An example of wide vs. narrow barriers is shown in Figure 13. The proximal section (Figure 13c) contains a large barrier that spans across roughly 30% of the entire width, while the distal section (Figure 13d) contains more barriers that are much smaller. These trends are discussed in more detail in the following section.

| DISCUSSION
Our finding that large-scale connectivity structure is far more sensitive to the history of channel location, coupled to net sedimentation, than to channel bathymetry is encouraging because it suggests that we are not limited by the variable that is typically the most difficult to measure or predict. A general understanding of the expected range of channel depths appears sufficient if paired with a knowledge of channel kinematics and network evolution. Previous work by Salter (1993) highlights the importance of understanding fluvial processes -such as channel type, scour, incision -in controlling dimensions of preserved sandbodies and cautions against the use of simple thickness-to-width ratios in fluvial reservoirs. Our results support these claims that kinematics play a critical role in subsurface connectivity structure (Heller & Paola, 1996;Sahoo et al., 2020;Salter, 1993). The methods we describe can readily be applied to other experimental data sets, particularly where stratigraphic sections were not obtained. More importantly, the basic situation in the experiments applies to the field as well: it is much easier, and much more common, to obtain records of channel location and width than channel bathymetry. With the advent of historical imagery engines such as that in Google Earth, such records are now easily accessed on a global scale. The important implication of this study for field data is that the channel kinematics dominates overall stratigraphic architecture and connectivity, so it is feasible and worthwhile to estimate stratal geometry from imagery alone, even using only crude depth estimates to provide the vertical dimension.
F I G U R E 1 1 Schematic diagram illustrating low aggradation rates (a) resulting in thinner floodplain packages which allow channels do incise into underlying sandbodies. High aggradation rates (b) make these connections less likely, leading to sandbody connectivity that is more sensitive to channel depth estimates

| Channel kinematics and sandbody connectivity
Here, we used an experimental data set that includes imagery with extremely high temporal resolution, which effectively allows us to constrain kinematics with high precision by measuring channel locations through time. A natural next step for this work is to link it to predictive understanding of channel kinematics. For the purposes of this discussion, we refer to channel mobility as any channel movement, which can include migration (continuous lateral movement that reworks sediment in its path) and avulsion (abrupt shift in channel location that leaves un-reworked sediment between the previous and new channel paths). The full range of controls on channel mobility and the interplay of migration and avulsion are still incompletely understood, but recent research has shown how channel migration is influenced by sediment and water supplies, and incision history (e.g. Bufe et al., 2016Bufe et al., , 2019Constantine et al., 2014;Wickert et al., 2013).
Experimental and field studies suggest that channel mobility increases with increased sediment supply and discharge (Dunne et al., 2010;Wickert et al., 2013). These changes in channel mobility are intimately connected to sea-level cycles: rising sea-level drives floodplain aggradation, which reduces downstream sediment flux (Bufe et al., 2019). Recent analysis of XES10 channel kinematics using new Particle Image Velocimetry techniques showed that the 15% reduction in sediment supply 132 h into the experiment (starting in Phase 4) led to lower migration rates (Chadwick et al., 2021). This decrease in lateral migration is in contrast to avulsion frequency, which tends to increase during sea-level rise (Chadwick et al., 2020;Liang et al., 2016). In XES10, Phase 4 experienced constant sea-level with lower sediment supply and slower channel migration compared to the previous constant sea-level Phase 1 (Figure 10; Chadwick et al., 2021). With the exception of a very small reduction in connectivity in the most proximal reaches of the delta (<100 cm downstream), the reduced sediment supply does not appear to F I G U R E 1 2 Boxplots for the maximum barrier width (W*) in synthetic stratigraphy, where W* is calculated as the width of the widest barrier divided by the width of the entire cross section. Each box-and-whisker set represents 100 iterations of the synthetic stratigraphy for a given location and sea-level phase (see Figure 9 for an example of 3 iterations). Boxes extend from quartile one to quartile 3 of the data, line is at the median and whiskers extend across the range. Circles plotting outside the whisker range are considered outliers. The lateral extents of all barriers in a given section were measured and the largest value is considered the maximum barrier width. Barrier width generally decreases downstream in all sea-level phases influence sand connectivity (Figure 10a,d). However, the combined influence of rapid sea-level rise and decreased sediment supply resulted in more variability in sandbody connectivity (Figure 10e,g). Classic fluvial models show that in 2D stratigraphic panels, channel connectivity should be lower during high rates of floodplain aggradation (Allen, 1978;Bridge & Leeder, 1979;Leeder, 1978). Aggrading systems drive avulsions through channel backfilling and superelevation (Mohrig et al., 2000), which can lead to reworking of fresh floodplain surfaces (Wickert et al., 2013). Rising sea-level drives high rates of fluvial aggradation, and we find that the synthetic stratigraphy generated during rapid sea-level rise shows more variability in sand connectivity ( Figure 10).
We propose that increased variability and sensitivity to channel depth estimates during rapid sea-level rise and in regions of high subsidence are tied to more frequent avulsions and reworking of fresh topset surfaces in the XES10 experiment ( Figure 10). Channels that are reworking recently occupied surfaces are likely to incise into older channel sands, whereas the ability for avulsed channels to form vertical connections with older channel sands relies on their ability to incise through fine-grained floodplain sediment. Furthermore, if channel mobility is driven by avulsion rather than lateral migration, there is an increased potential to preserve isolated floodplain barriers in the delta plain. This variability may become less apparent in a 3-dimensional analysis because channels are more highly connected in the downdip direction.

| Fine-grained barriers
The sandy and highly connected nature of the deposits developed in our experimental case would be highly conductive to groundwater flow at field scales. In cases like this, the key to understanding aquifer properties may instead be to characterize the barriers in the system. In sand-dominated stratigraphy, subsurface barriers add uncertainty to flow and transport behaviour. Scattered, small and generally disconnected barriers are unlikely to substantially affect overall fluid flow in a three-dimensional sandy subsurface, but they can be important for transport by increasing effective dispersion. More extensive barriers perpendicular to the hydraulic gradient can increase tortuosity of flow, even in connected systems, reducing effective permeability in the direction of the gradient and increasing the spread of solutes in the direction perpendicular to the gradient. Thus, understanding when and where barriers become more connected and extensive ( Figure 13) is important for predicting subsurface fluid flow and contaminant transport. We chose the maximum width of barriers as a measurement of their maximum F I G U R E 1 3 Examples of high and low barrier width. Both examples are from Phase 1, which was a constant sea level that allowed slow progradation of the delta. (a) Boxplot of maximum barrier width (W*) during Phase 1 (see Figure 12 for further explanation of the plots). Locations of sections shown in (b, c) are indicated by the arrows. (b) Boxplot of sand fraction. Although W* decreases between 125 and 200 cm downstream, the sand fraction remains the same, implying that the decrease in maximum barrier width downstream is not due to increased sand content, but instead due to changes in the geometry of barriers. (c) Synthetic stratigraphy generated during Phase 1 at a location 125 cm downstream of the delta apex. This section has relatively wide barriers (~20%), which can be observed by the large connected barrier in the top right region of the stratigraphy. (d) Synthetic stratigraphy generated during the same phase but at a more downstream location has a narrower maximum barrier width (~9%). This narrower value can be seen as many small barriers scattered throughout the stratigraphy extent perpendicular to vertical flow transport. The barriers in the XES10 synthetic stratigraphy are generally elongate and wider than they are thick. It is possible that a more cohesive delta may have a barrier that cuts diagonally across stratigraphy; however, even in this case, the maximum width that the barrier spans still provides an estimate for the extent of flow diversion. Other metrics for barrier characterization could include average barrier width, orientation or barrier connectivity. Although there may be more variability in sandbody connectivity during rapid sea-level rise (phases 6 and 9), the scattered and disconnected nature of the fine-grained barriers under these conditions should produce more consistency in barrier dimensions, that is, less uncertainty. In these cases, the increased avulsion frequency and greater reworking of topset surfaces is likely to dissect fine-grained floodplain deposits, resulting in more narrow barriers and decreasing the likelihood of a large, laterally extensive barrier. Furthermore, the downstream decrease in barrier width seen in all sea-level phases indicates that overbank sediment becomes more dissected downstream (Figures 10 and 12). Interestingly, the ratio of channel fill to overbank material (net to gross) does not vary appreciably downstream, with the exception of slightly lower values in the most proximal reaches of the delta attributed to the fixed inlet (Figure 13b). This implies that the relative proportions of overbank material remain similar downstream, but the organization and size of individual overbank packages decrease. Channel kinematics therefore affects not only the relative proportions of channel fill to overbank deposits but also the geometries and preservation extent of fine-grained, overbank strata.

| Applicability to cohesive systems
Results from this study are most applicable to bedloaddominated systems with poorly developed floodplains. The XES10 experiment did not have cohesive banks, meaning that our floodplains are less analogous to vegetated lowland systems that resist channel encroachment and slow down lateral migration, and are instead an end-member example of coarse-grained channel systems with highly mobile banks. The subtle differences in connectivity and large-scale structure between different locations and different depth estimates (Figures 9 and  10) are best explained by the frequent reoccupation of channels and amalgamation of channel bodies. Potential errors in channel depth estimates may therefore be hidden in the amalgamated channel bodies, making them inconsequential to simple connectivity metrics. Deltas with more cohesive sediment or with vegetated banks have lower rates of lateral migration, greater floodplain relief and more single-thread channels (e.g. Caldwell & Edmonds, 2014;Edmonds & Slingerland, 2010;Hoyal & Sheets, 2009;Lauzon & Murray, 2018). A clear avenue for future research is to test the applicability of this method to less mobile or less sandy systems. We hypothesize that lateral connectivity will decrease due to less mobile channels; however, this does not preclude channel mobility as a key factor in subsurface sandbody connectivity. A potential modification to the code presented here would be to introduce floodplain topography or a greater degree of relief between channels and floodplains, which may change channel stacking. A comparison between floodplain relief from topographic scans and channel bathymetry may highlight the influence of our simplified topographic profiles in more cohesive systems.

| Field-scale analogues
Field-scale analogues may be found in the ancient rock record prior to the advent of land plants, or in foreland basin systems with rapidly uplifting source areas such as the Ganges-Brahmaputra-Meghna (GBM) basin. Holocene stratigraphy in the GBM basin is sand-dominated, with very little preservation of overbank fines due to the highly mobile nature of the river networks (Goodbred et al., 2003;Pickering et al., 2014;Umitsu, 1993). Understanding the distribution of fine-grained barriers and the connectivity of sandstone aquifers in the GBM basin is a major factor in predicting large-scale groundwater flow, and transport of contaminants, in this heavily populated region (Khan et al., 2019;Michael & Voss, 2009).
The GBM basin has further similarities to our experimental data set, in that our fixed delta apex, which resulted in lower sandbody connectivity and larger W*, is analogous to network locations that have remained relatively stable. In the GBM Delta, the confluence of the Ganges and Jamuna (Brahmaputra) rivers forms an avulsion node (Reitz et al., 2013;Wilson & Goodbred, 2015). Although major river confluences can be mobile (Dixon et al., 2018), these nodes serve as regions of relatively lower lateral channel mobility and are analogous to conditions seen in the most proximal regions of the XES10 delta. Our results highlight that channel connectivity may be lower ( Figure  10) and barrier width may be greater ( Figure 12) near these confluences, which has important implications for reservoir and aquifer extent in large distributary systems.
The complex relationships among channel mobility, allogenic and autogenic controls and the resulting subsurface structure remain an exciting avenue for continued research. Experimental data sets such as the XES10 experiment remain a useful tool for distilling these complexities into stratigraphic predictions. Our finding that connectivity is controlled primarily by channel kinematics indicates that future research should focus on quantifying and predicting kinematics (e.g. Chadwick et al., 2021;Jarriel et al., 2021) to better understand the key control on subsurface architecture and connection.

| CONCLUSIONS
Based on combining an experimental data set with highfrequency overhead imagery and much sparser measurements of channel topography, we find that: • Channel location as determined from overhead imagery is much more important for basin-scale subsurface channel architecture and connectivity than is the channel depth. This is useful because depth is typically difficult to constrain at high spatial and temporal resolution, both in experiments and in the field. Results highlight the importance of predicting channel kinematics to fill in the gaps in data sets with lower temporal resolution. • Estimates of bulk stratigraphic connectivity of channel sand bodies are only slightly affected by even large changes in the probability distribution of channel depth • In highly mobile, sandy systems (sensu Jerolmack & Mohrig, 2007), connectivity tends to be high regardless of external controls. Connectivity is most sensitive to depth estimates during periods of rapid relative sea-level rise, which we hypothesize to be driven by increased floodplain aggradation making vertical connections to underlying sandbodies more tenuous (Bridge & Leeder, 1979). • The expected main influence on subsurface flow in highly connected systems is preservation of overbank deposits. At basin scales, the distribution and extent of these barriers to subsurface flow are controlled by channel mobility and the relative influence of lateral migration versus channel avulsion. Although the ratio of channel fill to overbank deposits varies little downstream, overbank deposits are more likely to be wider and therefore more effective barriers to subsurface flow in the more proximal reaches of the delta.