Spatio‐temporal reconstruction of avulsion history at the terminus of a modern dryland river system

Fluvial depositional architecture in an unconfined environment is governed by sediment dispersal across the alluvial plain through river‐path switching by avulsion. Documented inter‐avulsion periodicity from modern rivers ranges from tens to over a thousand years. In this study, a quantitative spatio‐temporal reconstruction of avulsion history is presented of the non‐vegetated and pristine modern Río Colorado dryland river system in the semi‐arid Altiplano Basin (Bolivia), based on the integrated analysis of satellite imagery and absolute age dating using optically stimulated luminescence, complemented with sedimentological and geomorphological ground‐truth data. This approach enables us to reconstruct the chronological order of channel belts of the Río Colorado, to determine avulsion recurrence time and inter‐avulsion periodicity, to identify mechanisms for flow path changes, and to present a morphodynamic model for the spatio‐temporal evolution of fluvial deposits in a semi‐arid environment. In a maximum timespan of 12.71 ± 1.5 ka, successive avulsions of the Río Colorado created a sheet of interconnected fluvial deposits, consisting of diverging and juxtaposed alluvial ridges that formed by sediment aggradation in point bars, crevasse splays, levees, and on the channel floor. The ridges show lateral onlap and amalgamation as the result of repeated avulsion and compensational stacking, whereby the river avoided the positive alluvial‐ridge relief of its precursors. The resultant morphology is fan‐shaped, convex‐up with a surface area of approximately 500 km2 and a maximum observed thickness of 3 m. The results show inter‐avulsion periods of the river of up to 1.28 ± 0.34 ka. A paucity in fluvial activity around 2 ka BP, and at present, is interpreted as the result of low river discharge related to long‐term dry periodicity in the El Niño Southern Oscillation circulation system. Each river path started as a low sinuous, single‐thread channel in a narrow belt, and in time increased its width and sinuosity by point‐bar expansion and rotation.

. Avulsions are controlled by the interplay of autogenic processes such as river hydrodynamics, sediment composition and supply rates, floodplain morphodynamics, breaching of river levees, alluvial-ridge aggradation, differential compaction of mud-prone floodplain with respect to sand-prone alluvial-ridge sediment, and related changes in along-river and cross-valley gradient (Aslan et al., 2005;Edmonds et al., 2016;Hajek & Edmonds, 2014;Hajek & Wolinsky, 2012;Jones & Schumm, 1999;Mohrig et al., 2000;Nicholas et al., 2018;Slingerland & Smith, 1998;van Toorenenburg et al., 2018). He et al. (2020), Stouthamer and Berendsen (2007), and Timár et al. (2005) documented allogenic processes (sea-level change, tectonics) as a cause of river avulsion. Berendsen and Stouthamer (2000) and Cohen et al. (2005) documented differential subsidence across Roer Valley Graben faults in the North Sea Basin as an independent local control on the occurrence of avulsions in the Rhine-Meuse Delta. Slingerland and Smith (2004) presented a comprehensive review of the different types of avulsion based on Holocene and modern examples, one of which, termed 'progradational avulsion', is a response to alluvial-ridge aggradation and associated decrease in channel gradient and flow efficiency (Hajek & Wolinsky, 2012;Kleinhans et al., 2013;Slingerland & Smith, 2004;Törnqvist & Bridge, 2002). Alluvial ridges aggrade through deposition in levees and crevasse splays, and associated accumulation of channel-floor sediment, super-elevating the river thalweg relative to the surrounding floodplain surface until the river avulses to a hydrodynamically more favourable, steeper path (Ethridge et al., 1999;Hajek & Wolinsky, 2012;van Toorenenburg et al., 2016van Toorenenburg et al., , 2018. Mohrig et al. (2000) found that this process of avulsion occurs at a normalized super-elevation (i.e. the ratio between levee height above the floodplain divided by channel depth) of 1.0 AE 0.5. From a flume-experiment study of channel avulsion in a backwater-controlled deltaic setting, Ganti et al. (2016) concluded that in-channel sedimentation, rather than super-elevation, was the prime control on the occurrence of river avulsion, with a depositional depth of 0.3 AE 0.13 of the total channel depth as avulsion threshold.
Avulsions have been documented from modern rivers such as the Kosi River in India (Chakraborty et al., 2010;Sinha, 1996Sinha, , 2009Wells & Dorr, 1987), Brahmaputra River in Bangladesh , Rhine-Meuse Delta in the Netherlands (Stouthamer & Berendsen, 2001), Saskatchewan River in Canada (Slingerland & Smith, 2004;Smith et al., 1989), Yellow River in China (Ganti et al., 2014;Slingerland & Smith, 2004), and Tshwane River in South Africa (Larkin et al., 2017). Valenza et al. (2020) documented downstream changes in avulsion style in various remotely sensed modern rivers, from reoccupation of abandoned braided channels in the upstream part of the river trajectory, to flooding and sedimentation during avulsion-channel formation in the downstream meandering river part. Documented avulsion frequencies range from once every 0.028 ka (Kosi River; Slingerland & Smith, 2004) to four occurrences in 5 ka (Mississippi River; Aslan et al., 2005). Studies of preserved avulsion deposits in the rock record are few (e.g. Berendsen & Stouthamer, 2000;Coronel et al., 2020;Flood & Hampson, 2014;Jones & Hajek, 2007;Kraus & Wells, 1999;van Toorenenburg et al., 2016) and yield limited information on the mechanisms and frequency of avulsion. The cause of this shortcoming may be that: (1) avulsion frequencies in the rock record are expected to be beyond the resolution of existing geochronological methods, given the short avulsion recurrence time interpreted from modern settings; and (2) the lack of reliable correlation surfaces hampers the interpretation and quantification of river gradient, floodplain aggradation, superelevation, or avulsion geochronology. Palaeosols have been proposed as correlation surfaces in the rock record (Hajek & Wolinsky, 2012;Kraus & Wells, 1999). The analysis of the chronological order of river avulsions from the stratigraphical position of channel deposits depends on the assumption that palaeosols formed on a topographically flat and horizontal floodplain surface, and that post-depositional differential compaction has not altered the original geomorphology. These uncertainties may impair the use of palaeosols as a reference surface for quantification of avulsion parameters in the rock record.
In this study, we aim to: (a) quantitatively reconstruct the order, chronology, and recurrence time of avulsions of the Río Colorado over the last $12 ka; (b) discuss the mechanisms for the flow-path changes; and (c) present a morphodynamic model for the spatiotemporal evolution that captures the organization of aggrading fluvial deposits in a semi-arid environment. To this end, we integrate remote-sensing-based mapping of the modern Río Colorado dryland river system in the semi-arid Altiplano Basin (Bolivia) with absolute age dating of fluvial sediment samples using optically stimulated luminescence (OSL), and with sedimentological and geomorphological ground-truth observations and high-resolution differential-GPS (dGPS) surveys. The results of this study will broaden our understanding of fluvial morphodynamics and sediment transport processes in a pre-vegetation earth (cf. Ganti et al., 2019;Ielpi, 2019;Ielpi & Lapôtre, 2019Ielpi, Fralick, et al., 2018;, and in the context of climate-change-related aridification of alluvial basins.

| STUDY AREA
The study area is located in the southeast part of the endorheic Altiplano Basin (Bolivia), which is a high-altitude hinterland plateau (3650-4200 m above mean sea level) enclosed by the Eastern and Western Andean Cordilleras (Figure 1; Horton et al., 2001). The basin fill consists of Cretaceous to Holocene alluvial and lacustrine sediments (Elger et al., 2005). Alternating wetter and drier climate periods during the Pleistocene and Holocene resulted in short-duration, highamplitude base-level fluctuations in the basin (Baker et al., 2001;Chepstow-Lusty et al., 2005;Fornari et al., 2001;Placzek et al., 2006;Rigsby et al., 2005;Servant et al., 1995;Sylvestre et al., 1999). The fluvial deposits of the Río Colorado dryland river system accumulated after the last base-level highstand of 13.4-11.5 ka (Baker et al., 2001) when the Coipasa Palaeolake was up to 55 m above the present-day base level (Bills et al., 1994;Placzek et al., 2006). The total area of fluvial sedimentation is approximately 500 km 2 and attains a maximum observed thickness of 3 m (Donselaar et al., 2013). The river gradient is extremely low (averaging 8.3 Â 10 À5 ) with an alongchannel drop in thalweg elevation of 3 m over 36 km. At present, the basin is in a dry-climate period with an average annual precipitation of 0.20 m/a (Argollo & Mourguiart, 2000) and an evapotranspiration potential of 1.5 m/a (Grosjean, 1994;Risacher & Fritz, 2009). As a consequence, lakes on the southern part of the plateau are ephemeral or have dried up, and river systems such as the Río Colorado now terminate on the alluvial plain bordering the present-day salt lake ( Figure 2). The study area is extremely flat, in part due to deflation in the dust storm events of July 2009 and July 2010 (Gaiero et al., 2013), and is a non-vegetated, pristine area, without anthropogenic alteration.
The Río Colorado is an endogenic ephemeral dryland river (cf. Bull & Kirkby, 2002) in a semi-arid climate. Annual rainfall in the study area is less than 0.2 m/a (Argollo & Mourguiart, 2000), and the evapotranspiration potential is 1.5 m/a (Grosjean, 1994;Risacher & Fritz, 2009). The river drains part of the Eastern Cordillera and is characterized by prolonged periods (multiple months up to years) of drought alternating with short (less than 24 h), severe peak-discharge events (Donselaar et al., 2013;Garreaud et al., 2003;Li, 2014;Li et al., 2014). Peak discharge in the study area takes place in the austral summer months (December to March) with short periods of intense rainfall (Donselaar et al., 2013). Transmission losses are high in the dry-climate setting due to flood-outs onto the adjacent floodplain, evapotranspiration, and infiltration in the floodplain surface (Costa et al., 2012(Costa et al., , 2013Jarihani et al., 2015;Knighton & Nanson, 1994;Lange, 2005). In combination with the very low thalweg gradient in the terminal part of the Río Colorado, this results in a downstream decrease of the river width and depth (Donselaar et al., 2013). In peak-discharge periods, the flood water spills out of the channel confinement via multiple levee breaches and forms extensive, laterally amalgamated crevasse-splay lobes fringing the river (  The geomorphological expression of channel belts in the study area was captured with field measurements and Google Earth-Pro 2016 satellite imagery, with an improved resolution with respect to the 2004 images used by Donselaar et al. (2013). Laser rangefinders with 20 cm resolution were used to verify crevasse-splay radii and point-bar radii on the ground, and to measure river width and depth ( Figure 3). In addition, subtle surface topography-such as the height difference across alluvial ridges and to the adjacent floodplain-was recorded at centimetre-scale resolution in dGPS surveys. The sinuosity index (SI) of channels was established in Google Earth-Pro as the ratio of channel-thalweg length to channel-belt length (cf. Leopold & Wolman, 1957;Schumm, 1963). Satellite imagery and laser rangefinder measurements yielded the channel widths and point-bar radii.
The absolute age of each abandoned channel belt was deter- black tape to prevent exposure to daylight. OSL measurements and analyses were carried out at the Netherlands Centre for Luminescence dating. Equivalent doses were obtained by measuring OSL signals on small aliquots of sand-sized quartz grains (Murray & Wintle, 2003), and palaeodoses were then derived by applying the bootstrapped minimum age model (BSMAM) to the equivalent dose distribution (Cunningham & Wallinga, 2012) (see online Supplementary Information for details). Sample OSL age was defined as the ratio of palaeodose to dose rate. Note that all time values in this paper are denoted as best estimate AE standard error (68% confidence interval).
The inter-avulsion period (i.e. the absolute time between consecutive avulsions) was estimated as the difference between OSL abandonment ages of consecutive channel belts. These interval times are presented as a normal difference distribution with its lower tail truncated at zero, given that the relative order of the sampled channel belts is known (hence, the difference cannot be negative). Note that this interval time includes any period of coeval activity (e.g. during gradual avulsion).
The active migration time of each channel belt was estimated as the quotient of maximum point-bar radius over estimated migration rate, based on the empirical power-law correlation between average channel width and migration rate (Ielpi & Lapôtre, 2020): depression (e.g. draining the inundated floodplain) was assumed to be negligible.

| Fluvial sedimentology and geomorphology
The improved resolution of the Google Earth-Pro imagery allowed for the expansion of the number of different channel belts from 8 ( figure 13 in Donselaar et al., 2013) to 12 ( Figure 6). The mapped channel belts formed by deposition of fine to very fine sand and silt in point bars, levees, and as channel-floor deposits (Figures 3, 7, and 8).
Each of the channel belts defined in this study comprises a single-channel sinuous river with SI ranging from 1.25 for channel 7 to 2.56 for channel 2. Bars in the inner bends of the channels show that lateral accretion of fluvial sediment occurred by expansion of lateral bars in the immature low-sinuosity channels, and by expansion and downstream translation and rotation of point bars in welldeveloped meandering channels (cf. Brice, 1974;Bridge, 2006;Daniel, 1971;Jackson, 1976).

| Channel-belt geochronology
The uppermost lacustrine sediments in sample 'Lac' are OSL dated to 12.7 AE 1.5 ka (Table 1, Figures 5 and 11a), which is in accordance with the most recent Coipasa Palaeolake highstand that lasted until 11.5 ka (Baker et al., 2001). This provides a maximum age for the fluvial deposits of the Río Colorado that are directly overlying the lacustrine sediments.
Channel belt 1 is the oldest abandoned river path captured in this study and displays a westward palaeoflow direction ( Figure 6). The intermittent exposure of this channel belt is due to partial reworking by posterior river trajectories. Three OSL samples (1a, 1b, and 1c), taken from preserved river sections, returned similar river-

F I G U R E 7 Visualization on a Google
Earth-Pro image of the formation of alluvial ridges by accumulation of sand in point bars (reddish brown) and crevassesplay and levee deposits (grey, outlined in white). Flow direction towards the northwest (arrow). Lines A-A 0 and B-B 0 correspond to the dGPS profiles in Figure 10 [Color figure can be viewed at wileyonlinelibrary.com] abandonment ages, the youngest of which is 3.63 AE 0.29 ka ('Age' in Table 1 and Figure 11a). As the onset of this river path has not been dated, it is not possible to establish a duration of activity for this channel belt.
Channel-belt complex 2 partly truncates channel belt 1 and flowed to the northwest ( Figure 6). The relative geochronology of the sampled river paths cannot be established from satellite imagery or dating because of the associated OSL age uncertainties. The youngest  Table 2 and interval time in Figure 11b). This does not necessarily represent the lifespan of the channel-belt 2 complex, as there is no evidence of a direct avulsion from channel belts 1 to 2 (i.e. the avulsion point is not exposed, hence the possibility that one or more unidentified channel belts have been active intermediately cannot be excluded).
Channel belt 3 originated from avulsion point A and shows a palaeoflow direction to the north, bypassing the alluvial-ridge geomorphology of channel belts 1 and 2 ( Figure 6). The river path of channel  Figure 11b).
Channel belt 4 shows a palaeoflow direction to the southwest, bypassing previously active river tracts in the northwest quadrant of the study area ( Figure 6) and changing to a westerly direction beyond the study area. The avulsion point from which the river path originated is not exposed because of truncation by posterior river paths 5 and 8, hence it is not possible to constrain the onset of this channel belt. The abandonment age captured in sample 4 is 2.24 AE 0.32 ka (Table 1 and Figure 11a). Note that the mean abandonment age is older than that of channel belt 3. The associated 51% left truncation of the inter-avulsion period indicates a high probability of these channel belts having been (partly) coeval.
Channel belt 5 diverted from channel belt 4 at avulsion point B with a palaeoflow direction to the west ( Figure 6). Upstream of the avulsion point, the new channel developed its own meanders that eroded the channel deposits of its precursor. The abandonment age in sample 5 is 2.26 AE 0.24 ka (Table 1 and Figure 11a), which yields an inter-avulsion period of 0.02 AE 0.40 ka (Table 2 and Figure 11b).
Again, the negative mean and corresponding 52% left truncation are indicative of (partly) concurrent activity. T A B L E 1 OSL dating results (values as best estimate AE standard error). Quality assessment based on assessment of the equivalent dose distribution. Sigmab denotes expected overdispersion for well-bleached samples and differs between samples due to grain size used for analysis and measurement protocol (Cunningham et al., 2011). OD denotes the overdispersion as determined for the sample using the central age model (Galbraith et al., 1999). Unreliable data are in italics. One outlier in sample 10b was removed F I G U R E 1 1 (a) OSL age for each of the dated samples, with 1-sigma standard error. (b) Interval time between consecutive channel belts based on OSL ages. Colour coding of sample numbers as in Figure 6 [Color figure can be viewed at wileyonlinelibrary.com] T A B L E 2 Interval time and metrics for river sinuosity (SI) and alluvial-ridge aggradation. Interval type denoted 'av' (inter-avulsion period), 'max' (maximum value), or 'err' (unreliable data: in italics) Channel belt 6 originated from avulsion point C. The river flowed to the west, then turned southwest ( Figure 6) and again to the west beyond the study area. The abandonment age in sample 6 is 2.35 AE 0.25 ka (Table 1 and Figure 11a). Note that the mean abandonment age is older than that of channel belt 5, hence resulting in a negative interval time and likeliness of coevality.
Channel belt 7 is a narrow (0.19 km wide), poorly developed channel belt formed by a low-sinuosity channel (SI = 1.25) that branched off from meandering channel 6 at avulsion point D, to flow towards the west (Figure 6). Upstream from the avulsion point, the channel has hardly modified the deposits of the former channel.
Sample 7 shows an abandonment age of 1.71 AE 0.18 ka (Table 1 and  8b and 9). The channel has been active for 0.58 AE 0.08 ka (sample 10b; Table 1 and Figure 11a). Sample 11 is indicative of the most recent lateral accretion, dated at 0.53 AE 0.11 ka (Table 1 and Figure 11a).
The present-day river forms channel belt 12 and diverted from channel belt 11 at local avulsion point H ( Figure 6). The river trajectory has an anastomosing pattern, rejoining river path 11 to the north-northwest. Sample 12 captures the first accretion of the channel near its avulsion point H but yields unreliable results (Table 1 and Figure 11a). Beyond avulsion point H, the previous river path 11 is reactivated during peak discharge (Figures 8d and 9), which shows that the avulsion is gradual and ongoing.

| Avulsion frequency, inter-avulsion period, and active migration time
The interval time between successive avulsions in the study area was quantified by using the reconstructed absolute geochronology of channel belts in the study area.
Channel belts 3, 5, 7, and 8 have a clearly defined period between their abandonment age and that of their precursor (Table 2  The active migration rate (Table 3) in the different channel belts (i.e. the point-bar expansion per year) ranges from 0.7 to 2.5 m/y.
Compared with the OSL inter-avulsion timescale, the migration times show that the channels are active only in part of the avulsion period, which is related to the dryland climate setting of the Altiplano Basin, in which point-bar expansion is limited to low-frequency, highmagnitude peak-discharge periods (Donselaar et al., 2013;Li, 2014). In unconfined prograding dryland river systems such as the Río Colorado, the vertical aggradation of channel-lag sand is interpreted to be directly proportional to levee and crevasse-splay sedimentation (van Toorenenburg et al., 2016). As the alluvial-ridge aggradation continued, the channel levees were increasingly perched above the surrounding floodplain (i.e. channel super-elevation). The along-river gradient decreased in the process, whereas the cross-floodplain gradient increased (Figure 10), thus creating a gradient advantage for the river to break out of its course by levee crevassing and switching to a hydrodynamically more favourable, steeper path. The minimum T A B L E 3 Inferred active migration time of the channel belts based on the power-law correlation of (average) channel width and migration rate for unvegetated rivers (figure 2a in Ielpi & Lapôtre, 2020) and point-bar morphometrics (maximum radius). Activity is expressed as a percentage of total OSL-based interval and cannot exceed 100%. Note that only activity-% in bold is based on type 'av' ( aggradation rates to reach the threshold of super-elevation in this study are of the same order as presented by Ganti et al. (2016).

| Gradual avulsion
Several channel belts show a clustering and overlap in OSL sample age. Notably channel belts 3 to 6, 8 and 9, and 10 to 12 (Figure 11a).
In part, the overlap can be explained by the standard error of the sam-  Fielding et al., 2012), but rather as a fluvial system with an anastomosing pattern, caused by the gradual abandonment of the avulsed channels in a low-gradient floodplain setting (Makaske, 2001).
The migration rates of the channel belts in the present study have a smaller range than those presented by Ielpi and Lapôtre (2020; see online Supplementary Information) for the Río Colorado (Table 4). The difference in the upper limit of the migration rate is explained by the smaller amount of data points in the present study, and by the differences in measuring methods of river-channel widths (on-the-ground measurements in this study, and satellite-imagery analysis in Ielpi & Lapôtre, 2020). The OSL sample ages show a time period with a minimum duration in the order of 1 ka between the end of activity of channel belt 6 (at 2.35 AE 0.25 ka) and the start of channel belt 8 (at 1.01 AE 0.17 ka; Figure 11a and Table 1). In this period, the only recorded channel is the small channel belt 7, characterized by a low SI of 1.25 (Table 2), by the development of only a narrow channel belt ( Figure 6) in its long inter-avulsion life of 0.64 AE 0.31 ka (Table 2), and by the short active migration time of 0.112 ka (Table 3). In combination, these data suggest a prolonged drop in river discharge (activity <20% of the total inter-avulsion time, Table 3), possibly related to long-term dry periodicity in the El Niño Southern Oscillation (ENSO) circulation system (Bräuning, 2009). This hypothesis is supported by the interpretation of Baucom and Rigsby (1999) rate of 0.9 m/y, is interpreted as the result of low river discharge related to long-term dry periodicity in the ENSO circulation system.
Multiple avulsions are concentrated in the same area, suggesting a near-nodal avulsion form for the development of the Río Colorado dryland river system. The most recent river paths (channel belts 9-12) suggest an anastomosing pattern for these channels with moderate sinuosity (SI = 1.29-1.75). In contrast, the older large channels are highly sinuous (SI = 1.70-2.56). The change of pattern from highly sinuous channels to anastomosing channels with a lower sinuosity suggests a change of conditions in time, either in the gradient of the area or in river discharge.

CONFLICT OF INTEREST
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available on request from the corresponding author.