Reconstructing backwater hydrodynamics from fluvial-deltaic deposits using stratigraphic inversion: An example from the Tullig Sandstone, Western Irish Namurian Basin, County Clare, Ireland

Fluvial-deltaic stratigraphy archives information used to reconstruct quantitatively paleohydraulic conditions. By using these reconstructions, one can gain insight into how both internal dynamics and changing environmental conditions affected surface processes, lowland landscapes, and sediment transport of ancient systems over time (Blum and Törnqvist, 2000; Hajek and Heller, 2012; Foreman et al., 2012; Lynds et al., 2014; Lyster et al., 2021; Barefoot et al., 2022). The accuracy of paleohydraulic reconstruction methods is limited, however, as natural systems demonstrate variable flow conditions that are not easily discernable based on evidence from the stratigraphic record (Leary and Ganti, 2020). Moreover, many paleohydraulic relations rely on the assumption of steady and uniform flow (Lynds et al., 2014; Mahon and McElroy, 2018), which is not suitable for fluvial-deltaic deposits influenced by non-uniform (backwater) flow. As rivers approach their terminus, the water surface and channel bed slopes diverge, giving rise to channel deepening and decreasing flow velocity in the down-dip direction. This divergence lowers sediment transport capacity (Lane, 1957; Lamb et al., 2012; Nittrouer, 2013) and triggers morphodynamic feedbacks (Nittrouer et al., 2012; Ganti et al., 2014, 2016), in particular, bed sediment aggradation and enhanced downstream fining (Smith et al., 2020), as well as adjustments in bedform morphology (Wu et al., 2020b). A first-order approximation of the upstream extent of backwater influence (L_bw), measured upstream of the river outlet, is estimated as the ratio of bankfull flow depth (H_bf) and water surface slope for the normal (steady and uniform) flow region (S_w): L_bw = H_bf/S_w (Paola and Mohrig, 1996). Water surface slope S_w and channel-bed slope S are subequal at reach scale in the condition of uniform flow so that they are interchangeable when calculating the backwater length L_bw (Paola, 2000). Non-uniform flow can extend from tens to hundreds of kilometers for low-sloping rivers. The large extent of non-uniform flow can cause changes in depositional patterns and impact stratigraphic architecture. However, the ability to assess the paleohydrology of fluvial-deltaic systems from the rock record is limited, due to the lack of tools for reconstructing first-order channel-bed slope and bankfull flow depth for non-uniform flow.

Stratigraphic signatures of backwater hydrodynamics have been inferred from several recent rock-outcrop studies. Petter (2010) recognized a down-dip fining in median channel-bed grain size and increase in flow depth estimated from an increase in fluvial bar thickness in the Castlegate Sandstone (Utah, USA). Fernandes et al. (2016) showed similar spatial trends in bar thicknesses of Holocene-aged Mississippi River deposits, coinciding with a downstream narrowing of the channel belt. Similar observations regarding channel-belt width were made in the Ferron Sandstone of Utah (Kimmerle and Bhattacharya, 2018) and in the subsurface of the Mungaroo Formation of Australia (Martin et al., 2018). Spatial variability in the dimensions of autogenic scours (Ganti et al., 2019) was documented in the Castlegate Sandstone and is attributed to non-uniform flow conditions (Colombera et al., 2016; Trower et al., 2018). These studies all assess rock-outcrop strata and link observed trends in fluvial-deltaic strata to hydrodynamic conditions evolving over short time scales (10⁰–10² yr), assuming that the extent of the backwater zone is relatively static.

An important and rarely investigated aspect of the linkages between stratigraphy and backwater hydrodynamics is to understand how backwater conditions and channel morphodynamics are preserved in the rock record for intermediate timescales (10²–10⁴ yr). While modern fluvial-deltaic systems bear the hallmarks of changing hydrodynamic conditions, they offer only a snapshot in time in the context of, for example, a delta lobe cycle. Experimental studies, on the other hand, effectively expand time and provide the opportunity to examine backwater-mediated stratigraphy constructed over temporal scales that include multiple channel avulsion cycles (Ganti et al., 2016; Chadwick et al., 2020); however, experiments lack the ability to resolve changes to strata that are observable at the scale of rock outcrops (e.g., extent/thickness of sandstone beds, and/or dune-scale cross sets). Numerical models have been used to predict how backwater sedimentation varies over centuries to millennia due to transgression/regression, lobe progradation, and shoreline movement (Parker et al., 2008b; Moran et al., 2017; Moodie et al., 2019); however, the transferability of findings from these studies to the rock record is incompletely understood. For example, the common practice of estimating paleo-channel depth by using bar and/or cross-set thickness measurements from channel stories is problematic if channel depth evolves in time and space, as is expected of rivers affected by backwater hydrodynamics and variable shoreline conditions driven by transgression and regression.

To address (1) the lack of paleohydraulic reconstruction methods for non-uniform flow, and (2) the impact of evolving backwater condition on stratigraphic architecture for intermediate time scales (10²–10⁴ yr), this study develops a novel stratigraphic inversion method that assesses how the transition between normal and backwater flow conditions evolves and impacts stratigraphic development. The inversion method provides estimates of channel-bed slope and bankfull flow depth using first-order geometric properties of fluvial-deltaic stratigraphy. Additionally, the outcomes are compared to measurements of channel-bed grain size, dune-scale cross-set thickness, and channel sandstone body dimensions, collected systematically across a series of outcrop exposures of a system spanning the fluvial to deltaic transition. The inversion method is used to estimate the extent of influence of non-uniform flow and thus provides a basis for choosing an appropriate location to apply paleohydraulic models intended for steady and uniform flow conditions.

Compared to standard paleohydraulic models, the inversion method provides comparable outputs of channel-bed slope and channel depth. Moreover, field data presented in this study document spatial variabilities in channel-bed aggradation, which are interpreted to arise due to backwater conditions. The magnitudes and locations of these measured changes agree well with the results from the inversion model. The analyses presented in this study support the use of stratigraphic features observed in outcrops to interpret backwater conditions of ancient fluvial-deltaic systems, and furthermore, show how depositional features can be used to infer paleohydraulic conditions in general.

The present study collected field measurements from the Tullig Sandstone of the Western Irish Namurian Basin (Co. Clare, Ireland, Fig. 1A; Stirling, 2003; Best and Wignall, 2016). The Western Irish Namurian Basin developed from the Late Devonian (395 Ma) to the early Carboniferous (317.5 Ma) in association with extensional reactivation of pre-existing structural lineaments associated with the Iapetus suture, the location of which coincides with the present-day Shannon Estuary (Fig. 1B; Dolan, 1984; Besly and Kelling, 1988; Haszeldine, 1988; Nenna and Aydin, 2011; see Supplemental Material¹ for additional details). Stratigraphic relationships, and dominantly NE paleo-transport indicators, suggest that the Western Irish Namurian Basin filled progressively from south to north (Wignall and Best, 2000; Nauton-Fourteu et al., 2021). Basin fill mainly consists of the thick siliciclastic sequences of the Shannon and Central Clare groups, which preserve deep-water (Ross Sandstone), slope (Gull Island Formation), and shelf and fluvial-deltaic (Tullig, Kilkee, Doonlicky, and younger cyclothems) deposits of the Western Irish Namurian Basin (Fig. 2). These cyclothems typically comprise alternating marine and continental successions and are equivalent to formations (Wanless and Weller, 1932). The Gull Island Formation and Tullig Cyclothem are bounded by two successive marine bands (R. dubium and R. stubblefieldi), which represent a duration of 65–109 k.y. (Best and Wignall, 2016; Fig. 2). After the basin was filled, a series of east–west-trending, kilometer-scale, thrust-related folds developed in association with the Variscan orogeny; the tectonic event deformed the upper Paleozoic sedimentary succession within the Western Irish Namurian Basin (Figs. 3F and S1; see ^{footnote 1}; Cooper et al., 1986).

The present study focuses on the multistory fluvial channel belts in the Tullig Sandstone of the Tullig Cyclothem (Pulham, 1989; Stirling, 2003) (Figs. 2 and 3A–3E). Five main exposures of the entire Tullig Sandstone stratigraphic column are located along the coast of western Ireland at Trusklieve, Pulleen, Killard, Carrowmore Point, and Furreera (Figs. 1 and 3). These locations span over 40 km along the coastline and are generally aligned in a southwest-to-northeast direction, which is also parallel to the dominant paleocurrent direction (Stirling, 2003). While also exposed at Tullig Point and Castle Point (Fig. 3F, Table S1; see ^{footnote 1}), the Tullig Sandstone is inaccessible at these locations due to the nature of the cliff exposures.

The Tullig Sandstone consists mostly of sandstone lithofacies from bar deposits (Stirling, 2003). Gravel is not present in these exposures of the Tullig Sandstone. Channel stories were identified based on erosional relationships and fine-grained (mudstone to siltstone) overbank or channel-fill facies by Stirling (2003). The number of channel stories present at each of the five key sites listed above decreases from southwest to northeast: there are five stories at Trusklieve, three at Pulleen, two at Killard and Carrowmore Point, and one at Furreera (Best and Wignall, 2016; Stirling, 2003) (Figs. 3A–3E and Table S1). The maximum channel-story thickness of all exposures ranges from 1.6 m to 22 m, with an average of ~9 m (Stirling, 2003). The thickness of the Tullig Sandstone also tapers basinward (i.e., northeastward; Fig. 2), from 31.2 m at Trusklieve, 22.5 m at Pulleen, 20 m at Killard, 22 m at Carrowmore Point, to 5 m at Furreera (Pulham, 1989; Stirling, 2003). The outlet of the paleo-river or the location of the paleo-shoreline is thus interpreted to be close to Furreera, since measured paleocurrent directions are dominantly NNE (Pulham, 1989; Stirling, 2003; Nauton-Fourteu et al., 2021) and Tullig Sandstone is absent in Western Irish Namurian Basin exposures several kilometers to the north, specifically at the Cliffs of Moher (Best and Wignall, 2016).

The sequence stratigraphic framework of the study area shows that the Tullig Sandstone is overlain by the transgressive succession of the Tullig Cyclothem (Obrock, 2011; Fig. 2), which suggests a history of base-level rise. Moreover, the lack of transgressive deltaic deposits overlying the Tullig Sandstone and the presence of an extensive flooding surface indicates rapid non-deltaic shoreline retreat (i.e., sediment-starved autoretreat; Muto and Steel, 2001; Muto, 2001), where the delta foreset is abandoned and marine deposits directly overlie the fluvial deposits (Wu et al., 2020a).

Stratigraphy of the Western Irish Namurian Basin is first restored, and then an inversion method is developed to reconstruct channel-bed slope and flow depth using the structurally restored geometry of the entire Tullig Sandstone. This method also predicts the extent of backwater influence on stratigraphic architectures. Then, field data (i.e., measurements of sandstone grain size, cross set, and bed thicknesses) were collected to assess their spatial variabilities in relation to backwater. Finally, previous paleohydraulic reconstruction methods developed for normal flow conditions were applied to suitable field locations to test the results of the inversion method.

Our structural restorations use systematic measurements of stratigraphic boundaries and geological mapping of tectonic deformation (i.e., folds and faults) to characterize accurately the geometries and orientations of pre-Variscan basin stratigraphy (see Supplemental Material and Figs. 3 and S1). The results of this analysis provide the spatial control necessary to reconstruct pre-deformation distances between the key bedrock exposures and the down-dip length of the Tullig Sandstone.

Geologic map data were used to compile six geologic cross sections oriented perpendicular to the regional structural trends (northeast) (Fig. S1). Line-length restorations of these cross sections suggest that 5–7% total horizontal shortening of the stratigraphic succession was accommodated by thrust faults and fault-related folding in a NW–SE direction. A parallel analysis of grain-scale rock fabrics conducted by Lopez et al. (2017) on samples collected from the Ross Formation indicate that the contribution of penetrative strain to regional, layer-parallel shortening of strata within the Western Irish Namurian Basin is negligible. These analyses suggest minimal shortening on the downslope direction (northeast). Therefore, the resulting geologic map depicts structures and stratigraphic contacts in sufficient detail to compare observed stratigraphic architectures and model results from the inversion method.

A stratigraphic inversion method was developed by evaluating the length and thickness scales of fluvial sandstones that characterize the typical basinward-tapering wedge geometry of fluvial-deltaic sandstones (Fig. 4; e.g., Armstrong, 1968). The inversion method utilizes a backwater morphodynamic model to simulate one-dimensional and dip-view stratigraphic development through the evolution of a channel longitudinal profile. The first-order length and thickness of simulated stratigraphy are then compared with observed Tullig Sandstone geometry to refine estimates of channel-bed slope and flow depth.

The backwater morphodynamic modeling framework of Parker et al. (2008a, 2008b) is used to simulate stratigraphic development. First, a backwater equation (Chow, 1959; Parker, 2004) provides a means to model bankfull flow depth (H_bf) in the along-stream distance (x):

where S is channel-bed slope, C_f is friction coefficient (0.0016; Parker et al., 2008b), Fr is Froude number, determined by Fr = U(gH_bf)^−0.5, with U the depth-averaged flow velocity and g gravitational acceleration. Depth-averaged flow velocity may be estimated by: U = q_w / H_bf, where q_w is water discharge per unit width and can be calculated by: (Parker, 2004).

Here, the bed material load equation of Ma et al. (2017) is used:

where R is the submerged specific gravity of sediment particle and τ_* is the non-dimensionalized shear stress, determined by τ_* = C_fU² / RgD₅₀.

The along-stream variation in sediment flux (∂q_s/∂x) is calculated and used to assess the rate of change in channel-bed elevation through time (∂η/∂t) using a modified Exner equation (Paola and Voller, 2005):

where A is determined by A = (1 + Λ)ΩI_f(1−λ_P)⁻¹r_B⁻¹, and λ_p is the mean porosity of the channel-floodplain complex, Λ is the mud/sand ratio, Ω is the channel sinuosity, I_f is the flood intermittency, and r_B is the ratio between channel width and width of the flood plain. Typical values for these parameters, as measured from modern lowland systems (Parker et al., 2008b; Moran et al., 2017), are listed in Table 1. Equations 1–3 are solved iteratively to simulate the channel-bed elevation at each time step over the along-stream distance (x), which thus produces longitudinal fluvial-deltaic stratigraphy through time (Wu, 2020; Wu and Nittrouer, 2020).

The strategy of the stratigraphic inversion method is to run a large number of stratigraphic simulations with possible ranges of input parameters. The misfit between simulated stratigraphic geometries and the measured geometries of the Tullig Sandstone, in terms of first-order sandstone thickness, L_t, and length, L_l (Fig. 4), can then be calculated to determine the best-fit set of model parameters.

In this study, we evaluate input parameters including channel-bed slope, S, rate of base level change, R_bl, and ratio between channel width and width of flood plain, r_B. The input parameter space for S, R_bl, and r_B was defined as 6 × 10⁻⁵–4 × 10⁻⁴, 0.5–10 mm/yr, and 10–50, respectively. The parameter space for S, R_bl, and r_B was then discretized into a 50 × 50 × 50 grid, which corresponds to 125,000 unique combinations of S, R_bl, and r_B that are used as inputs for the backwater morphodynamic model. The upstream boundary condition of normal flow depth, H_n (i.e., bankfull depth of steady and uniform flow of the river reach located upstream of the backwater zone), is calculated using the Li et al. (2015) relation (combining their equations 7B and 8B):

where ν denotes kinematic viscosity of water and is assumed to take the value 1 × 10⁻⁶ m² s^–1. Flow depth—a downstream boundary condition—is calculated as the difference between channel-bed elevation and base level, which is adjusted for each time step based on the value of R_bl (Wu et al., 2020a).

Similar to the definition of the misfit function for the first-order shelf geometry from Yuan et al. (2019), a misfit function, μ, for the first-order geometry of fluvial sandstone is defined as:

where is the maximum simulated thickness of fluvial stratigraphy, and is the distance from the location of maximum stratigraphic thickness to the downstream end of the fluvial stratigraphy (Fig. 4C), and and are the observed thickness and length of the Tullig Sandstone. is set as the Tullig Sandstone thickness at Trusklieve (31 m), since this location is the most upstream and thickest section of the Tullig Sandstone (Pulham, 1989; Stirling, 2003). The structurally restored distance from Trusklieve to Furreera is 43 km. Since the Tullig Sandstone is absent at the Cliffs of Moher (~3.5 km to the north; Fig. 3F), it is assumed that the distal reach of the Tullig Sandstone is located between Furreera and the Cliffs of Moher. Therefore, the length of the Tullig Sandstone, , is set as 45 km. The log-likelihood function is calculated as:

An optimal solution of the input parameters that corresponds to least misfit (maximum likelihood) was considered to be the best-fit model for the observed stratigraphic geometry. Importantly, posterior probability distributions of the input parameters evaluated can be constructed using the following relation (Charvin et al., 2009):

where X represents the parameter evaluated (i.e., S or R_bl), P(X | D_obs) is the posterior probability of the parameter evaluated given the observed data (D_obs, i.e., thickness and length of fluvial sandstone in this study), P(D_obs | X) is the maximum likelihood estimate of the parameter evaluated given observed data and is calculated as the log-likelihood L, and P(X) is the prior probability of the parameter evaluated. Since the morphodynamic model is computationally efficient, the parameters evaluated were estimated using a grid-search method. Rather than assuming a uniform distribution for the prior, here the prior probability was calculated as the kernel density for the frequency of each input value of a parameter that corresponds to a minimum misfit model. For example, for a given input value of channel-bed slope, the prior probability of this slope was calculated as the kernel density of frequency of the slope that corresponds to a minimum misfit model, given the ranges of other parameters (i.e., R_bl and r_B). Credible intervals were then calculated based on the posterior probability of the inverted parameters, which can be compared with the results from paleohydraulic methods that assume normal flow conditions. The inversion method also provides a means to test flow depths, channel bed slopes, and aggradation patterns inferred from the coefficient of variation of cross-set thickness (Jerolmack and Mohrig, 2005; Wu et al., 2020b). It is inferred that since the Tullig Sandstone is overlain by a marine transgressive succession, abandonment of deltaic foresets resulted from shoreline retreat (i.e., autoretreat; Muto, 2001; Wu et al., 2020a). As such, only the simulations where shoreline autoretreat developed through the model run were evaluated.

The upstream extent of backwater influence, considered to be the location of normal flow to backwater transition, is recorded for all simulations. Time (f_ti) and thickness (f_th) fractions of the modeled stratigraphy, developed within the backwater zone and set by the backwater length scale, L_bw, were then calculated. Specifically, the time fraction, f_ti, at a vertical section in the stratigraphy is calculated as the ratio between the total time when this section is located downstream of the backwater transition to the total time when there is an active channel present at this section. The thickness fraction, f_th, is calculated as the ratio between the total thicknesses when a section is developed within the backwater zone to the total thickness of this section. Time and thickness fractions from the best-fit simulation provide predictions of the spatiotemporal impact of backwater through the stratigraphy.

To evaluate the impact of in-channel deposition (channel bed aggradation) on stratigraphy (Parker et al., 2008a; Ganti et al., 2016), particularly dune-scale cross-set development (Wu et al., 2020b), we calculated the time-averaged depositional front from simulation results. The depositional front is identified as the location of maximum aggradation at each time step in the morphodynamic model (Equation 3). The time-averaged depositional front is taken as the average location of the depositional front over the duration of each model run. The coefficient of variation of cross-set thickness measured from stratigraphic sections can infer the location of elevated channel-bed aggradation (Jerolmack and Mohrig, 2005; Leary and Ganti, 2020; Wu et al., 2020b), which represents a maximum in average channel aggradation rather than an instantaneous depositional front from a particular time. The depositional front alone is insufficient for characterizing the average aggradation rates across an entire stratigraphic section, since it is not static in nature (Moran et al., 2017; Wu and Nittrouer, 2020). Therefore, the time-averaged depositional front from simulations provides a means to test the inferred pattern of channel aggradation from field observations, for example, cross-set thicknesses or bar thicknesses.

Stratigraphic sections were measured at each of the five key exposure locations (Trusklieve, Pulleen, Killard, Carrowmore Point, and Furreera) in the Tullig Sandstone to document sedimentary structures and thicknesses. A total of 67 rock samples from bar deposits were collected across all channel stories identified in a previous study (Stirling, 2003). Two to three samples were collected for each channel story to characterize the lower, middle, and upper portions of each story and account for different sandstone lithofacies. For example, the base, middle, and top of a channel story typically consist of massive sandstones, cross-bedded sandstones, and ripple cross-laminated sandstones, respectively. Thin sections prepared from these rock samples were analyzed for grain-size distributions using the Gazzi-Dickens point counting method (Folk, 1980; Ingersoll et al., 1984) in the software MicroPet™. An evenly spaced grid was generated in MicroPet™ to provide a minimum of 300 points of grain size measurements for each thin section. The grain-size distribution for the top story of the Tullig Sandstone was also analyzed separately to test the hypothesis that grain size of fluvial sandstone overlain by a regional flooding surface fines up-dip due to autoretreat (Wu et al., 2020a).

Outcrops of the Tullig Sandstone were surveyed by an unmanned aerial vehicle (UAV, DJI Mavic Pro Platinum™). Three-dimensional digital outcrop models were generated from images collected from UAV using Agisoft Metashape™ (Wu, 2022) (Supplemental Material). The thicknesses of cross-set and sandstone beds were measured in the field and from the 3-D outcrop models (Supplemental Material). A single measurement for cross-set thickness is collected for each set that is laterally continuous and uniform in thickness, which is typical for the study area. For sets that are laterally variable in thickness, multiple measurements were made, and an average value was determined. A “bed” in this study is defined as the coset of cross sets, which represents a group of cross sets that are separated by second-order bounding surfaces (Miall, 1988). Therefore, a bed or coset is typically several times the thickness of an individual cross set.

A fixed-effects analysis of variance (ANOVA) with unpooled variance was fit to the data using Markov Chain Monte Carlo methods (Venables and Ripley, 2002) to estimate the mean and standard deviation values of cross-set/bed thickness at each outcrop location using uninformative priors. The posterior distributions of the mean and standard deviation were used to interpret changes in cross-set/bed thickness, as well as the coefficient of variation in thickness down-dip.

Paleohydraulic reconstructions were conducted for the individual channel stories at suitable sites where normal flow conditions are expected, based on the result of the inversion analysis. Set thickness measurements for each channel story were grouped into two categories: cross-bedded and planar/low-angle bedded sets. The distribution of mean set thicknesses was approximated by random resampling with replacement (i.e., bootstrapping method). The average grain size for each channel story is calculated as the mean of median grain size for samples collected from that story. The bankfull flow depth, H_bf, and channel bed slope, S, were assessed for each individual set category and then for all set measurements combined using paleohydraulic reconstruction methods for normal flow conditions.

We used the method of Bradley and Venditti (2017) to estimate flow depth by using a scaling relationship to mean dune height (h_d):

Dune height was estimated as (Leclair and Bridge, 2001):

where s_m is the average cross-set thickness based on measurements from the field.

Slope for the normal flow reach was then estimated using the Lynds et al. (2014) method (their equation 4):

where R is the submerged specific gravity of sediment (1.65), is the bankfull Shields number (1; Lynds et al., 2014), and D₅₀ is the median grain size of bed material.

For comparison, slope was also calculated using a method from Trampush et al. (2014):

where α₀, α₁, and α₂ are empirical coefficients taking the mean values of −2.08, 0.254, and −1.09.

Equations 8–11 are intended for normal flow conditions and cannot be applied where backwater flow persists, in particular near the river outlet (Wu et al., 2020b). Therefore, we first estimated flow depth and channel-bed slope using the inversion method, then calculated backwater length, L_bw. The locations of Tullig Sandstone exposures in reference to the backwater length can be determined with the location of Furreera treated as the distal reach of the river outlet, since Tullig Sandstone is absent north of Furreera. Equations 8–10 were then applied to exposures outside the backwater zone, where normal flow persists.

The inversion analysis shows: (1) the best-fit model (i.e., least-misfit model run) corresponds to a channel-bed slope of 2.13 × 10⁻⁴ (with a normal flow depth of 4.2 m), and a rate of base-level rise of 5.9 mm/yr (Fig. 5). The best-fit model suggests a backwater length scale of 20 km. (2) The 95^th percentile credible interval for estimated channel-bed slope is 1.6–3.1 × 10⁻⁴, with corresponding normal flow depths ranging from 3.6 m to 4.8 m. (3) The 95^th percentile credible interval for the estimated rate of base-level rise is 1.1–9.6 mm/yr. Based on the inverted flow depth and channel-bed slope (the 95^th percentile credible intervals), the backwater length scale is calculated to be 11–30 km.

The analysis shows that the normal flow to backwater flow transition, determined by L_bw, varies spatiotemporally and is mostly located downstream of Trusklieve (Figs. 6B–6C). The location of the exposures of Tullig Sandstone at Trusklieve (43 km upstream from the downstream exposure at Furreera) would have experienced minimal impact of non-uniform flow compared to the rest of the Tullig exposure located within the range or near the limit of the backwater length scale (Figs. 6B–6C and 7A). Both the predicted time and thickness fractions at Trusklieve with respect to backwater influence are less than 0.3, based on the best-fit model run (Fig. 7A). The predicted times and thickness fractions of backwater influence for the rest of the exposure locations are all above 0.5 (Fig. 7A): Pulleen (0.61 and 0.54), Killard (0.73 and 0.62), Carrowmore Point (0.91 and 0.89), and Furreera (1 and 1). This result indicates that the exposures of the Tullig Sandstone in this study are influenced by backwater hydrodynamic conditions.

Using the inversion analysis, the time-averaged depositional front is quantified across the simulated stratigraphy for a group of low-misfit model runs (Fig. 7B). The simulated locations of time-averaged depositional front cluster around 20 km, near Killard, and skew upstream to Trusklieve. This distribution is reasonable as the shoreline advances from Trusklieve (0 km) at the start of simulation and reaches Furreera (43 km) before retreating to Trusklieve (Fig. 6). Therefore, the total amount of time for active in-channel deposition decreases downstream (Fig. 7A).

Channel fills in the Tullig Sandstone comprise five lithofacies (LF): (LF-1) massive sandstone, (LF-2) cross-bedded sandstone, (LF-3) planar to low-angle bedded sandstone, (LF-4) ripple-cross laminated sandstone and siltstone, and (LF-5) laminated or massive mudstone (Fig. 8). Massive beds of sandstone (LF-1) that occasionally contain rip-up clasts of mud are usually found at the base of channel fills, while the tops of channel fills are typically characterized by bioturbated, ripple-laminated beds of sandstone (LF-4). Mudstone lithofacies (LF-5) are sparsely distributed among channel fills with horizontal laminations or massive appearance.

The Tullig Sandstone is dominated by sandy lithofacies with limited mudstone content (Fig. 9A–C). The most common lithofacies in the Tullig Sandstone are cross-bedded sandstone (LF-2) and planar to low-angle bedded sandstones (LF-3). Together with massive sandstone, sandstone facies (LF-1–3) typically comprise more than 70% of the total volume of Tullig Sandstone. Fine-grained lithofacies (LF 4) and mudstone facies generally decrease in abundance downstream (Fig. 9E). However, we observe a distinct abundance of ripple-cross laminated sandstones (LF-4) at Carrowmore Point (~25 km downstream of Trusklieve; Fig. 9D).

The grain-size distributions of samples collected (Fig. 10) possess median sizes ranging from 100 μm to 200 μm. While the average of the median sizes of all samples at each location (~150 μm) does not demonstrate a spatial trend, the coarsest fraction (99^th percentile) of the grain-size distributions decreases down-dip, from Trusklieve (362 μm) to Furreera (268 μm). To characterize the degree of sediment sorting, the sorting index, σ_g, is calculated using (Folk, 1980): , where ϕ₈₄, ϕ₁₆, ϕ₉₅, and ϕ₅ represent the 84^th, 16^th, 95^th, and 5^th percentiles of grain size. The sorting index decreases from Trusklieve (1.296) to Furreera (1.2584) (Fig. 11).

The grain-size distribution for the top story of the Tullig Sandstone was also analyzed separately. Samples collected at both the flooding surface (top of the story, Figs. 8F and 12A) and base of the top channel story (Fig. 12B) show an up-dip fining trend from Killard to Tullig Point. In particular, the median grain size of samples from the base of the channel story fines from 300 μm at Killard to 150 μm at Tullig Point. Bioturbation at the contact between the flooding surface and the top of the youngest channel story precludes the identification of any trends of size distribution (Fig. 12A).

Thicknesses of dune-scale cross sets in the Tullig Sandstone range from less than 10 cm to 90 cm. The average cross-set thickness (s_m) increases downstream, from 12.5 cm at Trusklieve to 41 cm at Carrowmore Point, and then decreases to 18 cm at Furreera (Fig. 13A). This spatial trend in estimated average cross-set thickness is accompanied by an opposing trend in the estimated coefficient of variation in cross-set thickness (C_v–s = s_m/s_std–s, where s_std–s is the standard deviation of cross-set thickness), which decreases from Trusklieve (0.83) and reaches a minimum at Carrowmore Point (0.39) (Fig. 13B). Bed thickness of the Tullig Sandstone shows a spatial trend similar to that of the dune cross sets: mean bed thickness (s_b) increases from 0.55 m at Trusklieve, peaks around 1.87 m at Killard, and decreases to 0.43 m at Furreera (Fig. 13C). The coefficient of variation in bed thickness (C_v–b = s_b/s_std–b where s_std–b is the standard deviation of bed thickness) also decreases from Trusklieve (0.47) and reaches 0.32 at Killard, which is ~4 km up-dip from Carrowmore Point (Fig. 13D). The 95^th percentile credible intervals of cross-set and bed thickness are significantly different for each location (Figs. 13A and 13C). The 95^th percentile credible interval for coefficient of variation in cross-set thickness also shows a spatial trend that reaches a minimum at Carrowmore Point, although the credible intervals partly overlap, especially at Carrowmore Point and Furreera (Fig. 13B). The 95^th percentile credible interval for coefficient of variation in bed thickness still shows a minimum at Killard, but the credible intervals largely overlap at each location (Fig. 13D).

Trusklieve is the least impacted exposure in terms of non-uniform flow because its location from the channel terminus—43 km—is greater than the backwater length calculated by the inversion analysis (i.e., 11–30 km). Paleohydraulic reconstruction methods (i.e., Equations 8–10) are therefore only applied using measurements collected from the Tullig Sandstone exposure at Trusklieve, due to its minimum backwater influence, based on calculated time and thickness fraction values that were determined to be 0.29 and 0.24, respectively (Fig. 7A).

The paleohydraulic reconstruction method for normal flow conditions was then applied based on measured cross-set thickness at Trusklieve (range: 2–51 cm, average: 12.5 cm). The method provides estimates of mean dune height of 28–44 cm and the bankfull flow depth of 1.9–3 m (Equations 8 and 9). Using these flow-depth estimates, measured grain sizes, and the methods described by Equations 10 and 11, estimated channel slopes range from 0.8 × 10⁻⁴ to 1.4 × 10⁻⁴ and 2.7 × 10⁻⁴ to 4.5 × 10⁻⁴, respectively. The thickness of planar to low-angle beds can also be used to invert for the formative bedform geometry (Bridge, 1997). This technique typically provides higher estimates of flow depth (3.2–7.8 m, see stories 3 and 5 in Fig. 14) and correspondingly lower estimates of channel-bed slope (0.3 × 10⁻⁴–0. 8 × 10⁻⁴ and 0.9 × 10⁻⁴–2.5 × 10⁻⁴; using Equations 10 and 11, respectively). Considering all cross sets and planar to low-angle bed measurements together, flow depth is estimated to be 3.3 ± 0.3 m, and channel bed slope is estimated to be 0.7 × 10⁻⁴ (Equation 10) and 2.4 × 10⁻⁴ (Equation 11).

The conundrum of paleohydraulic reconstruction for regions of non-uniform flow lies in the fact that standard methods for calculating bankfull flow depth, H_bf, and channel-bed slope, S, are only applicable where uniform flow conditions exist (Equations 8–11); yet, to determine the extent of backwater flow (L_wb), flow depth and bed slope need to be known. Moreover, base-level adjustments can shift the shoreline and backwater region, and as a consequence, the zone of non-uniform flow is not fixed in space or time (Parker et al., 2008a; Moran et al., 2017; Chadwick et al., 2020).

The stratigraphic inversion method developed in this study allows for a quantitative reconstruction of paleo-flow conditions (i.e., where normal and non-uniform flow conditions prevailed) of ancient fluvial-deltaic systems, and thus resolves this conundrum. The method is advantageous because it makes use of the first-order geometric properties that can be measured directly from the stratigraphic record, specifically, thickness and length of the basinward-tapering sandstone wedge. Using fluid-flow and sediment-transport relations, the method generates estimates of regional channel-bed slope and normal flow depth.

Once bankfull flow depth and channel-bed slope for normal flow conditions are calculated using the inversion method, a backwater length scale can be estimated to determine where it is appropriate to use paleohydraulic reconstruction techniques. It is assumed that the sedimentological and stratigraphic data collected identify the distal extent of channel deposits, and more specifically, that the pinch-out of the sandstone wedge coincides with the channel terminus. Comparing the independent estimates of bankfull flow depth and channel-bed slope using the inversion method and standard methods will further refine paleohydraulic reconstruction results. For example, both the inversion method and paleohydraulic reconstruction method for normal flow generate similar results for estimates of channel depth for the Tullig Sandstone (i.e., 3–4 m). The reconstructed channel-bed slopes predicted using the inversion method and the Trampush et al. (2014) method also agree (2.13 × 10⁻⁴ and 2.4 × 10⁻⁴, respectively). On the other hand, the Lynds et al. (2014) method shows discrepancy with the other two methods, producing a slope value that is a factor of three lower (i.e., 10⁻⁵), a value that is typical for large, lowland continental-scale rivers. Given the estimated size of the Western Irish Namurian Basin and the nature of its opening (transtensional basin; Best and Wignall, 2016), a slope value in the range produced using the Trampush et al. (2014) slope method is likely a more reasonable value.

It is recognized that the Lynds et al. (2014) method for estimating channel-bed slope possessing fine- to medium-grained sand (i.e., ~250–1000 μm) could produce an underestimated value due to the choice of the value for the Shields number. Lynds et al. (2014) proposed using a lower bound to minimize errors, although Wilkerson and Parker (2011) indicate that the Shields number for suspended load-dominated rivers ranges from 1 to 10. We followed Lynds et al. (2014) (i.e., lower bound of Shields number), and due to the fine grain size of Tullig Sandstone, Equation 10 generated a relatively low channel-bed slope value. The Trampush et al. (2014) method, on the other hand, is an unconstrained and empirical relation without prescribed scaling between S, H_bf, and D₅₀. Thus, this method avoids the need to estimate a Shields number for fine-grained channel beds.

The spatial variability in the architecture of the Tullig Sandstone supports the hypothesis that the hydrodynamic influence of backwater flow can be interpreted from the stratigraphic record. The amalgamated channel bodies at Trusklieve indicate that the Tullig Sandstone was composed of multiple channels, and this architecture could reflect the effects of spatiotemporally varying backwater hydrodynamics (Wu and Nittrouer, 2020). Trusklieve, the most up-dip exposure, is located 43 km upstream of the inferred distal paleoshoreline, a distance that is greater than the backwater length scale (11–30 km) estimated by the inversion analysis. Therefore, Trusklieve is likely to be situated near the upstream avulsion node, as previous studies have shown that avulsion length scale L_a, which measures the river distance from shoreline to the avulsion location, is related to the backwater length: L_a = 0.5 L_bw –2 L_bw (Chadwick et al., 2020). The presence of multiple cut-and-fill deposits at Trusklieve further supports the supposition that this location experienced multiple avulsions (Ganti et al., 2014, 2016; Chadwick et al., 2019). Alternatively, the cut-and-fill deposits at Trusklieve may be the result of enhanced lateral mobility of the paleochannel, which is also common for natural fluvial-deltaic systems near their backwater transition, a location that typically coincides with the delta apex (Hudson and Kesel, 2000; Kleinhans et al., 2011; Nittrouer et al., 2012).

This interpretation is also supported by the evolution of the channel profiles, which was modeled for the condition of base-level change (Fig. 6). For example, the channel-bed slope in the backwater reach adjusts as the shoreline migrates basinward over multiple avulsion cycles (Moodie et al., 2019; Gao et al., 2020; Wu and Nittrouer, 2020), creating a concave-up bed profile. For modern river systems, including the Mississippi (Nittrouer et al., 2012), Po (Maselli et al., 2018), and Trinity (Smith et al., 2020), as well as has been observed in physical experiments (e.g., Ganti et al., 2016), bed profiles possess a concave-up appearance, with a downstream decrease in slope nearing the receiving basin. For example, near the outlet, the bed slope is typically around half the value of that of the normal flow region upstream (Nittrouer et al., 2011; Smith et al., 2020). This geometry promotes upstream migration of the depositional front, which leads to channel backfilling (Ganti et al., 2016) and eventually a channel avulsion. Moreover, base-level rise typically results in upstream migration of the avulsion node (Moran et al., 2017; Chadwick et al., 2020) and drowns the river mouth, abandoning the delta and leading to rapid shoreline retreat (i.e., sediment-starved autoretreat, Fig. 6C; Muto, 2001; Wu et al., 2020b).

The overall grain-size statistics (Figs. 10–11) from the Western Irish Namurian Basin demonstrate downstream fining by way of a reduction in the coarsest fraction of bed sediment, which may be reflective of decreasing bed shear stress as a consequence of backwater effect, as is found for modern systems (Nittrouer et al., 2012; Smith et al., 2020). This is consistent with numerical simulations, which have shown that an overall down-dip decrease in median grain size can be preserved in stratigraphic sections despite shoreline movement and base-level change (Wu and Nittrouer, 2020).

We interpret the up-dip fining trend for the top story of the Tullig Sandstone from Killard to Tullig Point (Fig. 12) to be a consequence of base-level rise and progressive upstream shoreline retreat (i.e., autoretreat; Fig. 6C; Muto, 2001; Parker et al., 2008a). The up-dip fining of sediment for the top story of the Tullig Sandstone therefore reflects the cumulative effect of backwater-induced downstream fining and shoreline movement via autoretreat (Wu et al., 2020a).

Cross-set thickness is hypothesized to vary systematically within the backwater zone, which can serve as a stratigraphic signature for backwater hydrodynamics (Wu et al., 2020b). As measured for the Tullig Sandstone, cross-set thickness varies down-dip. From Trusklieve to Carrowmore Point, thickness increases; however, from Carrowmore Point to Furreera, thickness decreases (Fig. 13A). This pattern of down-dip increase to decrease in cross-set thickness is diagnostic of backwater hydrodynamics (Wu et al., 2020b). Specifically, cross-set thickness is a function of multiple variables, including dune wavelength (l), rate of aggradation (r), dune celerity (c), mean dune height (h_d), and a coefficient of variation of dune height (C_v–d) (Paola and Borgman, 1991; Bridge and Best, 1997; Cardenas et al., 2019):

Backwater hydrodynamics modulate cross-set thickness mainly by affecting, r/c, h_d, and C_v–d. For example, dune height (h_d) is affected by the shear velocity (u*) of the flow and settling velocity (w_s) of bed material (Best, 2005; Bradley and Venditti, 2019; Cisneros et al., 2020), where both values are expected to vary down-dip due to backwater flow and reducing channel-bed grain size. A maximum dune height is reached when the suspension criterion is met (u*/w_s ≈ 1; Bradley and Venditti, 2019) and can result in a maximum cross-set thickness (Equation 12). Similarly, r/c can also reach a maximum within the backwater zone and promote the development of a maximum cross-set thickness (Equation 12). Specifically, dune celerity (c) should decrease downstream within the backwater zone as sediment flux, which scales with dune celerity (q_s = 0.5 hc(1–λ_p); Simons et al., 1965), is reduced (Nittrouer et al., 2012). Sediment aggradation (r) within the backwater zone is typically elevated (Nittrouer et al., 2012; Ganti et al., 2016; Moran et al., 2017), and so the combined impact of elevated aggradation and downstream decrease in dune celerity maximizes r/c in the backwater zone (Wu et al., 2020b). Additionally, it is worth noting that cross-set preservation within the backwater zone may simultaneously be influenced by non-uniform conditions such as flashy flood hydrographs or lateral migration of barforms (Ganti et al., 2020; Leary and Ganti, 2020; Lyster et al., 2022).

The coefficient of variation in cross-set thickness serves as an indicator of bed aggradation: enhanced aggradation reduces the variability in scour depth, thus promoting less variable cross-set thickness and a smaller coefficient of variation values (Paola and Borgman, 1991; Jerolmack and Mohrig, 2005). Therefore, the data suggest that channel-bed aggradation peaks at Carrowmore Point (Fig. 13B). This interpretation is also in accordance with the observation of a higher percentage of ripple cross-lamination, especially climbing ripples, at the same site (Fig. 9D), which suggests enhanced channel-bed aggradation (Ashley et al., 1982; Bridge, 1997). It has also been noted in physical experiments that the coefficient of variation of cross set C_v–s, with a value of ~0.8 (±0.3), is indicative of steady-state conditions (Leary and Ganti, 2020). Therefore, a C_v–s value of <~0.58 suggests that non-uniform conditions dominate the sections downstream of Killard (Fig. 13B), which is also supported by the estimated spatial extent and duration of backwater influence on stratigraphy (Fig. 7A).

The challenge associated with linking predicted cross-set thickness patterns and ancient stratigraphy lies in the temporal variability expected as a result of backwater hydrodynamics. Specifically, the nature of sediment deposits and ensuing stratigraphy expected to develop due to backwater conditions adjusts with a dynamic shoreline and evolving channel-bed profile (Wu and Nittrouer, 2020; Chadwick et al., 2020). The inversion analysis presented herein suggests that elevated aggradation can be recorded in the stratigraphy of channel deposits (Fig. 6B). Moreover, the results of model simulations indicate that the locations of time-averaged depositional fronts cluster over 7–22 km downstream of Trusklieve for simulations that use a range of channel-bed slopes and base-level conditions (Fig. 6B). Therefore, the spatial pattern in aggradation, which increases and then decreases down-dip, could be reflected in cross-set thicknesses, with a maximum value expected near the location of the time-averaged depositional front.

The observed pattern in cross-set thickness aligns with the modeled depositional front in our study. The best-fit simulations from the inversion analysis (Figs. 6 and 7) suggest that the time-averaged depositional front is located 4 km upstream of Carrowmore Point, where cross-set geometries (Fig. 13) and the high percentage of climbing ripples (Fig. 9D) also indicate channel-bed aggradation. The distance between Carrowmore Point and the predicted location of time-averaged depositional front from the best-fit simulation is ~9% of the total length of the Tullig Sandstone (43 km from Trusklieve to Furreera) analyzed in this study. Considering the limited geological information (i.e., and ) used in the inversion analysis, this result is consistent with field observations and highlights the high preservation potential for the signature of backwater hydrodynamics in stratigraphy.

The statistical and analytical theories underlying cross-set thickness and length are scale-independent (Paola and Borgman, 1991; Bridge and Best, 1997); as such, they can be applied to a range of bedform sizes (e.g., ripples to dunes; Cardenas et al., 2019). We hypothesize that the statistics of bed thickness and length measured for the Tullig Sandstone not only reflect the propensity for higher r/c values due to aggradation but are also affected by river kinematics, specifically lateral channel migration, which provides the sufficient local accommodation required for bed preservation (van de Lageweg et al., 2016; Ganti et al., 2020). A lower coefficient of variation in bed thickness indicates elevated aggradation, similar to the coefficient of variation of cross-set thicknesses. Lateral channel movement influences cross-set and bed thickness preservation potential by developing scour-fill successions and local aggradation (Lindsay and Ashmore, 2002; van de Lageweg et al., 2016). Lateral migration rates are elevated at the backwater transition of rivers and decrease downstream (Ikeda, 1989; Hudson and Kesel, 2000) due to decreasing reach-average velocity and bank erodibility as the mud content in the bank sediment increases (Li et al., 2020). The stratigraphic signature of elevated sandstone-bed thickness around the backwater transition could be common for ancient systems, as the increase in lateral migration of the channel around backwater transition is suggested by observations of an accompanying increase in river sinuosity from lowland river systems (Wu et al., 2021). Finally, sandstone bed lengths measured in a previous study (Stirling, 2003, their fig. 5.13) show an identical trend (i.e., maximum sandstone bed length at Carrowmore Point) to sandstone bed thickness measured from this study, which suggests consistent spatial variability in 3-D bed geometry that is possibly influenced by the local aggradation rate (Jerolmack and Mohrig, 2005).

The inversion method presented in this study compares the measured and simulated first-order geometry of the fluvial sandstone wedge. The latter presents a source of uncertainty for the method due to the natural variabilities in hydraulic parameters used in the morphodynamic model, including friction coefficient, C_f, mud/sand ratio, Λ, river sinuosity, Ω, flood intermittency, I_f, the ratio between channel width and floodplain width, r_B, and the rate of base-level rise, R_bl. However, the uncertainties associated with the predicted channel-bed slope and flow depth using the inversion method are limited despite the natural variabilities in these hydraulic parameters. This is because the first-order control on the geometry of a sandstone wedge is channel-bed slope. To demonstrate this, we analyzed the key model parameters and their impacts on stratigraphic development in our model (Supplemental Material).

First, mud/sand ratio, Λ, river sinuosity, Ω, flood intermittency, I_f, and ratio between channel width and floodplain width, r_B, appear in Equation 3 and can be treated as a single parameter, A. Parameter A is a scaling term and will affect the rate of change in channel bed elevation through time (∂η/∂t). Therefore, parameter A only determines how fast the stratigraphy develops and cannot affect the overall geometry of the stratigraphy (Parker, 2004). Simulations with different values of A can produce very similar stratigraphy contingent on the model run times and rate of base-level rise (Supplemental Material).

Second, the uncertainty in modeled stratigraphy is mainly due to the natural variability in friction coefficient and bankfull depth of normal flow (Figs. S2A–S2B; see ^{footnote 1}). However, sensitivity analyses using friction coefficient and bankfull depth suggest that the inversion method is robust even when faced with uncertainties (see Supplemental Material). For example, the slope prediction shows only a 9% increase with a C_f value that is a factor of three higher (Fig. S3A; see ^{footnote 1}). Depth prediction shows only a 5% increase for best-fit scenario from 4.2 m to 4.4 m, and a wider range of 95% credible interval (2.0–5.8 m; Fig. S3D) when the large natural variability in flow depth is taken into consideration (Figs. S2B and S2D).

Considering the limited information that the method requires from field measurements (i.e., length and thickness of sandstone wedge), the inversion method provides a quick and efficient way to reconstruct the flow depth and channel-bed slope of ancient fluvial-deltaic systems influenced by non-uniform flow. The inversion method thus represents an important quantitative method for paleohydraulic reconstruction of non-uniform flow that is independent of typical paleohydraulic reconstruction methods used for normal flow conditions.

Combining field observations with the inversion method further reveals how surface processes can be inferred based on measurements from rock outcrops and provides the basis for understanding the evolution of ancient fluvial-deltaic landscapes. The current inversion method assumes uniform subsidence; therefore, caution should be taken when applying the method to systems that experienced differential subsidence due to mechanisms such as tectonic tilting (Kim et al., 2010; Dong et al., 2016), differential compaction (Chamberlain et al., 2021), and subsidence-driven rotation (Wu et al., 2019). Moreover, the inversion method should be used to analyze dip-sections along the central depositional axis of a basin.

The present study demonstrates an approach for using variabilities in surface geomorphologic features when reconstructing backwater paleohydraulic conditions from the stratigraphic record (Fig. 15). This approach can be used to inform studies of ancient fluvial-deltaic settings by bolstering assessments of proximity to the marine terminus. This could be useful for studying ancient fluvial-deltaic systems on Earth and Mars, where depositional products of surface geomorphologic features can be readily identified from remote sensing data and rover data collection (DiBiase et al., 2013; Cardenas et al., 2018; Goudge et al., 2018; Wu et al., 2021).

A novel stratigraphic inversion method was developed to reconstruct paleohydraulic conditions for non-uniform flow. The method provides estimates of bankfull flow depth and channel-bed slope based on the first-order geometry of basinward-tapering, fluvial-deltaic sandstone wedge. The method also allows for evaluation of the impact of evolving backwater conditions on stratigraphic architecture over intermediate time scales (10²–10⁴ yr), especially in relation to channel-bed aggradation over the backwater zone.

By linking surface processes and the associated stratigraphic record using the inversion method and field observation, we demonstrated how stratigraphic features can record signatures of backwater hydrodynamics and how this record can be used to improve quantification of paleohydraulic conditions of ancient fluvial-deltaic systems. In particular, we focused on the Tullig Sandstone, Western Irish Namurian Basin, where the inversion method suggests a bed slope of 2.13 × 10⁻⁴ and a normal flow depth of 4.2 m, with 95^th percentile credible intervals of (1.6–3.1) × 10⁻⁴ and 3.6–4.8 m, respectively.

These estimates of paleohydraulic conditions allow interpretation of the downstream variation in sedimentary features that are indicative of non-uniform flow and observable in the Tullig Sandstone: (1) the number of cross-cutting channel stories decreases downstream, along with the total thickness of sandstone of the Tullig Sandstone interval; (2) the percentage of climbing ripples, an indicator of aggradation, peaks around Carrowmore Point; (3) sandstone facies show a reduction in the coarsest fraction of sediment grain size, and an increase in the degree of sorting, progressing downstream; (4) sediment grain size at the top of the Tullig Sandstone fines up-dip, marking the time-diachronous flooding surface that developed during autoretreat; and (5) cross-set and bed thicknesses show similar trends, increasing from Trusklieve, reaching a maximum at Carrowmore Point and Killard respectively, and decreasing downstream. A minimum of the coefficient of variation, for both cross-set and bed thicknesses, is found at Carrowmore Point and Killard respectively, and this is interpreted to reflect elevated channel-bed aggradation.

C. Wu acknowledges support from the American Association of Petroleum Geologists (David Worthington and Meckel Family Named Grants), the Geological Society of America (ExxonMobil/GSA Student Geoscience Grant), the International Association of Sedimentologists (Post-Graduate Research Grant), National Research Foundation of Korea (NRF-2017R1A6A1A07015374), and Yonsei University (Post-Doctoral Researcher Supporting Program, Yonsei University Research Fund #2021-12-0018). We acknowledge support from the University of the Pacific. We thank Simon Chan and Sam Zapp for field assistance. We thank Matthew Carter, who helped with the development of the study and provided field assistance. We thank Jim Best for introducing us to the geology of the Western Irish Namurian Basin. Chris Prince is thanked for providing access to MicroPet. Ordnance Survey of Ireland is thanked for providing access to the historical maps of Ireland. We thank associate editor Alex Whittaker, Sinéad Lyster, and an anonymous reviewer for their detailed reviews and constructive comments.