Down‐delta hydraulic geometry and its application to the rock record

Palaeodischarge estimation is largely undertaken within fluvial settings. However, there are limited palaeodischarge estimates specifically from delta deposits, despite their significance globally. Estimating water palaeodischarges for deltas using catchment‐based approaches developed using data from fluvial settings requires estimating parameters from the rock record (for example, palaeotemperature, palaeoslope and palaeorelief). These may be difficult to determine, leading to under‐estimation or over‐estimation of palaeodischarge values due to differences in process‐form relationships between alluvial rivers and deltas. When a sediment‐conveying fluvial channel enters a standing body of water, delta lobes develop through repeating mouth bar deposition due to flow deceleration, forming a deltaic morphology with distributary channel networks that differ morphologically from those developed in unidirectional flowing alluvial rivers. This study provides empirical relationships determined across five climate regions, using 3823 measurements of distributary channel width from 66 river deltas alongside the trunk river bankfull discharge that feeds into the entire delta, using a hydraulic geometry scaling approach. Empirical relationships are developed from the global delta dataset between bankfull discharge and catchment area (Qb–A), and bankfull discharge and median distributary channel width (Qb–Wmed). These empirical relationships produce very strong statistical correlations, especially between Qb and Wmed, across different climate regions (Qb = 0.34 Wmed1.48, R2 = 0.77). However, both Qb–A and Qb–Wmed relationships have outliers that may be explained by particular hydrological or geomorphic conditions. These new empirical relationships derived from modern systems are then applied to Cretaceous outcrops (Ferron Sandstone and Dunvegan Formation). The comparatively simple scaling relationships derived here produced palaeodischarge estimates within the same order of magnitude as palaeodischarge values previously obtained using existing, more complex approaches. This study contributes to source‐to‐sink investigations by enabling palaeodischarge estimates that intrinsically account for climate impacts on channel geometry at the time of deposition, using measurements of channel width or catchment area of a deltaic outcrop.

Delta channels have different flow conditions from the unidirectional flow of fluvial channels due to bidirectional flow in some parts of deltas and widespread influence of the backwater affect generated by the elevation of the receiving water body. The presence of a backwater zone decelerates and sometimes reverses the unidirectional flow of water influx from the alluvial rivers (Nittrouer, 2013;Wu & Nittrouer, 2019). When a sediment-conveying fluvial channel starts to debouch into a standing body of water, delta lobes initially develop through mouth bar deposition. Multiple successive mouth bars accumulated in front of a river mouth form the characteristic distributary deltaic morphology with channel networks that merge upstream at the delta apex (Edmonds & Slingerland, 2007). The split between two or more newly formed terminal distributary channels occurs due to this mouth bar deposition (Wright, 1977;Edmonds & Slingerland, 2007;Kleinhans et al., 2013). Hence, any estimate of flow discharge from a deltaic system will be incomplete if only a single distributary channel is considered. In addition, open water deltas are typically exposed to marine processes in the form of tides, wave energy, storm surges, lake level and sea-level change (Galloway, 1975;Wright, 1977;Hoitink et al., 2017). Due to these marine processes, the unidirectional fluvial channel flow becomes less prominent particularly closer to the shoreline .
In the proximal parts of deltas, channels are fluvially controlled and fluvial morphometric relationships apply; once the channel is influenced by marine processes, this scaling may change. In the distal part of a delta, the presence of large tidal, wave energy or backwatercontrolled flow regimes will significantly alter the geometry of delta distributary channels both in modern systems and in the rock record (Chatanantavet et al., 2012;Lamb et al., 2012;Nittrouer, 2013;Fernandes et al., 2016;Ganti et al., 2016;Martin et al., 2018;Chadwick et al., 2019Chadwick et al., , 2020. There is a need to understand morphometric scaling relationships in deltaic systems as these are likely to differ from fluvial systems because of their distributive pattern and marine processes that directly influence delta distributary channel geometry. This study uses globally-available satellite imagery, catchment area and river discharge datasets to build empirical scaling relationships between delta distributary channel widths, river catchment area and water discharge. The aims of this paper are thus to: (i) investigate the relationship between bankfull river discharge (Q b ) and catchment area (A) across different climates (Q b -A relationship), for modern river deltas; (ii) investigate the relationship between bankfull discharge and delta distributary channel widths (W) across different climates (Q b -W relationship); and (iii) outline and discuss how these relationships may be employed in deep-time stratigraphic successions, particularly where proxies for palaeoclimate can be retrieved. Such empirical relationships enable palaeogeographical reconstruction using modern systems, particularly where based on geometric properties that can be easily extracted from geological deposits. Empirical scaling relationships could aid our global understanding of delta hydraulic geometry, both for modern and ancient systems (Mikhailov, 1970;Andrén, 1994;Edmonds & Slingerland, 2007;Sassi et al., 2012;Gleason, 2015). equations relating, for example, the catchment erodibility (B), water discharge (Q), area (A), relief (R) and annual temperature (T) (the 'BQART' model - Syvitski & Milliman, 2007). Some of these approaches require measurements or estimates of parameters that are commonly challenging to obtain from rock record datasets (for example, palaeotemperature, relief, palaeoslope and catchment area) (Syvitski & Saito, 2007;Davidson & North, 2009;Holbrook & Wanas, 2014;Brewer et al., 2020). All available methods make assumptions, for example when using geometric scaling the channel geometry is assumed to be in equilibrium with the bankfull water discharge.
One of the most commonly used models, the Fulcrum model, assumes dynamic equilibrium where all sediment mass transported through a trunk channel is balanced by sediment mass eroded upstream and deposited downstream (Holbrook & Wanas, 2014). This model also assumes a fixed position and dimension of a rectangular palaeochannel geometry. Values of dimensionless bankfull Shields' stress and the Chezy friction coefficient are assumed, from which palaeoslope, velocity and bankfull depth, hence palaeodischarge, are calculated (Brewer et al., 2020;Lyster et al., 2021). The Shields' stress is challenging to estimate from ancient deposits, although it can be constrained using information on, for example, grain-size distribution (Ganti et al., 2019).
The second widely applied model for estimating palaeodischarge is the 'BQART' model, which utilizes catchment-scale parameters. Although the original goal of this model was to estimate the total suspended solid load (TSS) brought by the fluvial system to the ocean, it can be used to estimate discharge or palaeodischarge and is applicable to ancient sedimentary systems (Blum & Hattier-Womack, 2009;Sømme et al., 2011;Allen et al., 2013;Watkins et al., 2019). The 'BQART' model parameters can often be only partially constrained. For example, estimating palaeotemperature relies on proxy information (for example, biomes of flora and fauna, palaeosols, mineralogy) combined with plate tectonic reconstructions, which increase the uncertainty in 'BQART' sediment load estimates, especially in cooler climates (Nyberg et al., 2021).
Scaling between discharge and channel width and depth is an inevitable consequence of channel size adjusting to the volume of water being conveyed. Hydraulic geometry provides a theoretical basis for such scaling. Hydraulic geometry refers to empirical relationships relating channel width (W), depth (d) and velocity (v) to discharge (Q) (Leopold & Maddock, 1953). As discharge fluctuates at a single site, strong power relationships of the following form are found: with the coefficients (a, c, k) and exponents (b, f, m) derived empirically from repeat measurements (Leopold & Maddock, 1953). From the continuity equation The values of b, f and m are constrained by the hydraulics of water flow (Ferguson, 1986). For a discharge of specified recurrence interval, such as bankfull discharge, consistent downstream hydraulic geometry relationships exist, taking the same form as Eqs 1 to 3. In distributary deltas, the downstream relationships reflect abrupt reductions in discharge at bifurcations and also the increasing influence of bidirectional flow towards the downstream margin of the delta. Hence, 'down-delta' hydraulic geometry is complex but at any location along a distributary channel Eqs 1 to 3 applies consistently due to the continuity of discharge. This study investigates empirical relationships from 66 catchments feeding river deltas across different climate regions, that include 3823 distributary channel width measurements from 66 river deltas available at https://doi.org/10.6084/m9. figshare.19574938.v2 (Fig. 1A). Catchment areas and their associated bankfull discharges were related to the median channel width of each delta. As data for 66 rivers are available [via the Global Runoff Data Center (GRDC)] at the location close to the delta apex and not at individual distributary channels, this implies that the bankfull dischargemedian channel width relationship allows prediction of total system discharge/palaeodischarge that contributes to the whole delta plain, and not to individual distributary channels. The median is chosen for three reasons: firstly, it provides a more conservative estimate of central tendency than the mean in cases where there may be very wide channels close to the downstream limit of the delta; secondly, the preservation potential of delta channel deposits is greater away from the downstream limit and the median thus better represents channels that are likely to be preserved (Olariu & Bhattacharya, 2006); and the influence of outliers in ancient deposits is reduced by using the median. Moreover, the order of distributary channels, or number of channel bifurcations in the delta plain is not incorporated to simplify the scaling relationships built in this study. Assuming that the measured distributary channel widths are approximately bankfull widths, scaling relationships are determined between the measured median distributary channel widths and Q 2 (two-year recurrence flood as an estimate of bankfull discharge, Q b ) in the river, and between catchment area and Q 2 (Leopold & Maddock, 1953;Gleason, 2015). A re-arrangement of Eq. 1, Q b = αw β (Q b -W relationship) is used because this provides a basis for sedimentologists to estimate bankfull discharge from channel widths, measurement of which is often achievable in ancient deposits.

METHODS
Empirical statistical relationships were found between the median widths of delta distributary channels gathered from satellite imagery and their site-specific discharges. Although backwater effects in the form of wave and tidal influences may be present, other studies have demonstrated the effectiveness of this relationship in deltaic environments (Mikhailov, 1970;Andrén, 1994;Edmonds & Slingerland, 2007;Sassi et al., 2012;Gleason, 2015). Bankfull discharge has widely been considered as the flow that controls channel geometry in alluvial rivers (De Rose et al., 2008;Haucke & Clancy, 2011;Gleason, 2015), and is estimated here as Q 2 , where 2 is the recurrence interval (years) of the discharge, as also widely used by others (Eaton, 2013;Jacobsen & Burr, 2016;Morgan & Craddock, 2019). Distributary channel widths on the 66 river deltas were measured in ArcGIS software using annual composite Landsat 5 satellite images. Delta apex (i.e. valley exit) locations were obtained from digital elevation models (DEM) from the Shuttle Radar Topography Mission (SRTM) and Arctic-DEM (Tucker et al., 2004;Farr et al., 2007;Morin et al., 2016) (Fig. 1A). Satellite imagery from 1984 was used where available, with some imagery dated to more recent years. Using the older (1984) images reduces the impact of infrastructure and bank protection on channel widths. The satellite images and DEMs were projected using World Geodetic System (WGS) 1984 in ArcGIS to measure the channel widths and to extract valley exit locations (Hartley et al., 2017).
River deltas were identified based on their protrusion beyond the original lateral shoreline (Caldwell et al., 2019). Criteria for selecting river deltas includes any channel mouth that intersects with the open seawater, depositing sediment that protrudes beyond their lateral shoreline. Nonetheless, the river deltas studied here are not classified based on their dominant forces (for example, wave-dominated, tidedominated or river-dominated deltas) due to delta morphodynamics varying in time and space (for example, a tide-dominated delta could transform into a river-dominated delta or a wave-dominated delta into a river-dominated delta) and very few delta end-members exist in nature (Syvitski & Saito, 2007). The authors also note that some influence of tide and wave processes may exist in the dataset (Ta et al., 2002;Correggiari et al., 2005;Syvitski & Saito, 2007). However, as this paper focuses on the estimation of river discharge from distributary channel morphology, river deltas with clear wave and tidal morphologies were avoided (for example, abundant tidal creeks, deflected delta distributaries, elongated/parallel shoreline).
Channel widths were measured using a method, adapted from Sassi et al. (2012), in which a semicircular grid s/L is used to define a dimensionless distance from the delta apex to the shoreline, where s represents channelized distance from the delta apex and L is the channelized distance along the longest distributary channel (Fig. 1B). This grid allows measurement of the widths of multiple distributary channels located at the same dimensionless distance from the apex, hence allowing comparison across differently sized deltas. The apexes were defined as the valley exit points as recognized on DEMs (Hartley et al., 2017) or as the most landward avulsion node within the delta (Ganti et al., 2016). The semicircular grid has a resolution of about 10 times the width of the river channel at the first avulsion point to maintain consistent dimensionless distance and data frequency across deltas of varying size. As an example, the Mahakam delta, Indonesia, has a 500 m wide channel at the avulsion point which is ca 40 000 m following the longest channel from the shoreline (L). Channel widths are measured every 5000 m from the delta apex (i.e. s/L = 0) to the delta shoreline where s/L = 1 (Fig. 1C). Widths of distributary channels were included, and tidal creeks were omitted.
Catchment areas were delineated in ArcGIS using the watershed polygons available from the HydroBASINS dataset (Lehner & Grill, 2013). River discharge data for the closest measuring location to the delta apex were extracted from the Global Runoff Data Centre (GRDC) dataset (https://www.bafg. de/GRDC/EN/Home/homepage_node.html). The two-year recurrence interval flood (Q 2 ) was used to estimate the bankfull discharge, or the dominant channel-forming flow, and is referred to as discharge (Q) subsequently for simplification (Wolman & Miller, 1960;Phillips & Jerolmack, 2016Edwards et al., 2019;Dunne & Jerolmack, 2020;Rhoads, 2020). Q 2 was calculated from daily discharge data using the Flow Analysis Summary Statistics Tool ('fasstr') package in R (https://github.com/bcgov/fasstr). For some locations, only monthly discharge data are available. Thus, conversion of Q 2 from monthly to daily was applied for each climate region (Beck et al., 2018;Prasojo et al., 2022). The climate region for each delta is defined based on a Köppen-Geiger climate classification map (Beck et al., 2018). The predictive total system Q-W med relationships relate the median channel width measured for each delta as statistically representative values of right-skewed channel width distributions to the International Association of Sedimentologists., Sedimentology, 70, 362-380 total system discharge feeding into to the whole delta plain (Figs 2 and 3). The 66 median width values were obtained from 3823 individual measurements (mean number of width measurements per delta = 58; range from 15 to 177) (Fig. 2). Note that the data do not allow prediction of the discharge/palaeodischarge value of a single distributary channel, but enable calculation of the total riverine discharge that contributes sediment to build the delta plain. Ordinary least square (OLS) regressions were then used to calculate power-law scaling relationships between both median channel widths and catchment areas with bankfull discharge (Leopold & Maddock, 1953). OLS regression was used, which assumes error only in the dependent variable, as the aim is to produce predictive equations. The 95% confidence interval around the overall relationship for the 66 deltas is narrow, reflecting the statistical strength of the median channel width-bankfull discharge relationship across over three orders of magnitude of discharge. Using the regression equation to predict the discharge for an individual delta based on the estimate of the median channel width obtained from N width measurements yields a greater uncertainty (wider confidence interval) on account of the scatter in widths on individual deltas (blue shaded region in Fig. 2). The uncertainty in the median channel width estimate reduces as the number of width measurements increases since the uncertainty in the median decreases as a function of N −1/2 . OLS regressions were determined for each climate region to generate Q-A and Q-W morphometric scaling relationships.
The applicability of the power-law relationships determined from modern deltas was tested by applying the relationships to the channel widths and catchment areas derived from published outcrop data from Cretaceous formations in continental North America (Brownlie, 1983;Sageman & Arthur, 1994;Bhattacharya & MacEachern, 2009;Bhattacharya et al., 2016). Palaeodischarges were estimated in these studies using the Fulcrum method applied to outcrop and subsurface data of trunk rivers in incised valleys that fed downdip deltas. Due to the limited availability of measured trunk channel widths from the Ferron and Dunvegan formations, calculation of median channel widths is challenging. Hence, it may differ slightly from the methods that may be used to derive discharge values from the modern system. However, the power-law relationships the authors proposed from the study of modern deltas should be applicable to ancient deltas. This is because a similar approach was adopted by using total system discharge contributing to the whole delta plain that is analogous to discharge of a trunk channel feeding into a delta plain, as previously used in the literature.

Data distribution
Areas of the 66 catchments are log-normally distributed (Fig. 3A), similar to the global fluvial system dataset (pink and blue lines on Fig. 3A) (Milliman & Farnsworth, 2011). The fluvial catchment areas (Milliman & Farnsworth, 2011) are not significantly different from the delta catchment areas used in this study (t = 1.9; P < 0.06).
The median width of delta channels is almost one order of magnitude larger than the median in the global river channel width database (Allen & Pavelsky, 2018; Figs 3B and S1), although the range of widths are similar in both datasets. The channel widths in the present delta dataset are statistically significantly larger than in the fluvial data (t = −76.1; P < 2.2 × 10 −16 ). This difference suggests that scaling relationships from fluvial systems may not be able to be readily used for delta channels. Illustration of the use of median distributary channel widths to obtain the predictive relationships between median distributary channel widths and bankfull discharge from 66 river deltas measured in this study. The green arrow on the y-axis shows the uncertainty on the discharge estimation, while the green arrow on the x-axis shows the uncertainty on the width measurement.
Water discharge and catchment area scaling relationship (Q-A relationship) Globally, Fig. 4 shows a statistically significant (P = 3.3 × 10 −8 ; R 2 = 0.39; N = 66) power law relationship between catchment area and bankfull discharge, Q 2 = 50.1A 0.42 with 22 of the 66 deltas lying within the 95% confidence interval. Some of the more distant outliers are interpreted to be (D) boxplots of channel widths measured in this study and by Allen & Pavelsky (2018). In (A) and (B), N is the sample number, S sk is skewness, and t and P are the t-statistic and the associated probability from t-test comparison between the delta and fluvial datasets. The skewness values on (A) and (B) were calculated from the raw data, hence do not look skewed on log scales. present due to extensive river engineering (for example, embankments along riverbanks in Colorado, Nile and Ebro deltas) or due to wave and tide effects (for example, Orinoco, Mackenzie, Godavari, Ob and Irrawaddy deltas). In comparison to the global river Q-A relationship (Q = 0.075A 0.8 ), the scaling relationship for global deltas has a non-significantly lower regression slope (P = 0.1) using the significance of the difference, or slope test (Syvitski & Milliman, 2007).
The relatively low R 2 value for the global dataset can be explained in part by differences between climate regions. Separating the data into different climate regions produces significant relationships between A and Q 2 , except in arid and cold regions where the relationships are not significant (R 2 = 0.24 and 0.25; P = 0.13 and 0.069; N = 11 and 14, respectively).

Water discharge and median channel width scaling relationship (Q-W relationship)
In total, 66 paired measurements of discharge and median channel width were used to build the Q-W relationship. Overall, there is a statistically significant relationship with Q 2 = 0.34W med 1.48 (R 2 = 0.77; P = 2.2 × 10 −16 ; N = 66) (Fig. 5). The Q-W relationship produces a better fit globally than the Q-A relationship above (Fig. 4A). In comparison to the global river Q-W relationship (W = 17Q 0.45 ) (Moody & Troutman, 2002), the Q-W delta channel relationship has a statistically significant lower regression slope (P = 2.6 × 10 −5 ). As an example, predicting the discharge from a delta with median channel width of 300 m will result in Q 2 = 1576 m 3 /s, while the equivalent for a fluvial setting would be Q 2 = 589 m 3 /s. Deltas have multiple channels, hence using W med = 300 m will have maximum width of larger than 300 m near the apex (i.e. trunk channel), hence producing larger estimated bankfull discharge than the fluvial settings. These results suggest that predicting discharges from widths will produce different results if the channels are deltaic or fluvial.
When classified by climate region, Q-W relationships consistently show significant relationships (P < 0.05) with the strongest relationship for cold climates (N = 14) ( Fig. 5A and C). Polar, temperate and tropical regions also show strong relationships with R 2 values equal to 0.91, 0.88 and 0.63, respectively (Fig. 5D to F). Similar to the Q-A relationship, the Q-W relationship from arid regions (N = 11) shows the lowest R 2 although it is statistically significant (P = 1.2 × 10 −2 ) (Fig. 5B).
In summary, compared with the Q-A relationships on Fig. 4A to F, the Q-W relationships  (Syvitski & Milliman, 2007). The significance of the difference (slope) test between the gradient from delta Q-A OLS regression versus the global river from Syvitski & Milliman (2007) produces P = 0.1.
proposed in this study consistently show more statistically significant (P < 0.05) relationships that also have higher R 2 values (Fig. 5A to F). The Q-A relationship from the temperate region is the strongest (Fig. 4E) and the strongest Q-W relationship is for the cold climate region (Fig. 5C). The weakest relationships consistently come from the arid settings from both Q-A and Q-W (Figs 4B and 5B).

Application to the rock record
The scaling relationships obtained above from global modern river deltas are applied here to estimate palaeodischarges from several deltaic deposits. Data were compiled from palaeodischarge studies from well-exposed Cretaceous outcrops and subsurface datasets deposited in temperate-tropical climates. The data compiled from the literature used the Fulcrum approach to estimate palaeodischarge values (Table 1).
The Ferron Sandstone, exposed near Ivie Creek, south-west Utah, USA, is composed of Turonian (93.9-89.8 Ma) deltaic deposits from the western margin of the Western Interior Seaway (Bhattacharya & MacEachern, 2009;Braathen et al., 2018) (Fig. 6). The delta prograded north-east with an estimated drainage area of around 50 000 km 2 (Bhattacharya & Tye, 2004). Previous palaeodischarge studies on the Ferron Sandstone were based on trunk river characterization and estimation of palaeoflow velocity from its grain size, bedform size and estimated flow depth. The interpretation of a tropical palaeoclimate was obtained through facies analysis and catchment area is estimated from palaeogeographical reconstructions.
The Cenomanian (100.5-93.9 Ma) Dunvegan Formation was deposited in a temperate climate, and contains deposits from a large delta complex that are predominantly massive and cross-bedded non-marine and marine sandstones (Plint, 2002). The delta complex prograded 400 km north-west to south-east into the actively subsiding foreland basin of Alberta. It is estimated that the delta had a catchment area of around 100 000 km 2 (Bhattacharya & Walker, 1991;Sageman & Arthur, 1994;Plint, 2000;Bhattacharya & MacEachern, 2009;Hay & Plint, 2020).
Measured trunk channel widths and estimated catchment areas were obtained from the literature that compiled subsurface data with outcrop observations (Sageman & Arthur, 1994;Plint & Wadsworth, 2003;Bhattacharya & MacEachern, 2009;  Down-delta hydraulic geometry 371 Bhattacharya et al., 2016) and were used to calculate palaeodischarge using the equations calculated above from the modern systems. Four equations from this analysis of modern delta systems are used: (i) the global discharge-area relationship Q 2 = 50.1A 0.42 (Fig. 4A); (ii) the climate- ; and (iv) the climate-classified Q-W relationships, Q 2 = W med 1.4 for the tropical region and Q 2 = 0.07W med 1.66 for the temperate region ( Fig. 5E  and F). Palaeodischarges were calculated using these equations and channel widths measured from the rock record obtained from previously published work. The palaeodischarge values estimated using these equations were compared with previous palaeodischarge estimates (Fig. 7) (Sageman & Arthur, 1994;Bhattacharya & MacEachern, 2009;Bhattacharya et al., 2016).
These new estimates of bankfull discharges lie within one order of magnitude of the palaeodischarge values reported from the Fulcrum approach (Fig. 7) (Sageman & Arthur, 1994;Bhattacharya & MacEachern, 2009;Bhattacharya et al., 2016). Note that the climate-classified Q-W relationship provides a better fit to the previous estimates than the global Q-W relationship. Conversely, the global Q-A relationship estimates correspond better to previous estimates than do estimates from scaling relationships for individual climate zones (Fig. 4E and F). Overall, the statistical models proposed in this study perform similarly to the established Fulcrum method by producing values within the same order of magnitude as the palaeodischarge values derived from the literature.

Comparison to other palaeodischarge estimations
Analysis of river discharges, catchment areas and median channel widths from 66 river deltas has generated new global equations Q 2 = 50.1A 0.42 and Q 2 = 0.34W med 1.48 . These relationships have also been classified by five climate regions (Table 2). Fig. 7. Comparison of bankfull discharges estimated from previous studies with the estimated bankfull discharges calculated using the Q-W and Q-A relationships, both global and climate-specific, proposed in this study (Sageman & Arthur, 1994;Bhattacharya & MacEachern, 2009;Bhattacharya et al., 2016).
Applying these comparatively simple equations to the rock record produced palaeodischarge estimates within the same order of magnitude as the palaeodischarge values derived from existing, more complex approaches.
The new relationships proposed in this study allow quantification of palaeodischarge from the rock record based on measurements of channel width, estimates of palaeoclimate and morphometric scaling relationships derived from modern systems. The present approach uses fewer input parameters to estimate palaeodischarge than existing methods, the 'BQART' model or the Fulcrum model. Channel width is often measured from the rock record where channels are preserved and cross-channel exposures are available. The proposed morphometric scaling relationships simplify palaeohydrological calculations and enable more robust assessment of the uncertainties in the input parameter (channel width) to be accounted for when calculating palaeodischarges.
In comparison to the Fulcrum and 'BQART' models that implicitly include climate parameters, this work provides separate predictive equations for various climate regions. These proposed models show statistically significant correlations, especially between channel width and bankfull discharge across different climate regions, that have not previously been explicitly accounted for (Table 2). These climate-classified models will benefit source-to-sink studies by providing calculations tailored to individual palaeoclimates. Nyberg et al. (2021) provide a comprehensive overview of the uncertainties, sensitivities and practicalities of the 'BQART' model in estimating sediment load on geological timescales. They discussed in detail every parameter needed to estimate the palaeodischarge and palaeo-sediment load. For estimating the palaeodischarge, the 'BQART' model uses a global Q-A power law scaling relationship similar to this study but without explicitly allowing for climate. Eide et al. (2018) added runoff (Ro) parameters to take into account the impact of climate by applying a different multiplier value to discharges calculated for each climate region (for example, Ro = 0.0005 for arid and Ro = 0.0161 for humid regions). However, adding Ro constants shift the models, but does not change the models' gradients. In contrast, different equations were produced for each climate region, allowing the models' gradients to change, reflecting the role of soils or vegetation in controlling runoff. Also, the climate-classified models proposed in this study make palaeodischarge estimation more straightforward if the palaeoclimate can be deduced from the rock record.
Although the equations are statistically robust, defining palaeoclimate from the rock record is not straightforward due to the often sparse exposure of preserved channels, complexities in stratigraphic correlation and the need for palaeoclimate evidence. Reconstructing the relationship between evidence requires significant effort and may not always yield conclusive results (Shuman, 2014). Hence, it is reasonable to assess whether the present climate-specific equations significantly improve palaeodischarge estimates. Analysis of variance (ANOVA) tests were used to compare the global Q-W and Q-A regression equations to the climate-classified Q-W, Q-A relationships. Comparing the global and the climate-classified Q-W regression lines produced P = 0.62. While the comparison of the global Q-A and climateclassified Q-A regression lines produced P = 0.07. Both of the tests showed that both global and climate-classified Q-W and Q-A relationships are not significantly different, hence could be used interchangeably. The tests imply that when the palaeoclimate is challenging to be deduced from the rock record, the global Q-W or Q-A scaling relationship could be used instead.

Limitation of the proposed scaling relationships
For the Q-A and Q-W relationships, the standard errors of residuals are 1.23 and 0.76 in log units, respectively. Despite overall significance of the regressions, additional factors may affect both relationships, such as anthropogenic effects on channel width and/or river flows that may disrupt the dynamic equilibrium assumption that underpins the proposed scaling relationships (Aslan et al., 2005;Li et al., 2017;Ninfo et al., 2018), vegetation type and density (Huang & Nanson, 1997), sediment load (Hey & Thorne, 1986), grain size (Eaton, 2013), anabranches of multi-thread channel systems (Tabata & Hickin, 2003), material forming the channel boundary (Ellis & Church, 2005) and flood variability for each climate region (Rodier & Roche, 1978). Although the accuracy of predictive models can be improved by adding more variables (Mosley, 1981) this addition leads to models becoming increasingly less applicable to the rock record. For example, using the calculations herein, palaeodischarge can be determined from any dataset in which a catchment area or channel width can be determined (for example, outcrop or seismic). However, if other variables such as grain size or palaeoslope are needed these additional data may not always be readily available. Thus, keeping the variables as simple as possible (for example, catchment areas and channel widths) is beneficial in creating models that are applicable to the rock-record. Also, adding more variables does not necessarily result in an increase in model accuracy. Mosley (1981) showed that 88% variability of channel crosssectional area (for example, width, depth) is controlled by the bankfull discharge, bed sediment size and bank sediment character, with the rest of the variability being explained by morphological variables (for example, braiding and sinuosity index). For reconstruction purposes, there is merit in simplicity and careful examination of the contributing factors of channel cross-section and bankfull discharge should be undertaken before adding in more variables into the morphometric scaling relationships proposed in this study. The prediction intervals for palaeodischarge (Table 1; Fig. 2) are wide because of scatter in the observations and the small number of width measurements available for the prediction. These wide prediction intervals need to be acknowledged when using the proposed scaling relationships. Consequently, when applying the scaling relationships to the rock record they should be further constrained as far as possible using all contextual information gathered from the rock record (for example, grain size, bedforms interpreted from sedimentary structures, stratigraphic position) to justify the palaeodischarge estimation produced by this approach. The present source dataset of modern measurements is spatially distributed across the delta in one time horizon, from which a median width to use for prediction was determined. However, deltas are depositional systems and due to transgression/regression measurements made in outcrop or from subsurface imaging may produce biased samples across the delta, hence yielding a biased estimate of median channel width, or may aggregate measurements across time horizons with different external controls, such as changing Q 2 due to climatic fluctuations. Hence, as noted above in the context of climate interpretation, applying the new statistical models to the rock record requires interpreting the stratigraphic context of the measured distributary channels. As more data become available, larger datasets, modern and ancient, will be able to be used to constrain what are 'reasonable' palaeodischarge estimates. This constraint will be quantitative as more datasets such as those in Table 2, are obtained.
This study's approach uses width measurements from satellite imagery because width is the most readily obtained measure of channel scale on such images. In outcrop or subsurface datasets it is commonly easier to measure distributary channel depths (d) than widths. Channel widths and depths are very highly correlated empirically and theoretically (Ferguson, 1986), and depth and width measurements from distributary channels reported in the literature are summarized in Table S1. Although the depth and width exponents in Eqs 1 and 2 are consistent, variations in the multipliers mean that the W:d ratio cannot be taken as a global constant due to the influence of additional factors on channel geometry (for example, vegetation, bank sediment cohesion). Some studies have found that in deltas W:d varies with discharge (Wang & Li, 2011) or with the measurement location .
To accommodate the complexity of the relationship between channel width and depth in river deltas, it is assumed here that the flow was steady during the W and d measurements in Table S1, and that width was in equilibrium with depth and slope. The range of measured W:d ratios is from 10 to 200 with typical values of W:d = 30:1 [for example, Mississippi delta;  and Fly delta (Latrubesse, 2008)], to 100:1 [for example, Yellow River delta (Wang et al., 2008;Wang & Li, 2011), Amazon and Brahmaputra deltas (Latrubesse, 2008)] and the extreme values of W: d = 200:1 from Wax Lake and Lena deltas (Olariu & Bhattacharya, 2006) (Table S1). By assuming that delta channel W:d relationships globally lie within the range suggested by the available measurements, the Q-W relationship was rescaled to yield a novel discharge-depth (Q-d) scaling for W:d = 30, 100 and 200 (Table S2). Then the Q-d scaling was classified based on climate regions. All of the Q-d relationships are statistically significant, with different W:d ratios affecting scaling constants only.

Climate impacts on the proposed scaling relationships
Climate-classified Q-W relationships may produce more reliable palaeodischarge results than either the Q-A relationships or the global Q-W relationship due to the direct impacts of climatic factors on channel geometry. Most of the climateclassified Q-W relationships have higher R 2 values (0.52-0.94) than the global Q-W relationship (R 2 = 0.77). The Q-A relationships have R 2 values of 0.39 for the global data, and 0.24-0.85 for the climate-classified relationships. This does not necessarily mean that Q-A relationships should not be used, but depending on the data availability from the rock record, both Q-A and Q-W relationships remain useful for inferring palaeodischarge.
Climate-classified Q-A relationships should give more reliable predicted palaeodischarges than a single global Q-A relationship due to discharge being directly controlled by rainfall and runoff in each climate region (McCabe & Wolock, 2016;Eide et al., 2018). However, for palaeodischarge studies, catchment areas calculated from palaeogeographical reconstructions may contain significant uncertainties due to the assumptions and interpretations involved in building palaeogeographical maps. Hence, the ability to estimate palaeodischarge through regional hydraulic geometry scaling relationships (Davidson & North, 2009) supported by provenance analysis (Blum et al., 2017), remains constrained by scatter in the modern data and the need to supplement the calculations with further estimated variables. Errors of at least one order of magnitude are not uncommon (Bhattacharya et al., 2016), but may provide valuable information that cannot be obtained by other means, or that supplements independent reconstructions of palaeoenvironments.
Particular caution is required when estimating palaeodischarge in arid and cold climates. Arid climates have annual rainfall between 150 to 200 mm (Thornthwaite, 1948) and a highly episodic runoff regime with flood flows lasting for only a few hours or days in a year (Rodier & Roche, 1978). This regime makes the definition of bankfull discharge challenging in this climate (Shamir et al., 2012). As an example, it is common to have rapid intermittent high flood with low and steady flow periods throughout the year in an arid region (for example, due to snowmelt in Colorado river catchment and intermittently anabranching river during low flow) (Segura & Pitlick, 2010). Catchment area and bankfull duration are poorly correlated in arid regions (Dodov & Foufoula-Georgiou, 2005), and interannual runoff irregularity and downstream loss of water are very significant in these regions (Rodier & Roche, 1978).
In cold climate regions, flow may be noncontinuous or substantially reduced in winter so reducing how representative Q 2 is as the bankfull discharge, also resulting in a weak correlation between catchment area and Q 2 (Beltaos & Prowse, 2009;Stonevičius et al., 2014). Flood hydrology in this region depends on interactions between snow and ice cover, precipitation and air temperature, that may induce shifts in runoff over decadal timescales (Stewart et al., 2005;Shiklomanov et al., 2007). Consequently, bankfull discharge estimation from both modern and ancient systems in these two climate regions should consider hydrological conditions in the relevant climate zone.

Further developments
Although the relationships calculated herein produce realistic discharge estimates in Cretaceous deltas constrained by outcrop and subsurface data, there is a need to test these relationships across different aged systems across varying climate belts to understand the extent to which they can be applied. Also, despite scaling relationships being available from modern estuaries (Diefenderfer et al., 2008;Gisen & Savenije, 2015) and tideinfluenced river deltas (Sassi et al., 2012), development of similar rock-record focused scaling relationships for other systems (for example, tidal creeks or other delta types) remains an area for further study.
Finally, the method proposed here adopts metrics that are more easily extracted from the rock record and which is based on specific climate zones has potentially important implications with regard to assessment of hydrocarbon, hydrogen, geothermal and carbon capture and storage (CCS) sizes (Bhattacharya & Tye, 2004). In addition, it will help in deducing climate and tectonic forcing on systems and palaeohydraulics across various types of depositional systems in source-to-sink studies (Montgomery & Gran, 2001;Merritt & Wohl, 2003;Brardinoni & Hassan, 2006;Wohl & David, 2008;Davidson & Hartley, 2010;Eaton, 2013).

CONCLUSION
This study has obtained Q-A and Q-W med scaling relationships for 66 modern river deltas across different climate regions by extracting catchment areas for each delta, making 3823 distributary channel width measurements and calculating their associated total system discharges. These relationships are intended to provide quantitative information on source catchment properties from data typically available in the rock record. Applying the simple scaling relationships derived here from modern systems to the rock record, coupled with palaeoclimate information, produced palaeodischarge estimates within the same order of magnitude as palaeodischarge values derived from existing, more complex, approaches that require a larger number of parameters. These new relationships promise enhanced deduction of climate and palaeodischarges across various types of depositional systems in source-to-sink studies, assessment of hydrocarbon, hydrogen, geothermal and carbon capture and storage (CCS) sizes, and more accurate palaeogeography interpretations. The relationships have been tested against data from some Cretaceous deltas, applying these scaling relationships to other palaeoclimate regions, systems of different ages and to different types of deltaic environment, remain areas of further study.

Supporting Information
Additional information may be found in the online version of this article:

Table S1
In-situ measurement of width and depth of several river deltas collected from the literature.

Table S2
Global and climate-classified discharge: depth (Q-d) scaling relationships proposed based on three width:depth (W:d) ratios found from a number of modern river deltas.