University of Birmingham The role of discharge variability in determining alluvial stratigraphy

We illustrate the potential for using physics-based modeling to link alluvial stratigraphy to large river morphology and dynamics. Model simulations, validated using ground penetrating radar data from the Río Paraná, Argentina, demonstrate a strong relationship between bar-scale set thickness and channel depth, which applies across a wide range of river patterns and bar types. We show that hydrologic regime, indexed by discharge variability and flood duration, exerts a first-order influence on morphodynamics and hence bar set thickness, and that planform morphology alone may be a mislead-ing variable for interpreting deposits. Indeed, our results illustrate that rivers evolving under contrasting hydrologic regimes may have very similar morphology, yet be characterized by marked differences in stratigraphy. This realization represents an important limitation on the application of established theory that links river topography to alluvial deposits, and highlights the need to obtain field evidence of discharge variability when developing paleoenvironmental recon-structions. Model simulations demonstrate the potential for deriving such evidence using metrics of paleocurrent variance.


INTRODUCTION
Alluvial deposits are a key archive for reconstructing river morphology, hydrology, and paleoenvironments (Blum and Törnqvist, 2000;Miall, 2006). However, interpretation of deposits is commonly difficult due to the lack of unambiguous criteria linking fluvial processes to sedimentary product (Bridge, 2003;Ethridge, 2011), and because the stratigraphic record is incomplete (Strauss and Sadler, 1989). Quantitative theory has been used to relate bedform geometry and dynamics to bedset thickness (Paola and Borgman, 1991). However, studies have necessarily focused on small spatial scales, such as laboratory settings (Straub et al., 2012;van de Lageweg et al., 2013), or deposits associated with individual bedform trains (Bridge and Best, 1997;Leclair, 2011). Moreover, existing theory neglects the role of hydrologic variability (e.g., flood magnitude and duration), despite its importance as a control on river evolution and deposit reworking (Tamminga et al., 2015) and bar and bedform geometry (Wilbers and Ten Brinke, 2003), all of which determine the resultant stratigraphy. Recent work highlights a need to understand better the link between morphodynamics and sedimentology, particularly at bar and channel-belt scales, and for a range of river patterns and hydrologic regimes (Fielding et al., 2009;Ethridge, 2011;Plink-Björklund, 2015). Achieving this aim has proven virtually impossible to date due to a lack of suitable field data sets. However, recent advances in numerical modeling mean that it is now possible to simulate river morphodynamics over large temporal and spatial scales (Nicholas, 2013;Schuurman et al., 2013). Here, we aim to (1) evaluate the potential for models to generate spatially rich data sets quantifying alluvial architecture; (2) elucidate the roles of hydrologic regime and river pattern as controls on the resultant stratigraphy; and (3) identify some key limitations on the application of existing theory linking alluvial deposits to their formative flows.

APPROACH
The deposits of large sand-bed rivers were simulated using a physicsbased numerical model of hydraulics (for subcritical and supercritical flows), sediment transport, bank erosion, and floodplain formation. This model, described and evaluated elsewhere , is suitable for representing meandering, braided, and anabranching channels (Nicholas, 2013). Twenty-six (26) simulations were conducted herein using a range of bed slopes, sediment loads, and bank erodibilities to generate rivers (50 km in length) with contrasting channel patterns. All 26 simulations used the same hydrologic regime (flood hydrographs where discharge varied from a low of 10,000 m 3 s -1 to a peak of up to 30,000 m 3 s -1 ). In all simulations, the river evolved from a straight initial channel of constant width. Simulation duration (typically 175 floods, nominally equivalent to a scaled time period of 350 yr) was sufficient to rework deposits multiple times. Here, we focus on six simulations that generated low-sinuosity anabranching channels similar in form to the Río Paraná, Argentina ( Fig. 1), for which we have characterized the deposits of kilometer-scale bars using ground penetrating radar (GPR) and sediment cores up to 5 m in length (Reesink et al., 2014). The Río Paraná has a mean annual discharge of 17,000 m 3 s -1 at Corrientes, Argentina, where the geophysical surveys were based. To investigate the influence of hydrologic regime on the stratigraphy, these six simulations were also run with floods in which hydrograph duration was increased by factors of two and four, and additional simulations that used a constant discharge of 22,500 m 3 s -1 (the average peak discharge for simulated floods), yielding 44 simulations in total (see Table DR1 in the GSA Data Repository 1 ).
Modeled deposits were reconstructed from channel topography and flow conditions at 700 points in time over the course of the simulations. Bedsets, defined as depositional elements bounded by erosional surfaces (Straub et al., 2012), were identified from vertical profiles in each model grid cell (80 × 40 m in size). Modeled sets are associated with macroscale morphologic features (e.g., unit bars) represented by the model topography, rather than smaller bedforms (e.g., dunes) that are finer than the model grid resolution. Deposits were subdivided into three classes (termed "slackwater," "ripples," and "dunes") based on modeled flow conditions. Slackwater deposits were classified as those that form below a velocity threshold of 0.1 m s -1 , and the criterion of van Rijn (1984, see his figure 1) was used to define the ripple-dune transition. To account for the existence of non-equilibrium dunes (e.g., on the falling limb of a flood), deposits were only classified as ripples where the threshold for dune formation was not exceeded at any point during the hydrograph. Simulations are characterized by zero net aggradation, and hence total deposit thickness scales with maximum thalweg depth (typically 25-30 m). Analysis of deposits is restricted to sediment below the vegetation that is established on bar tops that are inundated infrequently.

RESULTS
Simulations that yield low-sinuosity anabranching rivers (e.g., Fig.  1A), similar in form to the Río Paraná (Fig. 1B), are characterized by kilometer-scale sand bars that grow by vertical stacking of unit bars, and lateral accretion of bar wings that wrap around the bar head. Modeled bar deposits (Fig. 1C) are composed of stacks of four to eight bar sets (similar to the three to seven bar sets reported by Bridge and Lunt [2006]). Cross-bar channel fills and slackwater sediments deposited in the lee of the bar are common in simulations. Truncation of deposits by unit bar migration is common, as is reworking to depths of 5-10 m below the bar surface (Fig. 1D). Comparison of model results with GPR data from the Río Paraná (see Table DR2) indicates that the model reproduces both the vertical dimensions of bar sets and the tendency for sets to thin toward the bar surface (Fig. 1E). Modeled deposits comprise 1%-3% slackwater sediments (predominantly composed of silt) and 5%-30% ripples, compared to 30% ripples and 3% silt/clay (deposited in slackwater areas), on average, for bars from the Río Paraná near Corrientes (Reesink et al., 2014).
Previous studies have applied existing theory (e.g., Paola and Borgman, 1991) to relate set thickness to formative flow depth for largescale strata generated by migrating bars (Bridge and Lunt, 2006;van de Lageweg et al., 2013). Such analysis commonly involves the assumption that the spatial distribution of bed topography at an instant in time is representative of the temporal distribution of topography at a point in space (i.e., that morphology is a reliable measure of morphodynamics). We demonstrate below that this assumption may be unjustified. Despite this, we observe a strong positive relationship between mean channel water depth, calculated as the average depth at all channel locations and model time steps, and mean bar set thickness for all 26 simulations conducted using the same variable hydrological regime ( Fig. 2A). These simula-  tions are associated with a wide range of channel patterns and widths (total channel width varies from ~1.5 km to 7 km). We find no statistically significant difference in the ratio of mean bar set thickness to mean flow depth between channels with low and high width:depth ratios (n = 13 for both groups). Moreover, the transition from wider and shallower (anabranching) to narrower and deeper (meandering) channels is associated with a transition from unit bar-dominated to scroll bar-dominated deposits. These results imply a near constant ratio of mean bar set thickness to mean flow depth irrespective of bar type and channel pattern.
The six simulations with constant discharge plot well above the regression line in Figure 2A, with bar set thickness for these simulations increasing by a factor of 1.6, on average, compared to equivalent simulations where discharge varies. This cannot be explained fully by differences in morphology (e.g., channel width, depth, or pattern). Moreover, simulations run with constant and variable discharge experience similar average rates of deposition and bed reworking over decadal time scales. Despite this, bar set thickness exhibits a clear relationship with channel morphodynamics, as defined by measuring the thickness (Dz) of packages of continuous erosion or deposition in all individual model grid cells throughout the simulations, in order to derive a probability density function of morphodynamic event magnitude (Fig. 2B). Simulations with constant discharge experience an increase in both small-scale (|Dz| < 0.025 m) and large-scale (|Dz| > 2.25 m) erosion and deposition events but a reduced frequency of intermediate-scale events, and an overall increase in the variance of vertical change increments. We attribute this to two factors. First, bars aggrade until reaching the water surface, and hence when discharge is constant and water level changes are small, many bar surfaces experience lower rates of vertical change. Second, cut-and-fill cycles driven by flood hydrographs are absent under constant discharge because temporal changes in flow velocity (at any given location) are limited. This allows the duration and magnitude of continuous deposition events to increase, thus promoting thicker sets. Similarly, where flood hydrograph duration increases, periods of continuous deposition are also longer. This promotes a positive relationship between flood duration and relative bar set thickness (Fig. 2C). Significantly, the standard deviation of Dz values is an excellent predictor of mean bar set thickness across all 44 simulations (Fig. 2D). Thus, while mean bar set thickness is a function of mean channel depth, morphodynamics rather than morphology is the dominant control on stratigraphy.
Further insight into these relationships can be derived by analysis of the variability in paleocurrent directions associated with the deposits (defined by the modeled velocity vectors at the time of sediment deposition) and by the ratio of the downstream and cross-stream dimensions of facies units (defined as deposits characterized by similar proportions of dunes, ripples, or thick sets [see metrics used in Figure 3 and Table 1, and in Table DR3]). Simulations that use a variable discharge regime (flood duration, T = 2 yr) are characterized by distinct values of these metrics for both low-and high-sinuosity channels. Moreover, low-sinuosity anabranching channels generated by variable and constant discharge regimes exhibit marked differences in deposit characteristics, despite having similar morphology. Channels formed by constant discharge ex-  Note: Columns 2 and 3 show results for two simulations with contrasting morphology (see Fig. 3). Columns 4 and 5 show the mean of results for six simulations of anabranching channels that use variable discharge (hydrograph duration, T = 2 yr) and constant discharge, respectively. λ is the mean set thickness. Lxy is the ratio of the average downstream and cross-stream lengths of contiguous model grid cells classified by deposit type as: dunes (cells where >90% of sediment is classed as dunes); ripples (cells where >10% of sediment is classed as ripples); and large sets (cells where >50% of sediment comprises sets thicker than twice the mean set thickness). σ V90 is the 90 th percentile of the probability density function of the standard deviation of paleocurrent direction. ψ is the mean thickness of contiguous vertical packages of each deposit type. hibit lower variability in paleocurrent direction and facies units that are preferentially elongated in the downstream direction. Vertical packages of each deposit type (notably dunes) exhibit marked differences in thickness between contrasting channel planforms (see Table 1). Moreover, where discharge is constant, sediment packages are thicker on average compared to those generated under a variable discharge regime. This is consistent with the inverse relationship between unit bar set thickness and discharge variability suggested previously by Sambrook Smith et al. (2009) and indicates a tendency for bar overtopping to be inhibited where discharge is constant. This limits the occurrence of lateral bar-top flows and channels, and promotes flow streamlining that encourages the elongation of deposits.

SUMMARY
This study illustrates the potential for using physics-based modeling to link river morphodynamics to stratigraphy. Our results demonstrate that bar set thickness is a good predictor of channel depth, irrespective of river channel pattern and associated differences in bar type. However, depth estimates derived from bar set thickness data may be highly sensitive to the hydrologic regime. This suggests that paleoflow reconstructions should attempt to assess the nature of discharge fluctuations, for instance, as expressed by reactivation surfaces that are not associated with bedform superimposition.
Our simulations examine large rivers, such as the Rio Paraná, characterized by low discharge variability (ratio of annual range in discharge to mean discharge, Q var < 2) and gentle slopes, that are dominated by subcritical flows. Herein, we do not consider rivers characterized by high discharge variability or flash floods (e.g., Q var > 100; Fielding et al., 2009), where deposits formed under supercritical flow may be abundant and accretionary sets associated with bar migration may be poorly developed (Plink-Björklund, 2015). Consequently, our conclusions regarding the significance of hydrologic regime are almost certainly conservative.
Our results indicate that data quantifying paleocurrent variance and the downstream and cross-stream dimensions of facies units may be valuable for constraining hydrologic variability. However, such characteristics are also a function of channel pattern, in particular sinuosity. These results also indicate that physical and numerical models that impose a constant flow discharge may not be simulating correctly the alluvial architecture of natural channels that experience discharge variability.
When relating bed topography to set thickness, some studies (e.g., van de Lageweg et al., 2013) have assumed that the spatial distribution of bed heights (in bathymetric data) is representative of the temporal distribution at a point in space. Our results demonstrate that this need not always be true. Modeled rivers with similar morphology can be characterized by significant differences in temporal dynamics and hence stratigraphy. Moreover, while a positive relationship between topographic variability and set thickness is central to accepted theory (e.g., Paola and Borgman, 1991), we find that increased hydrologic variability suppresses bar set thickness due to its influence on morphodynamics. Hydrologic regime thus plays a key role in controlling stratigraphy that has yet to be incorporated within predictive theory. This implies that use of stratigraphic evidence to link paleoenvironment to morphology can only succeed by giving consideration to the essential role of dynamics as a control on sediment accumulation and preservation.