Reconstructing the morphologies and hydrodynamics of ancient rivers from source to sink: Cretaceous Western Interior Basin, Utah, USA

Quantitative reconstruction of palaeohydrology from fluvial stratigraphy provides sophisticated insights into the response, and relative impact, of tectonic and climatic drivers on ancient fluvial landscapes. Here, field measurements and a suite of quantitative approaches are used to develop a four‐dimensional (space and time) reconstruction of palaeohydrology in Late Cretaceous palaeorivers of central Utah, USA – these rivers drained the Sevier mountains to the Western Interior Seaway. Field data include grain‐size and cross‐set measurements and span five parallel fluvial systems, two of which include up‐dip to down‐dip transects, across seven stratigraphic intervals through the Blackhawk Formation, Castlegate Sandstone and Price River Formation. Reconstructed palaeohydrological parameters include fluvial morphologies (flow depths, palaeoslopes, palaeorelief and planform morphologies) and various hydrodynamic properties (flow velocities, water discharges and sediment transport modes). Results suggest that fluvial morphologies were similar in space and time; median flow depths spanned 2 to 4 m with marginally greater flow depths in southerly systems. Meanwhile palaeoslopes spanned 10−3 to 10−4, decreasing downstream by an order of magnitude. The most prominent spatio‐temporal change is an up to four‐fold increase in palaeoslope at the Blackhawk–Castlegate transition; associated alluvial palaeorelief is tens of metres during Blackhawk deposition and >100 m during Castlegate Sandstone deposition. This study observed no change in unit water discharges at the Blackhawk–Castlegate transition, which argues against a climatically driven increase in palaeoslope and channel steepness. These findings instead point to a tectonically driven palaeoslope increase, although one limitation in this study is uncertainty in palaeochannel widths, which directly influences total water discharges. These reconstructions complement and expand on extensive previous work in this region, which enables the efficacy of quantitative reconstruction tools to be tested. Comparison of results with facies‐based interpretations indicates that quantitative tools work well, but inconsistencies in more complex reconstructions (for example, planform morphologies) highlight the need for further work.

However, palaeohydrology is limited by incomplete (or absent) records of palaeorivers (Sadler, 1981;Jerolmack & Sadler, 2007), uncertainty as to what information fluvial stratigraphy actually preserves (Castelltort & Van Den Driessche, 2003;Jerolmack & Paola, 2010;Romans et al., 2016;Straub et al., 2020), and uncertainties associated with data type, data measurement and reconstruction tools (e.g. Bridge & Tye, 2000). Where it is possible to overcome these challenges, the ability to decipher palaeohydrological information with high fidelity can enable sophisticated insights to be drawn about the sensitivity and response of ancient fluvial systems to tectonic and climatic drivers.
Here, a quantitative framework is used to reconstruct the palaeohydrological evolution of well-known source-to-sink systems of Late Cretaceous central Utah, USA. The focus of this study is the Blackhawk Formation-Castlegate Sandstone-Price River Formation fluvial succession as outcrops are extensive and welldocumented (Kauffman, 1977;Kauffman & Caldwell, 1993;Cobban et al., 2006). These strata represent eastward flowing palaeorivers that drained the Sevier orogenic fold-and-thrust belt to the Western Interior Seaway (WIS). Previous work has primarily focused on qualitative inferences of palaeohydrology in these systems (Miall, 1994;Miall & Arush, 2001;Adams & Bhattacharya, 2005;McLaurin & Steel, 2007;Hampson et al., 2012;Flood & Hampson, 2014), which are sometimes complimented by simple quantitative reconstructions (e.g. Hampson et al., 2013). Meanwhile, quantitative work has mostly focused on architectural-scale elements in these systems, including preservation of channelized bodies and bars and associated autogenic processes, such as avulsion and backwater dynamics (Hajek et al., 2010;Hajek & Wolinsky, 2012;Flood & Hampson, 2015;Trower et al., 2018;Chamberlin & Hajek, 2019;Ganti et al., 2019a). The palaeohydrological evolution of these rivers at the system scale has not been comprehensively addressed using quantitative toolsthis study addresses this outstanding research challenge to shed new light on these ancient systems.
Palaeohydrological field data were collected for five parallel transverse fluvial systems (spaced ca 20 to 25 km apart) across seven stratigraphic intervals within the Campanian stage (83.6 AE 0.2 to 72.1 AE 0.2 Ma) of the Late Cretaceous, which spanned 11.5 Myr (Figs 1 and 2). These data allow for high resolution spatiotemporal reconstructions of these systems, both up-dip to down-dip and along depositional strike (Fig. 1). Reconstructed palaeohydrological parameters include: flow depths; palaeoslopes and palaeorelief (specific to the alluvial domain); hydrodynamic properties, including flow velocities, water discharges and sediment transport modes; and planform morphologies. First and foremost, results show how the morphologies and hydrodynamic properties of these palaeorivers varied in space and time. Moreover, reconstruction of palaeoslopes and palaeorelief in the alluvial domain enable evaluation of the competing roles of tectonic and climatic drivers on the evolution of these ancient rivers. Finally, the results provide new insights regarding the extent to which quantitative palaeohydrological methods (which are increasingly borrowed from the field of engineering) can be reconciled with sedimentological observables. The black outlined box indicates the study area (i.e. part A), and the two highlighted drainage routes (shaded blue) represent the northern and southern depositional-dip transects defined in this study (see part A). (C) The location of Utah relative to the modern North American continent (left) and the Late Cretaceous North American continent (right), which features the Western Interior Seaway (blue). Utah is highlighted as a red box. Jerolmack & Paola, 2010;Hajek & Wolinsky, 2012;Ganti et al., 2013;Ganti et al., 2014;Reesink et al., 2015;Romans et al., 2016;Ganti et al., 2020;Straub et al., 2020) and a number of these studies have focused on Late Cretaceous fluvial strata in central Utah (Flood & Hampson, 2015;Trower et al., 2018;Chamberlin & Hajek, 2019;Ganti et al., 2019a). Meanwhile, other quantitative work has applied fluid and sediment transport models to stratigraphic field data, with an overarching goal of constraining the characteristics of catchments, regional systems or entire fluvial landscapes in the geological past (Ganti et al., 2019b;Lapôtre et al., 2019), or even on other planetary bodies (Lamb et al., 2012;Buhler et al., 2014;Hayden et al., 2019;Lapôtre et al., 2019). This includes using quantitative palaeohydrological tools to reconstruct water and sediment discharges within mass balance frameworks (Holbrook & Wanas, 2014;Lin & Bhattacharya, 2017;Sharma et al., 2017), decipher local palaeogeographies (Bhattacharyya et al., 2015;Li et al., 2018), characterize pre-vegetation rivers (Ganti et al., 2019b), and reconstruct fluvial response to climatic perturbations for well-preserved fluvial strata straddling events such as the Palaeocene-Eocene Thermal Maximum (PETM) (Foreman et al., 2012;Foreman, 2014;Colombera et al., 2017;Chen et al., 2018;Duller et al., 2019).
Despite the breadth of quantitative palaeohydrological tools available, previous applications to fluvial stratigraphic field data have typically centred on individual catchments and instantaneous or short-period intervals (i.e. individual discharge events and mean annual discharges) (Holbrook & Wanas, 2014;Lin & Bhattacharya, 2017;Sharma et al., 2017), or reconstructions across stratigraphic boundaries and short-period tectono-climatic events, such as the PETM (Foreman et al., 2012;Foreman, 2014;Colombera et al., 2017;Chen et al., 2018;Duller et al., 2019). Far fewer studies have focused on long-period intervals, such as the evolution of source-to-sink systems across geological timescales (>10 6 yr). This outstanding opportunity can be exploited in Late Cretaceous fluvial systems of central Utah, where outcrop availability supports a fourdimensional (space and time) study in a region subject to active tectonics, spanning both Sevier and Laramide deformation.
3. An overview of fluvial strata from which palaeohydrological field data were collected. Data were collected for five parallel palaeorivers in Late Cretaceous central Utah, USA. These five palaeorivers cropped out in canyons on the eastern front of the Wasatch Plateau -(A) and (B) show typical exposure of the Blackhawk Formation, Castlegate Sandstone and Price River Formation in these canyons. Specifically, (A) shows strata in Salina Canyon, in which the elevation difference between the road and the top of the cliff is ca 400 m, and (B) shows strata in Straight Canyon, in which the elevation difference between the road and the top of the cliff is ca 450 m (see Fig. 1), and dashed white lines indicate lithostratigraphic boundaries. For two of these five palaeorivers, data were additionally collected upstream to downstream along defined depositional-dip transects (see Fig. 1 (Lawton et al., 2003;Jinnah et al., 2009;Szwarc et al., 2015) (Fig. 1B) and led to down-system sediment mixing (Bartschi et al., 2018;Pettit et al., 2019). Detrital zircon data (Bartschi et al., 2018) indicate that these fluvial systems were dominated by a thrust-belt source in close proximity to the Sevier thrust front, but that more southerly transverse systems may have additionally featured a longitudinal component of drainage (Bartschi et al., 2018;Pettit et al., 2019). Herein, focus is on transverse fluvial systems that predominantly drained the Sevier mountains ( Fig. 1). Tectonic forcing in this region is well-studied (DeCelles, 1994(DeCelles, , 2004DeCelles & Coogan, 2006) and palaeoclimate has been reconstructed from a variety of palaeontological, geochemical-proxy and modelling studies (e.g. Wolfe & Upchurch Jr., 1987;Fricke et al., 2010;Miller et al., 2013;Sewall & Fricke, 2013;Foreman et al., 2015). In central Utah, eastward propagation of the Sevier thrust belt (due to eastward subduction of the Farallon plate) resulted in thin-skinned deformation and movement on the north-south trending Canyon (ca 145 − 110 Ma), Pahvant (ca 110 − 86 Ma), Paxton (86 − 75 Ma) and Gunnison (75 − 65 Ma) thrust systems (DeCelles, 1994(DeCelles, , 2004DeCelles & Coogan, 2006). Associated exhumation created substantial topographic relief in the Sevier mountains, which has been described as 'Andean' in scale with mean elevations approaching near 4000 m (Sewall & Fricke, 2013;Foreman et al., 2015). Modelling results and stable isotope evidence suggest a strong monsoonal precipitation along the eastern flank of the Sevier mountains and seasonal flooding across low-relief regions (Roberts, 2007;Roberts et al., 2008;Fricke et al., 2010;Sewall & Fricke, 2013). The tectono-geographic set-up of the Western Interior was particularly conducive to a monsoonal climatethe proximity of a warm sea to high elevation mountains commonly results in strong seasonal precipitation and convective circulation (e.g. Zhisheng et al., 2001). A seasonal temperate to subtropical climate therefore prevailed throughout Campanian deposition (Parker, 1976b;Kauffman & Caldwell, 1993;Roberts & Kirschbaum, 1995). The Campanian onset of thick-skinned deformation as the subducting Farallon plate transitioned to lower-angle, or flat-slab, subduction (DeCelles, 2004) began to manifest as basement-cored Laramide uplifts (for example, San Rafael Swell, central Utah, and Uinta Mountains, northern Utah), which partitioned the Sevier foreland basin and disrupted patterns of both regional subsidence and drainage (Bartschi et al., 2018;Pettit et al., 2019).

Stratigraphic framework
Establishing a consistent stratigraphic framework in space and time is crucial for system scale palaeohydrological reconstructions. Here, focus is on the Upper Cretaceous Mesaverde Group and up-dip equivalents (Figs 1 and 2) in central Utah, USA, specifically fluvial sediments situated less than ca 100 km from the Sevier orogenic front (DeCelles & Coogan, 2006) in the flexurally subsiding foredeep (Fig. 3). These sediments include the Blackhawk Formation, Castlegate Sandstone and Price River Formation along the eastern front of the Wasatch Plateau (Figs 1 to 3). Up-dip, on the western Wasatch Plateau, the Blackhawk-Castlegate-Price River succession is correlated with the Sixmile Canyon Formation (Indianola Group) and the Price River Conglomerate [following Robinson & Slingerland (1998);Horton et al. (2004); Aschoff & Steel (2011a,b)] (Figs 1 to 3). Up-dip to downdip, these sediments encompass the entire alluvial domain of these palaeorivers draining the Sevier highlands. A broad summary of field sites and the stratigraphic framework (Figs 1 and 2) is given belowextended information regarding regional stratigraphy and correlations is provided in the Supporting Information.
Down-dip field sites were grouped spatially into five field areas that represent five parallel transverse fluvial systems draining the Sevier thrust front: Price Canyon, Wattis Road, Straight Canyon (including Joe's Valley Reservoir), Link Canyon and Salina Canyon (Figs 1 and 3). These five field areas are approximately ca 50 km from up-dip alluvial fan lobes (Figs 1 and 3). Assuming typical outlet spacings of rivers draining orogenic fronts (ca 25 km) (Hovius, 1996), it is likely that these field areas represent five distinct palaeorivers and form a ca 125 km transect along depositional strike. For the two up-dip to down-dip transects (Fig. 1), the northern transect included four field areas: Dry Hollow, Lake Fork, Bear Canyon, and terminating at Price Canyon ( Fig. 3C to E), and the southern transect included three field areas: Mellor Canyon, Sixmile Canyon, and terminating at Straight Canyon ( Fig. 3D to F). These transects follow those widely implemented in previous work, both along-strike (Hampson et al., 2012;Hampson et al., 2013;Flood & Hampson, 2014Chamberlin & Hajek, 2019) and up-dip to downdip (Robinson & Slingerland, 1998;Horton et al., 2004;Aschoff & Steel, 2011a,b).
Down-dip, on the eastern front of the Wasatch Plateau, it is straightforward to assign sediments of the Blackhawk-Castlegate-Price River succession to the appropriate 'space-time' interval by facies associations, following extensive work that has been undertaken in this region (Lawton, 1983(Lawton, , 1986bMiall, 1994; van Wagoner, 1995;Yoshida et al., 1996;Miall & Arush, 2001;Lawton et al., 2003;Adams & Bhattacharya, 2005;Hampson et al., 2012;Hampson et al., 2013;Flood & Hampson, 2014;Hampson et al., 2014;Flood & Hampson, 2015). The lower-middle Campanian Blackhawk Formation represents deposition on coastal plains behind wavedominated deltaic shorelines which, up-section, pass landward into alluvial and fluvial plains (Hampson, 2010;Hampson et al., 2012;Hampson et al., 2013). The size and abundance of channelized fluvial sand bodies (deposited by both single-thread and multi-thread rivers) increase from base to top of the Blackhawk Formation (Adams & Bhattacharya, 2005;Hampson et al., 2012;Hampson et al., 2013;Flood & Hampson, 2015). The middle-upper Campanian Castlegate Sandstone is situated atop the Blackhawk Formation and is an extensive, cliffforming river-dominated deposit. The lower Castlegate Sandstone and upper Castlegate Sandstone (Bluecastle Tongue) comprise amalgamated braided fluvial channel-belt deposits, whereas the middle Castlegate Sandstone comprises less amalgamated, more meandering, fluvial channel-belt deposits with interbedded mudstones (Fouch et al., 1983;Lawton, 1986b;Miall, 1994; van Wagoner, 1995;Yoshida et al., 1996;Miall & Arush, 2001). The ledge-forming upper Campanian Price River Formation sits conformably atop the Castlegate Sandstone and comprises large channelized sand bodies with interbedded siltstones and mudstoneschannelized sand bodies form ca 75% of the formation (Lawton, 1983(Lawton, , 1986b. Fluvial sediments of the Price River Formation represent the end of Sevier thrusting; the late Maastrichtian-Eocene North Horn Formation unconformably overlies the Price River Formation. Up-dip, on the western Wasatch Plateau, correlative strata include more proximal sediments of the Indianola Group and Price River Formation, which is now known to not be timeequivalent with the down-dip Price River Formation exposed near Price, Utah (Robinson & Slingerland, 1998;Horton et al., 2004;Aschoff & Steel, 2011a,b). To avoid confusion, these updip strata are here referred to as the Price River Conglomerate, following Aschoff & Steel (2011a, b). Up-dip to down-dip correlations are limited by incomplete exposure on the western Wasatch Plateau and difficulty in dating conglomerates (see Supporting Information). Nevertheless, Robinson & Slingerland (1998) used palynology to correlate these strata across a variety of localities on the Wasatch Plateau (Fig. 2), which can be traced in seismic reflection data (Horton et al., 2004). The up-dip Price River Conglomerate is time-correlative with the down-dip lower, middle and upper Castlegate Sandstone, and Price River Formation (Robinson & Slingerland, 1998;Horton et al., 2004;Aschoff & Steel, 2011a,b), and is characterized by quartzitedominated synorogenic fanglomerates and few gravel-sand fluvial bodies (Robinson & Slingerland, 1998;Aschoff & Steel, 2011a,b). Of the Indianola Group, the upper Sixmile Canyon Formation is time-correlative with the Blackhawk Formation (Lawton, 1982;Fouch et al., 1983;Lawton, 1986b) and is predominantly characterized by synorogenic gravel-sand fluvial facies, spanning polymictic fluvial conglomerates to medium-coarse-grained sandstones (Lawton, 1982(Lawton, , 1986a. Here a conservative approach is taken to up-dip to down-dip correlations; the upper Sixmile Canyon Formation of the Indianola Group (intervals 1 to 3) is time-averaged, and the Price River Conglomerate (intervals 4 to 7) is also time-averaged, but exceptions were made where field sites were known to be situated at either the top of the upper Sixmile Canyon Formation or at the top/base of the Price River Conglomerate. A full description of these correlations, including new logging in Mellor Canyon, is presented in the Supporting Information.
Each depositional-dip transect is pinned at the most downstream location, i.e. it is assumed that the most down-dip sites in each transect (Price Canyon and Straight Canyon) are approximately parallel and at the same downstream distance. Transects then work upstream, such that the most up-dip field site (Dry Hollow; northern transect) is at a downstream distance of 0 km. Downstream distances follow Robinson & Slingerland (1998) post-depositional extension is not corrected for. Alternatively, when reconstructing along-depositional-strike transects, transects are pinned at the most northern location (Price Canyon) with an along-strike distance of 0 km, meanwhile southern locations have along-depositional-strike distances up to 125 km.

METHODS
Data were collected from channel-fill stratigraphy (cross-stratified sandstone and gravel deposits are interpreted as channel floor deposits) and were time-averaged across each stratigraphic space-time interval (field sites are listed in Table S2 in Supplementary S4). These field data, including uncertainties, were propagated through a quantitative framework to reconstruct the morphologies and hydrodynamics (flow depths, palaeoslopes, river long profiles, flow velocities and discharges, sediment transport modes and likely planform morphologies) of palaeorivers in both space and time.

Field observations
Grain size At each field site the coarse-fraction (>2 mm in diameter) and sand-fraction (<2 mm in diameter) grain sizes of channel-fill deposits were established ( Fig. 4A and B). For coarse-fractions, grain-size distributions were measured via Wolman point counts (Wolman, 1954) (Fig. 4A); this technique has been successfully used to decode spatio-temporal trends in grain size (e.g. Whittaker et al., 2011;D'Arcy et al., 2017;Brooke et al., 2018). For sand-fractions, scaled photographs were processed in ImageJ software and, similarly, the long axis of a minimum of 50 randomly selected grains was measured to recover grain-size distributions (Fig. 4B). From each measured grain-size distribution, the median grain size, D 50 , and 84 th percentile, D 84 , were extracted.
Where grain-size facies were disparate, for example gravel topped with sand, data were collected for each grain-size facies and the proportions of each were estimated (Fig. 4C).
In order to achieve representative sampling for spatio-temporal grain-size trends, multiple grainsize observations were collected at each field site. Not only were data collected for each grain-size facies ( Fig. 4A to C), but depending on overall outcrop extent Wolman point counts were repeated and/or additional scaled photographs were taken for ImageJ processing at intermittent stratigraphic intervals (for example, one count per 5 to 10 m of strata or per channelized body). The extent of each field site can be approximated as the extent of outcrop apparent in Fig. 3C to H. From these data an average grain size was produced for both the sand-fraction and gravelfraction at each field site. Because each spacetime interval includes multiple field sites, this results in multiple average sand-fraction and gravel-fraction grain sizes, capturing channel-fill deposits from several channelized bodies. Finally, a bulk-grain-size was produced for each space-time interval using the gravel to sand proportions at each field siteeach site within a space-time interval was assigned equal weighting. Further information regarding grain-size data collection, including axis selection, sample size sufficiency and weighting, is presented in the Supporting Information.

Cross-sets
Cross-set heights were measured as these data can be used to reconstruct original bedform heights and formative flow depths. Trough cross-bedding and planar cross-bedding, which are inherently indicative of bedload transport, were present at nearly all field sites. They occurred predominantly in sand-grade deposits, but also in granule to pebble-grade deposits ( Fig. 4D to F). To establish mean cross-set heights, the sampling strategy of Ganti et al. (2019b) was followed. Cross-set boundaries (i.e. the lower, asymptotic bounding surface and the upper, erosional bounding surface) were delineated and then heights were measured at regular intervals along the entire width of the cross-set dip-section ( Fig. 4G to I). Measurements were made to a precision of AE5 mm. This protocol was repeated for individual cross-sets within cosets to establish a mean cross-set height for each individual cross-set. Subsequently, maximum cross-set heights (i.e. the maximum distance between lower and upper bounding surfaces) were measured for a representative sample across the exposed outcrop (usually n = 25-50).
From cross-sets for which height distributions were measured (n = 470), the mean, 84 th percentile (P 84 ) and maximum heights of each individual cross-set were extracted. From these data, the relationship between maximum and mean cross-set heights was established. This new relationship was then used to estimate mean cross-set heights from all measured maximum crossset heights (n = 4053). This maximised the amount of field data that could be collected, and therefore analysed, at each field siteit is more efficient to measure maximum cross-set heights than height distributions of individual crosssets. These estimates of mean cross-set heights were propagated through subsequent calculations, as measurements of mean cross-set heights are more appropriate than maximum cross-set heights in reconstruction of palaeohydrological parameters (for example, Eq. 1).

Channel geometry and architectural element data
Above grain and bedform-scales, channel geometries and major architectural elements were also measured, where possible, using a Haglöf Laser Geo laser range finder (Haglöf Sweden AB, Langsle, Sweden) to a precision of AE5 cm. This included maximum channel body/storey thicknesses and bar-scale clinoform heights. Previous work in this region has documented the dimensions and distributions of fluvial architectural elements using high-resolution imagery and 3D outcrop models (Hajek & Heller, 2012;Rittersbacher et al., 2014;Flood & Hampson, 2015;Chamberlin & Hajek, 2019). Field data collection therefore focused on grain-size and cross-set measurements, with compilation of published secondary data (alongside new data from this study) to augment field data and evaluate the palaeohydrological reconstructions herein (see Tables S4 and S5 in Supplementary  S6).

Channel geometries
To calculate original bedform heights from cross-set measurements, the relation of Leclair & Bridge (2001) was used, which is based on theoretical work by Paola & Borgman (1991). Leclair & Bridge (2001) showed that mean bedform (i.e. dune) height, h d , can be approximated as a function of mean cross-set height, h xs , as: where 2.9 is the mean and 0.7 is the standard deviation. Given that exact error margins of palaeohydrological inversion methods cannot be known, a Monte Carlo uncertainty propagation method is used in this study to estimate uncertainty and offer plausible spreads of values for each reconstructed palaeohydrological parameter. In Eq. 1, uncertainty is represented as the mean (μ) and one standard deviation (σ). As such, 10 6 random samples were generated between bounds defined by μ − σ and μ + σ. Samples were generated from a uniform distribution as the shape and the scale of the full distribution of the data is not knownthis approach avoids introduction of additional assumptions. These values are then propagated through subsequent calculations. While bedform height generally scales with flow depth, the mechanistic explanation for this is not fully resolved. As such, many scaling relations simply relate bedform height and flow depth (e.g. Yalin, 1964), whereas some incorporate additional parameters such as Froude number, D 50 , and transport stage (e.g. Gill, 1971;van Rijn, 1984a); however their incorporation does not improve predictive power. Bradley & Venditti (2017) revisited previous bedform height−flow-depth scaling relations and derived a new relation between h d and median formative flow depth, H, based on >380 field observations, where: In detail, Bradley & Venditti (2017) derived two relations to reconstruct H from h d . Their first relation was derived from regression analysis and recovered μ and σ, however the authors argued that this relation is not useful as the data are not normally or log-normally distributed, and that the tails of the distribution are not fully represented (Bradley & Venditti, 2017). The authors additionally presented a non-parametric relation to derive median H (Eq. 2) with a probabilistic uncertainty estimator in which the first and third quartiles of H are given by H = 4.4h d and H = 10.1h d , respectively (Bradley & Venditti, 2017). Bradley & Venditti (2017) noted that this probabilistic uncertainty estimator better represented their data, because it does not assume an underlying distribution. This relation is more appropriate in palaeohydrological reconstructions as, with a larger uncertainty estimate, it offers a broader spread of possible H values. As such, 10 6 uniformly distributed random samples were generated between 4.4 and 10.1, and these model parameter values were used to generate likely palaeoflow depths in these ancient systems. Where cross-bedding was absent (i.e. the most up-dip field sites), channel-body thicknesses were used as a proxy for flow depth.
Similar to H, channel width, W, can be estimated using scaling relations as direct measurement is not normally possible from outcrop. Bridge & Mackey (1993) proposed the relation W = 8.8H 1.82 for single-thread channels. Alternatively, widths of fully-braided channel systems can be approximated as, for example, W = 42H 1.11 (Leopold & Maddock Jr, 1953). However, estimates of W from outcrop data and scaling relations are particularly tentative and, where systems are braided, subject to further uncertainty pertaining to the number of threads. As such, results in this study are reported per unit width.

Palaeoslopes and palaeorelief
Palaeoslopes were estimated using two independent methodologies, adapted from Ganti et al. (2019a). First, Shields stress, τ*, was estimated using the bedform stability diagram of Carling (1999), which expresses bedform stability in terms of τ* and D 50 (for D 50 < 33 mm). Minimum and maximum bounds of τ* for the stable existence of dunes were then identified for a range of D 50 values. Then, 10 6 uniformly distributed random samples of τ* were generated between these grain-size-dependent bounds. Where D 50 exceeded 33 mm, and in the absence of bedforms, possible τ* values of 0.03 to 0.06 were assigned. To reconstruct palaeoslope, S, the bed shear stress, τ b , was approximated as the depth-slope product (τ b = ρgHS) and then S can be given as: where R is the dimensionless submerged specific gravity of sediment in water (1.65 for quartz) and H is the flow depth (ρ is density and g is acceleration due to gravity). For the second approach, the method of Trampush et al. (2014) was used, which is based on Bayesian regression analysis of bankfull measurements in modern alluvial rivers (n = 541); here slope is expressed as: where the constants are given by α 0 = −2.08 AE 0.036, α 1 = 0.254 AE 0.016, and α 2 = −1.09 AE 0.044. Again, 10 6 values of α 0 , α 1 and α 2 were generated (uniformly distributed random samples between μ − σ and μ + σ). Having propagated 10 6 values of τ*, H, α 0 , α 1 and α 2 into these calculations, 10 6 values of S were recovered for both Eqs 3 and 4, which can then be contrasted.
Along up-dip to down-dip transects, palaeoslope estimates can be used to infer the shape of the river long profile, and therefore palaeorelief, in the alluvial domain. Palaeorelief was reconstructed using estimates of S from Eqs 3 and 4. For simplicity, median S was extracted from these values and used to derive palaeorelief. The local slope at downstream position x, S x , can be related to its upstream contributing catchment area, A x , (Hack, 1973;Flint, 1974;Whipple, 2004) as: where k s is the steepness index and θ is the concavity, typically between 0.4 and 0.7 (Tucker & Whipple, 2002). Given that the palaeo-concavity is unknown, a range of plausible concavities (0.4, 0.5 and 0.6) were tested to gauge the spread of possible results. Following Hack's law, local catchment length, where c H is the Hack coefficient, commonly taken as near 2 when L x and A x are in units of km 2 (Castelltort et al., 2009), and h is the Hack exponent, commonly taken as 0.5 (Hack, 1957). Using Hack's law, local slope can instead be estimated as a function of downstream distance, where: k s is calculated from field data using downsystem palaeoslope estimates and knowledge of catchment lengths at each downstream location. Because this study solely focuses on the alluvial domain, this means that up-dip fan apexes would have a catchment length of 0 km. Here, the most up-dip field sites are set as having a catchment length of 5 km to allow for additional up-dip fan length. Knowledge of distance to the coeval palaeoshoreline from the most down-dip sites (Price Canyon and Straight Canyon) is also required. Based on previous studies, approximate distances to the palaeoshoreline are set as ca 10 km for the lower Blackhawk Formation, ca 35 km for the middle Blackhawk Formation, ca 50 km for the upper Blackhawk Formation, ca 110 km for the Castlegate Sandstone (Hampson et al., 2012;Hampson et al., 2013) and ca 200 km for the Price River Formation (Hettinger & Kirschbaum, 2002;Aschoff & Steel, 2011a). A nonlinear least squares regression was used to find best fit palaeoslope profiles (Eq. 6) for both the northern and southern transects at each time interval. Palaeoslope profiles were then transformed into river long profiles by summing elevation increments along the downstream length to the palaeoshoreline. This elevation decrease is indicative of the likely relief in the alluvial domain of these palaeorivers.

Hydrodynamics
In subsequent calculations, values derived from Monte Carlo uncertainty propagation were used, i.e. 10 6 estimates of H, S, etc. Specifically, estimates of S derived from the Shields stress inversion (Eq. 3) were carried forward. Flow velocities, U, were calculated following Manning's Equation, where: and n is Manning's constant, set as 0.03. Water discharges were then estimated by multiplying flow velocity by flow depth, to obtain discharge per unit width (Q = UH).
To determine dominant mode of sediment transport, the Rouse number, Z, was calculated as: where β is a constant that correlates eddy viscosity to eddy diffusivity, typically taken as 1, κ is the von Karman constant, taken as 0.4, and u * is the bed shear velocity (gHS 0.5 ). Sediment settling velocity, w s , was calculated as a function of grain size following Ferguson & Church (2004): where v is the kinematic viscosity of water (1 × 10 −6 m 2 s −1 for water at 20°C) and C 1 = 18 and C 2 = 1 are constants associated with grain sphericity and roundness. With Z, dominant mode of sediment transport is typically wash load for Z < 0.8, 100% suspended load for 0.8 < Z < 1.2, 50% suspended load (i.e. mixed load) for 1.2 < Z < 2.5, and bedload for Z > 2.5. To corroborate inferred sediment transport modes, the particle Reynolds number, Re p , was additionally calculated in line with previous work (cf . Parker, 2004) as: and plotted as a function of τ*, following Dade & Friend (1998). This enables field results to be contrasted with data that are typical of either suspended, mixed or bedload sediments (Leopold & Wolman, 1957;Schumm, 1968;Chitale, 1970;Church & Rood, 1983;Andrews, 1984), and to identify where these data are positioned among characteristic flow regimes (no sediment transport; ripples and dunes; upper plane beds) following Allen (1982a,b).

Fluvial style
Fluvial style (i.e. planform morphology) of Blackhawk-Castlegate rivers has been described qualitatively from outcrop architecture (Miall, 1994;Miall & Arush, 2001;Adams & Bhattacharya, 2005;Hampson et al., 2013). Here, a quantitative approach is implemented to decipher fluvial style to complement these works, check for consistency, and interpret the interplay between different planform morphologies and the tectonogeographic setting. This is carried out for field areas along the eastern Wasatch Plateau. First, Froude number, Fr, is calculated as: and, then, depth/width ratios were plotted against palaeoslope/Froude ratios (Parker, 1976a). Various flow widths were assigned to determine what depth/width ratios are required such that the data fall within the theoretical stability fields for single-thread and multi-thread fluvial planform morphologies. These flow widths are then contrasted with estimates of apparent maximum flow width from architectural analysis of channelized sandstone bodies (e.g. Flood & Hampson, 2015) and field interpretations of fluvial style (Miall, 1994;Miall & Arush, 2001;Adams & Bhattacharya, 2005;Hampson et al., 2013). For all palaeohydrological parameters the median (second quartile) result is presented. Where minima and maxima are presented, these bounds reflect the full spread of recovered values. These are offered as plausible minimum and maximum values for the median, derived from propagation of uncertainty margins. In instances where a first to third interquartile range is additionally presented, specifically in 'box and whisker' plots, this is the first to third interquartile range that has been extracted for each parameter from the 10 6 values recovered by Monte Carlo error propagation. The whiskers in these plots effectively describe the minimum and maximum values of the data, and can also be considered as plausible minimum and maximum values for the median.

Channel geometries
Linear relationships between maximum cross-set height and both the mean and the P 84 cross-set height were established from measured cross-set distributions (n = 470) for the studied field area ( Fig. 5A and B). Maximum and mean cross-set heights are well-correlated (R 2 = 0.88) and 95% of observed mean cross-set heights fall within ca 3 cm of the predicted mean cross-set height. Using these new relationships, mean cross-set heights were estimated for all (n = 4053) measured maximum cross-set heights (Fig. 5C to E; Table S3 in Supplementary S4).
Maximum cross-set heights typically span 0.1 to 0.35 mthese field data are comparable to the results of previous work (e.g. Adams & Bhattacharya, 2005). From maximum cross-set heights, mean cross-set heights spanning 0.07 to 0.25 m are estimated, which correspond with original bedform heights of 0.2 to 0.75 m. Flow depths for the along-depositional-strike transect suggest that, in both space and time, these five transverse fluvial systems maintained median flow depths of 2 to 4 m, with a range of 1 to 7 m (Fig. 6). Overall, flow depths do not change across the Blackhawk-Castlegate transition but exhibit a marginal decrease during middle Castlegate Sandstone deposition of <0.5 m. Flow depths are also projected to be overall <1 m . In (C) to (E), 'n' indicates the number of maximum cross-set heights used to predict mean and P 84 cross-set heights. Full cross-set data for each field site, through each stratigraphic interval, are located in Table S3 in Supplementary S4. greater in southern fluvial systems (Fig. 6). However, these observed differences all lie within the uncertainty margins of calculations, suggesting that these systems were similar to one another.
Reconstructed palaeoflow depths are consistent with independent palaeoflow depth proxies (Table S4 in Supplementary S6), which demonstrates applicability of cross-set scaling relations in the absence of well-preserved macroforms. Bar heights, where available, are consistent with projected flow depths of 2 to 4 m across field sites. For instance, Chamberlin & Hajek (2019) reported mean bar heights of 2.6 m, 3.6 m and 3.9 m for the entire Castlegate Sandstone at Price Canyon, Straight Canyon and Salina Canyon, respectively. At Price Canyon, both Lynds & Hajek (2006) and Hajek & Heller (2012) reported greater mean bar heights of 4.1 m specifically for the lower Castlegate Sandstone, with a typical span of 1 to 8 m (Lynds & Hajek, 2006;McLaurin & Steel, 2007) the authors note that the full range of our reconstructed palaeoflow depths is typically 1 to 7 m and therefore agrees with this range. Meanwhile, channelized fluvial sandstone bodies are more extensively documented for the Blackhawk Formation and their heights offer a maximum limit on palaeoflow depths. Flood & Hampson (2015) recovered mean apparent heights for channelized sandstone bodies of 6 to 8 m across the entire Blackhawk Formation between Straight Canyon and Salina Canyon. As maximum bounds on palaeoflow depth, these values are also in good agreement with the upper bounds of estimated palaeoflow depths.

Palaeoslopes and river long profiles
Palaeoslope estimates for the northern (Fig. 7A to F) and southern (Fig. 7G to M) transects and results from each method (Eqs 3 and 4) were compared (Fig. 7). Palaeoslopes are presented as y/xa palaeoslope of 0.001 results in an elevation decrease of 1 m per 1000 m and is equivalent to 0.057°. Maximum (up-dip) palaeoslopes of 5 × 10 −3 are equivalent to slopes of ca 0.3°; these magnitudes of palaeoslope are comparable with the slopes of modern rivers, including middle-upper reaches of the Colorado (USA) and upper reaches of the Niger (west Africa) (Roberts et al., 2012;Paul et al., 2014;Fernandes et al., 2019). Minimum (down-dip) palaeoslopes of ca 5 × 10 −5 are equivalent to slopes of ca 0.003°; palaeoslopes in the range 10 −5 to 10 −4 are characteristic of lowland/low-slope rivers, such as lower reaches of the Mississippi (USA), Ebro (Spain), Nile (north-east Africa) and Murray--Darling (Australia) (Carlston, 1969;Rudge et al., 2015;Fernandes et al., 2019;Roberts et al., 2019;Soria-Jáuregui et al., 2019).
Up-dip, palaeoslopes are consistently of order 10 −3 (Fig. 7), with the exception of the Blackhawk Formation in the southern transect where first to third interquartile range of recovered palaeoslope values extends down to 7 × 10 −4 (Fig. 7K to M). Importantly, an order of magnitude decrease in palaeoslope is reconstructed between a down-system distance of 10 km and 25 km; this occurs in all stratigraphic intervals, at the same downstream distance, for both the northern and southern transects (Fig. 7). Downdip, from ca 25 km onward, palaeoslopes are flatter and typically span 5 × 10 −5 to 5 × 10 −4 . In these lower gradient regions, there is an apparent down-dip increase in palaeoslope in Fig. 7B, C and I to M. However, this apparent increase is within the first to third interquartile range of values and may not be significant. Updip to down-dip palaeoslope estimates derived from Eqs 3 and 4 are broadly consistent with one anotherthey are the same order of magnitude and the first to third interquartile ranges either overlap with, or are within a factor of 2 to 3 of, one another. However, Eq. 3 overpredicts and underpredicts palaeoslope relative to Eq. 4, such that palaeoslope estimates derived from Eq. 3 imply higher topographic relief and estimates derived from Eq. 4 imply lower topographic relief (Fig. 7).
To constrain temporal changes in palaeoslope, the palaeoslope evolution at the most up-dip locations of both the northern and southern transects can be compared (Fig. 8). Palaeoslopes increase at the onset of Castlegate Sandstone deposition (intervals 4 to 6) and the magnitude of this increase differs between the north and the south (Fig. 8). In the north, the initial palaeoslope is higher (ca 2 × 10 −3 ) and increases by a factor of 1.5 to ca 3 × 10 −3 (Fig. 8A), whereas, in the south, the initial palaeoslope is lower (ca 1 × 10 −3 ) and increases by a factor of up to 4, to ca 4 × 10 −3 (Fig. 8B). This implies a coeval increase in palaeoslope at the onset of Castlegate Sandstone deposition which was more pronounced in the south. Again, estimates derived from Eq. 4 dampen this increase relative to estimates derived from Eq. 3.
With up-dip to down-dip palaeoslope estimates for both the northern and southern transects, best-fit palaeoslope profiles were derived as a function of downstream distance (Eq. 7; Table S6 in Supplementary S7). Palaeoslope profiles generally fit reconstructed palaeoslopes well, with typical R 2 values >0.85, and it is noted that of three reference concavities, θ, used, the higher value of θ = 0.6 typically recovered the best fits (Table S6 in Supplementary S7). A notable exception to this is palaeoslope profiles reconstructed from Shields stress palaeoslope estimates for the Castlegate Sandstone in the northern depositional-dip transectthe lower θ = 0.4 value generates the best fit and this fit is relatively poor (R 2 of 0.35-0.6). However, palaeoslope profiles for these same space-time intervals derived from alternative palaeoslope estimates (Eq. 4) fit well (R 2 > 0.9; Table S6 in Supplementary S7).
In reconstructing palaeoslope profiles steepness index, k s , values were recovered for each stratigraphic interval (for θ = 0.5), which were mostly between ca 5 m and 35 m (Table S6 in Supplementary S7). There is an increase in reconstructed k s values across the Blackhawk--Castlegate transition for both methods of palaeoslope estimation. For estimates derived from Eq. 3, k s values increase across the Blackhawk--Castlegate transition by a factor of ca 2 to 3 in the northern transect, and by a factor of ca 4 to 5 in the southern transect. In contrast, for estimates derived Eq. 4, k s values increase across the Blackhawk-Castlegate transition by a factor of <1.5 in the northern transect, and by a factor of ca 2 in the southern transect (Table S6 in Supplementary S7).
Palaeoslope profiles were transformed into river long profiles, which are indicative of the palaeorelief in the alluvial domain, or depositional reaches, of Blackhawk-Castlegate-Price River fluvial systems only (Fig. 9). Given that the concavities of these ancient rivers are not known, implementing plausible concavities of 0.4, 0.5 and 0.6 enabled a likely spread of values for palaeorelief to be constrained (Fig. 9). Results indicate that different concavities recover similar values for palaeorelief; total estimates vary within a factor of ca 2, between a concavity of 0.4 and 0.6 ( Fig. 9). Using palaeoslope estimates derived from Eq. 3, palaeorelief during Blackhawk deposition was estimated as ca 40 to 60 m in the northern transect ( Fig. 9E and F) and 15 to 25 m in the southern transect (Fig. 9K to M). During Castlegate Sandstone deposition, palaeorelief increased by a factor of 1.5 to 2.5 in the northern transect, to an estimated 65 to 145 m of palaeorelief, whereas it increased by a factor of 5 to 6 in the southern transect, to an estimated 90 to 130 m of palaeorelief. Alternatively, using palaeoslope estimates derived from Eq. 4, palaeorelief during Blackhawk Formation deposition was estimated as ca 30 to 50 m in the northern transect ( Fig. 9E and F) and 15 to 25 m in the southern transect ( Fig. 9K to M). During Castlegate Sandstone deposition, palaeorelief increased by a factor of ca 1.8 in the northern transect, to an estimated 55 to 90 m of palaeorelief, whereas it increases by a factor of 2 in the southern transect, to an estimated 30 to 50 m of palaeorelief. In detail, palaeorelief implied by Eq. 3 (Shields) is up to a factor of 2 greater than the palaeorelief implied by Eq. 4 (Trampush). This higher palaeorelief during Castlegate Sandstone deposition is sustained into Price River Formation times. It is stressed that these estimates refer to the alluvial domain only.

Hydrodynamics and sediment transport
Median flow velocities of 0.8 m s −1 , with a median plausible range of 0.4 to 1.6 m s −1 , are deduced across all field data (Fig. 10A), as well as median unit discharges of 2.5 m 2 s −1 with a median plausible range of 1 to 10 m 2 s −1 (Fig. 10B). Using plausible single-thread channel widths of 100 to 500 m at down-dip locations (see Planform morphologies section), this would imply median total discharges between 250 and 1250 m 3 s −1 , which is comparable with total discharges of well-known North American rivers such as the Platte, Hudson, Colorado, Arkansas and Susquehanna. However, if multi-thread rivers are assumed to possess >1 branch/braid, total discharges would have been several times greater. With a reconstructed increase in palaeoslope at the Blackhawk-Castlegate transition, a coeval increase in flow velocities and unit water discharges is expected analytically. Here, across all up-dip field areas, flow velocities are overall greater during Castlegate Sandstone deposition, up to a factor of 2 to 3 (Fig. 10C), relative to Blackhawk Formation deposition, whereas down-dip flow velocities are broadly the same through time (Fig. 10D). Both up-dip and down-dip, unit water discharges overall do not change at the Blackhawk-Castlegate transition ( Fig. 10E and F). To offer a specific example for the Blackhawk--Castlegate transition (intervals 3 and 4), at Mellor Canyon, median flow velocity, U, increased from 1.9 to 3.0 m s −1 , and median unit water discharge, Q, only increased marginally from 4.4 to 4.6 m 2 s −1 .
Reconstructed Rouse numbers, Z, indicate that dominant transport modes of bed-material varied in space and time (Fig. 11). Up-dip field sites consistently exhibit high Z values for both the median and first to third interquartile range, indicating predominant bedload transport (Fig. 11). Median Z values then decrease by a downstream distance of 30 km, indicating local transition to predominantly mixed load systems, however the likely spread of values indicated by the interquartile ranges implies that dominant transport modes at this downstream distance may have spanned both mixed load and a near entirely suspended load (Fig. 11). A crucial exception to this observation is for Castlegate Sandstone deposition in the southern transect (intervals 4 to 6) where, at a downstream distance of 30 km, median Z values suggest that bedload remains the most important transport mode (Fig. 11G to I). At downstream distances associated with the most down-dip field sites, median Z values have further decreased; however first to third interquartile ranges mostly still span both the mixed load and entirely suspended load domains.
The inferred dominant sediment transport modes are corroborated with results in Fig. 12, in which Shields stress, τ*, is plotted as a function of particle Reynolds number, Re p , for each field site. These data are plotted alongside observed data that are characteristic of suspended load, mixed load and bedload regimes (Leopold & Wolman, 1957;Schumm, 1968;Chitale, 1970;Church & Rood, 1983;Andrews, 1984). Up-dip field sites (Dry Canyon, Lake Fork and Mellor Canyon) plot among secondary data that are typical for bedload rivers, meanwhile all other field sites plot in the mixed-load realm (Fig. 12). Of field sites dominated by a mixed load, data from Sixmile Canyon and Straight Canyon plot closest to the bedload realm, which is consistent with observations in Fig. 11, where results suggest that bedload transport remained important in the southern transect during Castlegate Sandstone deposition (intervals 4 to 6). Overall, results in Fig. 12 suggest that, down-dip, field sites are firmly in the mixed load rangeit is unlikely that bed-material loads were predominantly suspended. In contrast, the first to third interquartile ranges in Fig. 11 suggest that dominant sediment transport modes may have spanned the mixed load/predominantly suspended domain. Downdip, all field sites straddle the bounds between the stability fields for ripples and dunes and upper-stage plane beds (Fig. 12), which implies unidirectional flow and high sediment transport rates (both suspended transport and bedload transport).

Planform morphologies
Finally, these data provide insights into the implied planform morphology of these ancient fluvial systems. However, to do this effectively estimates of palaeochannel widths are needed. Widths are difficult to constrain with confidence from field observations, and estimates from empirical scaling relations are tentative. Assuming single-thread channels, reconstructed median flow depths of 2 to 4 m might suggest channel widths of order 30 to 110 m and, using the upper bound of the 1 to 7 m range, widths up to ca 300 m (following Bridge & Mackey, 1993). In contrast, if multi-thread channel belts are assumed, then channel belt widths of order 90 to 200 m, and up to ca 400 m, might be expected (following Leopold & Maddock Jr, 1953).
For a range of possible widths, palaeoslope/ Froude ratios were plotted against channel depth/ width ratios (cf. Parker, 1976a;Ganti et al., 2019b) (Fig. 13). Results imply that, for Blackhawk--Castlegate-Price River fluvial systems, singlethread planforms would be stable at channel widths <1 km; channel and channel-belt widths >1 km would have been required to instigate formation of bars and support transition to multithread systems, forming vast channel-belt complexes ( Fig. 13A to D). However, planform reconstructions are very dependent on grain size, a factor which is often not evaluated systematically. Bulk grain sizes were used in initial calculations ( Fig. 13A to D; see Methods section). However, when using gravel-fraction grain sizes, which can be associated with tectonic or climatic perturbations (for example, increased palaeoslope or highmagnitude low-frequency discharge events), the results show that multi-thread planforms were more likely (Fig. 13E to H). For gravel-fraction grain sizes, results imply that single-thread planforms were likely stable at channel widths <500 m, and that channel and channel-belt widths >500 m would have supported transition to multithread systems (Fig. 13B).
Of the Blackhawk-Castlegate-Price River fluvial systems, field results for the Castlegate Sandstone plot closest to the singlethread-multi-thread transition, whereas field results for the Price River Formation plot furthest from this transition (Fig. 13). This indicates the relatively high propensity of Castlegate fluvial systems to braiding, relative to Blackhawk and Price River systems.

What did Campanian palaeorivers look like?
These analyses provide new insights that build on previous work characterizing ancient rivers in the Campanian of central Utah as a series of distinct parallel transverse systems draining the Sevier front (Robinson & Slingerland, 1998;Bartschi et al., 2018;Chamberlin & Hajek, 2019;Pettit et al., 2019). These rivers traversed a lowgradient landscape; alluvial relief was tens of metres to ca 100 m, and the length scale of the alluvial domain (i.e. the distance from fan apexes to the palaeoshoreline) varied from as little as ca 70 km during lower Blackhawk Formation deposition, up to and in excess of 250 km during Price River Formation deposition (Hettinger & Kirschbaum, 2002;Aschoff & Steel, 2011a;Hampson et al., 2012;Hampson et al., 2013). Relief was tens of metres during Blackhawk deposition, when the length scale of the alluvial domain was at its narrowest. At the onset of Castlegate Sandstone deposition an increase in palaeoslope is documented, with palaeorelief increasing to ca 100 m, which persisted into Price River deposition (Figs 7 to 9). For comparative purposes, such values of palaeoslope and palaeorelief are characteristic of the Mississippi River and downstream reaches of its principal tributaries, for example the Missouri, Tennessee, Arkansas and Red rivers (Carlston, 1969;Fernandes et al., 2019).
Results imply that palaeoriver morphologies were similar in space and time, with palaeoflow depths of order 2 to 4 m (Fig. 6). Previous detrital zircon results suggest that northerly field sites (Price Canyon and Wattis Road) represent smaller transverse systems and that southerly field sites (Straight Canyon, Link Canyon and Salina Canyon) represent larger systems that include a longitudinal drainage component (Bartschi et al., 2018;Pettit et al., 2019). These results indicate that size disparities between Fig. 12. Shields stress, τ*, plotted as a function of particle Reynold's number, Re p , for all field sites and for each stratigraphic interval (1 to 7), where possible, using bulk grain-size data. Colour-filled circles indicate field results from this study for Bear Canyon (BC), Dry Hollow (DH), Lake Fork (LF), Link Canyon (LC), Mellor Canyon (MC), Price Canyon (PC), Salina Canyon (SC), Sixmile Canyon (SmC), Straight Canyon (StC) and Wattis Road (WR). *For comparison, this plot includes secondary data, originally compiled by Dade & Friend (1998), from Leopold & Wolman (1957); Schumm (1968); Chitale (1970); Church & Rood (1983); Andrews (1984), for characteristic dominant transport modes. Black squares indicate bedload, white circles indicate mixed load and black circles indicate suspended load. Solid black lines indicate stability fields of different flow regimes: no sediment transport (NT), ripples and dunes (R&D) and upperstage plane beds (UP), in line with Allen (1982aAllen ( , 1982b. these five systems were not statistically significantreconstructed variations in palaeoflow depths are within the full range of plausible values. However, palaeoflow depths appear to have been marginally greater in southerly systems (Fig. 6). If true, this may be attributed to the possible longitudinal drainage component (Bartschi et al., 2018;Pettit et al., 2019).
Comparisons with modern rivers suggest that these five parallel palaeorivers (being ca 25 km apart) were substantial systems. Reconstructed hydrodynamic properties, such as flow velocities and unit water discharges, are consistent with the ranges of values of modern systems with similar outlet spacings and similar distances to range fronts (Perry et al., 1996;Schulze et al., 2005; Fig. 13. Theoretical stability fields of fluvial planform morphologies, i.e. single-thread and multi-thread planforms, for both bulk grain sizes (A) to (D) and gravel fraction grain sizes (E) to (H), where present (not all field localities possessed a gravel fraction). For both bulk and gravel grain-size fractions, a range of river widths are assumed (500 m, 1 km, 2 km and 3 km) and used to calculate the depth/width ratio. Data points are for all localities, in space and time, along the defined along-depositional strike transect, i.e. these data points represent the five parallel fluvial systems and do not consider up-dip localities. Data are further subdivided into the Blackhawk Formation (intervals 1 to 3), Castlegate Sandstone (intervals 4 to 6) and Price River Formation (interval 7). Coloured markers indicate the median value, and error bars represent plausible minimum and maximum values for the median, derived from propagated uncertainty margins. Solid black lines indicate the bounds of each stability field, and therefore the predicted transition from single-thread (straight/ meandering) to multi-thread (anabranching/braided) planform morphology. Dashed black lines indicate a potential transition from 1 to 10 threads to >10 threads, based on modern data ). Milliman & Farnsworth, 2013;Global Runoff Data Centre). Notably, unit discharges are overall constant in timethere is no apparent increase in unit discharge at the Blackhawk-Castlegate transition (coeval with palaeoslope increase). This raises questions as to the nature of down-system width evolution and has implications for total discharge plausible single-thread river widths of 100 to 500 m at down-dip locations would imply median total discharges of 250 to 1250 m 3 s −1 . Bedload transport was dominant at graveldominated up-dip localities, as expected, and suspended-load and mixed-load systems prevailed further down-dip, with some localized variations (Figs 11 and 12). For example, results highlight the importance of bedload transport during Castlegate Sandstone deposition in the southern transect (Figs 11 and 12). With this information it is possible to map out how river behaviour varied spatially within catchments, and this informs best practices when it comes to reconstructing sediment discharges. This is especially important when interested in reconstructing the entire sediment load of an ancient system. For instance, channel palaeohydrological approaches are often used to reconstruct sediment discharges in ancient source-to-sink systems (Holbrook & Wanas, 2014;Lin & Bhattacharya, 2017;Sharma et al., 2017); however these reconstruction tools solely reconstruct the bedload fraction and the suspended fraction of the bed material load (van Rijn, 1984b;Wright & Parker, 2004), i.e. the portion of the suspended load that interacts with the bed. As such, these reconstruction tools are not appropriate, by themselves, for reconstructing the total sediment load of a wash load-dominated system, for example. Knowledge of prevailing sediment transport modes is important for evaluating whether different sediment discharge reconstruction methods are consistent with one another, because studies that reconstruct sediment discharges often corroborate results with an independent approach (Lin & Bhattacharya, 2017;Watkins et al., 2018;Zhang et al., 2018;Brewer et al., 2020;Lyster et al., 2020).
Here, reconstructions of planform morphology, following Parker (1976a), and assuming channel widths <1 km, imply that single-thread rivers would have prevailed throughout Blackhawk--Castlegate-Price River deposition. Localized or intermittent transitions to braided planforms may have been associated with tectonic or climatic perturbations, such as increased palaeoslope or high-magnitude, low-frequency discharge events (Fig. 13). In detail, these perturbations (which can be associated with the gravel-fraction grain size) can support braiding at narrower channel/channel-belt widths of order 500 m. Of these fluvial systems, Castlegate systems had a higher propensity to braiding. At this point, it is important to highlight that traditional bipartite classification of fluvial systems aims to define fluvial systems as either straight/ meandering or braided/anabranching end members (Leopold & Wolman, 1957). However, these are not mutually exclusive; both straight/meandering and braided/anabranching planforms can co-exist at reach scales. These reconstructions can be contextualized by field evidence; however, field observations point to a discrepancy and this topic is returned to in the final Discussion section.
To create a holistic view as to the nature of these ancient fluvial landscapes, various modern analogues can be considered. In the Amazon basin, several of the most up-system tributaries axially drain the central and eastern Andean cordillera. For example, the Huallaga River, Peru, is an axial river fed by transverse systems draining the eastern Andean range front. These transverse rivers have regular outlet spacings, channel-belt widths of order hundreds of metres (up to 1 km), and combine both single-thread and multi-thread planforms which vary at reachscales. In the eastern Himalayas, transverse systems draining the range front into the axial Brahmaputra (Assam Valley) provide another modern analogue for the pattern and style of these ancient fluvial systems, despite the larger scale of this system.

What drove spatio-temporal changes in morphological properties?
A key result in this study is quantification of an increase in palaeoslope at the Blackhawk-Castlegate transition by a factor of 1.5 to 4.0, as well as the associated increase in palaeorelief (Figs 7 to 9). Increased palaeoslopes have implications for the morphological and hydrodynamic properties of these palaeorivers, including their flow velocities and unit discharges. In this study, the increase in palaeoslope and palaeorelief implies that rivers were actively responding to changes in uplift rate in the hinterland region.
At the Blackhawk-Castlegate transition, palaeorelief increased from tens of metres to ca 100 m (Fig. 9). An important point to remember is that these estimates are specific to the alluvial domain only. Behind the Sevier front, existence of a highelevation plateau known as 'Nevadaplano' is inferred (Allmendinger, 1992;DeCelles, 1994DeCelles, , 2004DeCelles & Coogan, 2006), which has been likened to the modern high-elevation plateau, Altiplano, of the central Andes. Palaeoelevations in the Sevier highlands and Nevadaplano are argued to be 3 to >4 kmthese values have been deduced from a combination of climate modelling studies (Sewall & Fricke, 2013;Foreman et al., 2015), kinematic reconstructions (DeCelles, 1994(DeCelles, , 2004DeCelles & Coogan, 2006) and other data, including palaeoflora (Chase et al., 1998). Here, alluvial palaeorelief of order 100 m is reconstructed. Given that the low-lying alluvial domain of these palaeorivers has a length scale of order 70 to 250 km, and given the proximity to high-elevation Sevier highlands, the entire river long profile is inferred to have likely been highly concave. This is supported in part by the fact that, in reconstructing palaeoslope profiles, the best fits were recovered when using a higher reference concavity of 0.6 ( Table S6 in Supplementary S7). If best-fit palaeoslope profiles were projected up-dip into the Sevier hinterland, palaeoslopes of 10 −1 might be reached within as little as 10 km of the most up-dip field area, and therefore elevations in excess of 1 km might be reached within a further 10 km. To again use the modern Andes as an analogue, a longitudinal river profile from the Peruvian shoreline to the western Andean cordillera and Altiplano would have a length scale of 50 to 150 km, with 0.5 to 1.0 km of relief in the alluvial domain and elevations >3 km in the western cordillera and Altiplano. With a similar tectonogeographic setting in Late Cretaceous Utah, this comparison can also be used to highlight the potential high concavity of these ancient river profiles.
In reconstructing palaeorelief, steepness indexes, k s , were also recovered for northern and southern transects (Eqs 5 and 6) ( Table S6 in Supplementary S7). While k s was solved for using field data and a nonlinear least squares regression, k s values are often estimated (albeit tenuously) as a function of known uplift rate and erodibility in bedrock channels, but additionally (although less frequently) in downstream alluvial reaches (Kirby & Whipple, 2012;Pederson & Tressler, 2012;Stucky de Quay et al., 2019). Inversely, where k s can be measured, and where erodibility is known, first-order estimates of uplift rate can be made. Steepness indexes recovered in this study were typically ca 5 to 35 m (for a reference concavity, θ, of 0.5) and, despite unknown erodibility, global data compilations indicate that low uplift rates of order 0.01 to 0.1 mm yr −1 are generally associated with these kinds of values (Kirby & Whipple, 2012). Despite overall low k s values, it is important to note the relative increase in k s by a factor of <1.5 to 5.0 at the Blackhawk−Castlegate transition. While these are first-order estimates, and are derived solely for the alluvial domain, an increase in k s (and palaeorelief) can be attributed to a relative increase in uplift rate in the hinterland region. Here, this increase might be attributed to frontal thrust migration, or thrust initiation in the Sevier highlands (DeCelles, 2004;DeCelles & Coogan, 2006). This includes Sevier shortening in the Charleston-Nebo Salient (CNS), an eastward convex portion of the Sevier thrust front in northcentral Utah (Fig. 1B) (Bruhn et al., 1986;Bryant & Nichols, 1988;Constenius et al., 2003;Bartschi et al., 2018), which is commonly attributed to the influx of quartzite-dominated coarse-grained detritus associated with Castlegate Sandstone progradation (Robinson & Slingerland, 1998;Horton et al., 2004). For Castlegate Sandstone deposition in the northern transect, results show that palaeoslope profiles did not fit reconstructed palaeoslopes well and favoured lower concavities (which also did not fit well). The interpretation herein is that shortening in the CNS, which has been structurally linked with coeval basement Laramide uplifts in northern Utah (Bruhn et al., 1986;Bryant & Nichols, 1988;Constenius et al., 2003;Bartschi et al., 2018), may have significantly influenced river long profiles associated with northerly Castlegate fluvial systems near Price, and locally lowered their concavities. Whereas ca 60 km south in the southern transect, higher concavity values of 0.6 deliver best fitting palaeoslope profiles through all seven stratigraphic intervals (Table S6 in Supplementary S7).
While tectonic drivers are commonly attributed to variations in channel steepness (Kirby & Whipple, 2001;Kirby et al., 2003;Wobus et al., 2006;Boulton & Whittaker, 2009;DiBiase et al., 2010), climatic drivers, especially precipitation rates, also play a crucial role but are notoriously difficult to disentangle from their tectonic counterpart (Wobus et al., 2010;DiBiase & Whipple, 2011;Champagnac et al., 2012;Whittaker, 2012;D'Arcy & Whittaker, 2014). The role of climate is important to consider here, given the assumed monsoonal climate and, therefore, highly seasonal discharge variability (Roberts, 2007;Roberts et al., 2008;Fricke et al., 2010;Sewall & Fricke, 2013). Previous work shows that precipitation rates have a discernible role on steepness indexes D'Arcy & Whittaker, 2014); analytically, an increase in channel steepness and palaeoslope can be attributed to a decrease in precipitation rate (to maintain similar total water discharge) (D'Arcy & Whittaker, 2014). To reduce palaeoslopes by a factor of 2, precipitation rate must typically be quadrupled (D'Arcy & Whittaker, 2014). Despite the supposed warm and wet climate (Parker, 1976b;Kauffman & Caldwell, 1993;Roberts & Kirschbaum, 1995), few workers have argued for, or investigated, the possibility of increased aridity at the Blackhawk-Castlegate transition ( van Wagoner, 1995;Adams & Bhattacharya, 2005). In theory, increased palaeoslopes can be explained by decreased precipitation (D'Arcy & Whittaker, 2014); however, here, no decrease in either flow velocities or unit discharges is reconstructed at the Blackhawk-Castlegate transition (Fig. 10). Generally, in down-dip locations, flow velocities and unit discharges are constant across this interval ( Fig. 10D and F). At up-dip field sites, however, flow velocities are overall slightly greater during Castlegate Sandstone deposition relative to Blackhawk Formation deposition, but unit discharges remain similar for both.
With unit discharges constant in space and time, the crucial unknown is palaeochannel width. At minimum, channel widths can be considered as broadly the same across the Blackhawk-Castlegate transition. During Blackhawk Formation deposition, channelized sandbody widths of order 350 to 420 m offer a maximum limit on palaeochannel widths (Hampson et al., 2013;Flood & Hampson, 2015). Meanwhile, during Castlegate Sandstone deposition, bar package widths are between ca 60 to 180 m (Chamberlin & Hajek, 2019); assuming two or three threads, these bar widths might imply channel belt widths of order half a kilometre. However, planform stability estimates based on Parker (1976a) indicate that these rivers could have possessed anywhere between one and ten threads (Fig. 13), which could result in channel-belt widths up to and in excess of 1 km. At maximum, this implies increased channel widths at the Blackhawk-Castlegate transition. Unless a significant decline in river widths is projected, then field results do not directly support a climatic driver. Consequently, the interpretation here is that increased channel steepness and palaeoslope at the Blackhawk-Castlegate transition is due to tectonically driven uplift in hinterland regions.

Effectiveness of palaeohydrological and palaeomorphological reconstructions
While quantitative reconstructions have led to significant advances in both the quantity and level of detailed information that can be extracted from fluvial strata (e.g. Ganti et al., 2019a), it is unclear how accurately these tools characterize ancient systems. Addressing this question is particularly important as sedimentology becomes increasingly numerical and it becomes easier to apply quantitative tools to stratigraphy (Duller et al., 2010;Whittaker et al., 2011;Holbrook & Wanas, 2014;Ganti et al., 2019b). With extensive existing work on Late Cretaceous fluvial systems of central Utah, results in this study offer a unique opportunity to highlight consistencies and discrepancies between quantitative interpretations of fluvial palaeohydrology and more qualitative fieldbased facies and architectural interpretations.
To first order, whether point reconstructions of various morphological and hydrodynamic parameters agree with qualitative interpretations can be evaluated using independent proxies (derived from field measurements or facies interpretations). As previously mentioned, reconstructed flow depths agree with several secondary observations of bar heights (Adams & Bhattacharya, 2005;Lynds & Hajek, 2006;McLaurin & Steel, 2007;Hajek & Heller, 2012;Chamberlin & Hajek, 2019) (Table S4 in Supplementary S6), which can be used as a direct proxy for flow depth (Bridge & Tye, 2000;Hajek & Heller, 2012). This agreement indicates that the uncertainty estimator in Eq. 2 is reasonable, and that cross-set heights can therefore be used to reconstruct reasonable flow-depth constraints and are useful as a bedform-scale approach. Such an approach is particularly useful in core data, locations with limited outcrop exposure, or deposits where the degree of bar preservation is poor. It is noted that scaling relations that relate cross-set heights with original bedform heights (and subsequently formative flow depths) are derived from theory and experiments that assume statistical steady state, in which flow is constant (Paola & Borgman, 1991;Leclair, 2002;Jerolmack & Mohrig, 2005). As such, agreement of flow depth reconstructions with bar heights might therefore imply that these dunes were formed in steady flow conditions . This contrasts with literature that alludes to the preferential preservation of dunes in unsteady flow conditions (Reesink & Bridge, 2007;Reesink & Bridge, 2009;Reesink et al., 2015;, and merits further work regarding the kinematic controls on dune preservation in this region. For more complex palaeohydrological reconstructions, such as palaeoslopes and palaeorelief (Figs 7 to 9), it is not possible to directly corroborate estimates with independent proxies derived from field data. Nevertheless, it is still possible to evaluate reconstruction tools by contrasting commonly used methods. In this study the first approach used a theoretically-based Shields stress inversion (Eq. 3), whereas the second approach used the empirically-derived model (Eq. 4) of Trampush et al. (2014). Palaeoslope estimates derived from each approach are in broad agreement with one another. Each method typically recovers estimates of the same order of magnitudein many cases the interquartile ranges of values overlap, and, in all cases, the full ranges of plausible values overlap (i.e. the whiskers in Figs 7 and 8). These point comparisons between the two methods are promising, and in line with comparisons made elsewhere (e.g. Ganti et al., 2019a). However, there are implications when larger spatial scales are concerned, imparting uncertainty that must be carried forward in interpretation of palaeorelief in the depositional reaches of these systems. Along the northern and southern transects, Shields stress inversion estimates consistently show higher differences in palaeoslope (i.e. higher slopes up-dip and lower slopes downdip) relative to palaeoslopes derived from the Trampush et al. (2014) method. This difference is likely an outcome of the Trampush et al.
(2014) method using a continuous function to estimate slope, whereas the Shields stress inversion relies on a step-change empirical estimate for gravel or sand-bed rivers. Regardless of the method used, palaeoslope reconstructions are dependent on grain-size and flow-depth estimates. Because flow depths did not appreciably change in Blackhawk and Castlegate palaeorivers, variations in reconstructed slopes and derivative estimates (for example, water and sediment discharge) are largely driven by observed differences in grain size.
Despite the differences of the two methodologies on palaeorelief, estimates of palaeorelief can be compared with relief in modern systems possessing similar tectono-geographic set-ups. Palaeorelief estimates between 50 m and 100 m in depositional reaches of these ancient fluvial systems are reasonable when compared with relief in modern systems with a similar tectonogeographic setting. For example, one can return to the Andean analogue, but cross over to the eastern Andean cordillera and into the foreland basin and low-lying plains of the Amazon River. For most of its course, the Amazon long profile has a relief of less than 100 m (Milliman & Farnsworth, 2013)relief only exceeds 100 m in proximity to the range front (Milliman & Farnsworth, 2013).
Finally, these results complement field evaluation of the nature of Blackhawk Formation and Castlegate Sandstone planforms, but also raise new questions. Channelized sandstone bodies of the Blackhawk Formation are typically 350 to 420 m wide (Adams & Bhattacharya, 2005;Hampson et al., 2013;Flood & Hampson, 2015), although a small proportion are much larger and some exceed 1 km (Flood & Hampson, 2015). These sandstone bodies offer a maximum cap on palaeoflow width. The Blackhawk Formation is considered to mostly represent single-thread systems, which results in this study agree with. However there is significant field evidence that many channelized sandstone bodies of the Blackhawk Formation represent multi-thread systems with mid-channel bars, based on bar facies observations (Adams & Bhattacharya, 2005;Hampson et al., 2013;Flood & Hampson, 2015). Field observations of multi-thread Blackhawk fluvial systems of order hundreds of metres are inconsistent with the present results, which suggest that multithread systems would not have been stable (Fig. 13). Meanwhile, the Castlegate Sandstone is interpreted to be fully-braided from facies observations (Miall, 1993(Miall, , 1994Miall & Arush, 2001;McLaurin & Steel, 2007). Reported mean bar package widths of order 60 to 180 m for the Castlegate Sandstone (Chamberlin & Hajek, 2019) would imply total channel widths <1 km (assuming a few braids); this study's reconstructed planform stability estimates, which indicate that Castlegate systems should have been singlethreaded, are again inconsistent with sedimentological facies and architectural interpretations. Other quantitative reconstructions of planform have contradicted traditional field-based facies observations (Ganti et al., 2019a), and these inconsistencies must be treated carefully. The main limitation to reconstructing ancient channel planforms is a lack of reliable methods for estimating palaeochannel widths. Interpreting palaeochannel planforms from facies associations and stratigraphic-architectural data is not trivial, particularly where outcrop is limited or where observations are equivocal. However, in this case, a number of workers have concluded that braided conditions prevailed at the time of Castlegate Sandstone deposition (Lawton, 1986b;Miall, 1994; van Wagoner, 1995;Miall & Arush, 2001) and occurred at times during Blackhawk Formation deposition (Adams & Bhattacharya, 2005;Hampson et al., 2013;Flood & Hampson, 2014. As such, it can be argued that further detailed work to test and reconcile facies-based and hydraulically derived interpretations of channel planforms is a pressing research goal.

CONCLUSIONS
Here a four-dimensional reconstruction of palaeohydrology in Late Cretaceous palaeorivers of central Utah, USA, is presented, using field data and a well-established quantitative framework. Overall, fluvial morphologies were similar in space and time, although marginally greater reconstructions of flow depths in southerly systems likely reflect the contribution of a longitudinal drainage component. The most prominent spatio-temporal change is an increase in palaeoslope at the Blackhawk-Castlegate transition by a factor of 1.5 to 4.0; this reflects an increase in palaeorelief (for the alluvial domain) from tens of metres during Blackhawk Formation deposition up to, and in excess of, 100 m during Castlegate Sandstone deposition, which persisted into Price River Formation times. The observation that unit water discharges do not change at the Blackhawk--Castlegate transition does not support a climatically driven increase in palaeoslope and channel steepness. Results therefore point to a tectonically driven palaeoslope increase. In deciphering the relative role of tectonic and climatic drivers, the main limitation in this study is uncertainty in palaeochannel widths, which directly affect total water discharges. Palaeochannel width reconstructions therefore remain a prominent research challenge.
Results complement and expand on extensive facies-based interpretations of these systems, which offers a unique opportunity to evaluate the efficacy of quantitative palaeohydrological reconstruction tools. Bedform-scale palaeoflow depth reconstructions are in good agreement with observations of preserved barforms. Moreover, while different palaeoslope reconstruction methods produce results that broadly agree, the results show that at larger spatial scales they over-predict and under-predict relief relative to one another, which has implications for quantifying alluvial palaeorelief and, therefore, the magnitude of change in relief at the Blackhawk-Castlegate transition. Finally, quantitative hydraulic reconstructions of planform somewhat disagree with facies-based interpretations. While this discrepancy ties back to uncertainty in palaeochannel widths, these results highlight that further work is required to reconcile hydraulically-based and facies-based approaches in order to facilitate their application in the geological past.

Supporting Information
Additional information may be found in the online version of this article: Supplementary S1. Variables list. Supplementary S2. Field localities. Fig S1. Study area showing key localities mentioned in the supplement, reproduced from Fig. 1 in the main manuscript. Table S1. Field localities visited in this study, for each field area (e.g. Price Canyon, Wattis Road, etc). Supplementary S3. Regional correlation -further information. Fig S2. Measured section through the Sixmile Canyon Formation (Indianola Group) and (extrapolated) Price River Conglomerate at Mellor Canyon. Supplementary S4. Field data. Table S2. Grain-size data collected and used in this study. Table. S3. Cross-set data collected and used in this study.
Supplementary S5. Grain size sample sufficiency. Fig S3. The iterative convergence of median grainsize for (A) pebbles-cobbles, (B) medium-coarse pebbles, (C) granules-fine pebbles, (D) medium-coarse sand, and (E) fine-medium sand, as calculated from scaled photographs in ImageJ software. Fig S4. The iterative convergence of median grainsize for different outcrops of gravel-grade sediments, as calculated from field Wolman counts. Supplementary S6. Secondary field data. Table S4. A compilation of field measurements (secondary data from published literature) for architectural scale elements, e.g. bar heights, that are commonly used as palaeoflow depth proxies. Table S5. A compilation of field measurements (secondary data from published literature) for architectural scale elements, e.g. sandstone bodies, that are commonly used as a proxy to infer the magnitude of channel width. Supplementary S7. Goodness of fits on palaeoslope profiles including resolved steepness indexes. Table S6. Steepness indexes, k s , recovered for the defined northern and southern depositional-dip transects, through each stratigraphic interval (1-7), where possible. Supplementary S8. Additional results. Fig S5. Reconstructed palaeoflow depths for the five parallel fluvial systems, for each stratigraphic interval where possible.