Fast Response of Amazon Rivers to Quaternary Climate Cycles

Large alluvial rivers transport water and sediment across continents and shape lowland landscapes. Repeated glacial cycles have dominated Earth's recent climate, but it is unclear whether these rivers are sensitive to such rapid changes. The Amazon River system, the largest and highest‐discharge in the world, features extensive young terraces that demonstrate geologically rapid change temporally correlated with changes in runoff from Quaternary climate cycles. To test the plausibility of a causal relationship, we use a simple model to estimate from empirical measurements how quickly a river profile responds to changes in discharge or sediment supply. Applying this model to data from 30 gauging stations along alluvial rivers throughout the Brazilian Amazon, we find that many rivers of the Amazon basin can respond faster than glacially induced changes in runoff or sediment flux. The Amazon basin is unusually responsive compared to other large river systems due to its high discharge and sediment flux, narrow floodplains, and low slopes. As a result, we predict that the Amazon basin has been highly dynamic during Quaternary glacial cycles, with cyclical aggradation and incision of lowland rivers driving repeated habitat and environmental change throughout the region. This dynamic landscape may have contributed to the exceptional biodiversity of the region and patterns of ancient human settlement.

. Chronology of Quaternary climate and fluvial deposits in the Amazon. (a): One million year records of global benthic foraminifera oxygen isotopes (green curve, left axis), a proxy for ocean temperature and ice volume (Lisiecki & Raymo, 2005), and iron oxide mineralogy from the Amazon fan (goethite/[goethite + hematite], orange curve, right axis) (Harris & Mix, 1999). The iron oxide ratio reflects varying contributions of cratonic and Andean eroded material, which is thought to be determined by climatic factors. (b): Orange and green curves same as (a). Blue points and curve show 50 kyr record of hydrologic change in the eastern Amazon from a speleothem oxygen isotope record (Wang et al., 2017); higher values correspond to drier conditions. Purple, orange, and brown points show 14 C ages from different fluvial terrace levels as shown on block diagram; ages and terrace designations from Rossetti et al. (2014). Vertical arrangement of terrace ages is for visual clarity only, and does not reflect a second variable. Terraces ages and schematic diagram are shown for the Madeira River, which is typical of the Andean rivers in the basin.

The Amazon River Basin
The Amazon River is by far the largest in the world by discharge, five times larger than the next, and carries nearly one-fifth of global river discharge (Milliman & Farnsworth, 2011). The Amazon basin includes the tectonically active Andes mountains, the Brazilian and Guyana cratons, and lowland basins ( Figure 2). Although far from polar ice sheets, it has experienced substantial changes in precipitation during glacial cycles, with more precipitation during interglacial episodes and dry conditions during glacial episodes. This precipitation pattern is attested for during the most recent glacial cycle by speleothem oxygen isotopes (Wang et al., 2017) (Figure 1b) and the preservation of relict aeolian dunes (Carneiro Filho et al., 2002); offshore terrigenous sediments are interpreted to indicate a similar correlation of precipitation and glacial cycles for at least the past million years (Figure 1a) (Harris & Mix, 1999). The same precipitation trend is observed across the Amazon basin from west to east (Cheng et al., 2013;Wang et al., 2017). We use this climate history and the diverse array of rivers in the Amazon basin as a natural experiment to explore how a large, continent-scale river system responds to cyclical environmental forcing.
Unusually for a large lowland river system, most of the rivers of the lowland Amazon basin are entrenched by several tens of meters into alluvial deposits composing the surrounding landscape (Figures 2 and 3), and the high bluffs created by this entrenchment are often the most substantial topographic relief in this generally flat region. This landscape morphology exerts a strong influence on the environmental structure of the region by clearly separating the várzea and igapó (seasonally flooded forests) on the floodplain from the upland terra firme habitat (perennially subaerial) despite typical annual river stage variation of 5-10 m  (Lehner et al., 2008) showing widespread entrenchment of Amazonian rivers into the surrounding landscape. Red circles show locations of sediment flux measurements. Insets A and B show examples of an Andean river (Madeira), a foreland river (Purus), and a cratonic river (Branco), using MERIT DEM elevation data (Yamazaki et al., 2017). Outlined and shaded regions show extent of Cenozoic sediments (Bizzi et al., 2003;Schenk et al., 1999). Dashed lines in A denote intermediate alluvial terraces (Figure 1b). (Meade et al., 1991). This environmental diversity likely plays a role in the exceptional biodiversity of the region, given the ecological differences between the two habitats (Campbell et al., 1986;Gama et al., 2005). This structure may also have fostered ancient agriculture; the largest archaeological settlements with highly fertile anthropogenic soils are often located near the edges of bluffs at the interface of these two environments (Denevan, 1996;Kern et al., 2003;McMichael et al., 2014;WinklerPrins & Aldrich, 2010).
Much of the Amazonian lowlands are composed of recent (Neogene and Quaternary) fluvial sediments ( Figure 2) (Bizzi et al., 2003;Schenk et al., 1999). The entrenchment of modern rivers and floodplains into these sediments indicates past aggradation by fluvial systems, creating the terraces and surfaces preserved today above the modern floodplains, followed by incision of rivers into their former deposits. This switch from aggradation to incision profoundly reshaped the environment. When the present upper terrace adjacent to most Andean and foreland rivers (upper surface in Figure 2a) was an active depositional surface, the landscape was dominated by an extensive, seasonally flooded wetland. The subsequent incision of rivers and their floodplains caused a restriction of the várzea forest to the active floodplain and the emergence of the upland terra firme environment on the relict depositional surface as the areally dominant environment (Pupim et al., 2019). Hypothesized causes of this switch from aggradation to incision include neotectonic crustal deformation (Rossetti, 2014), climatic and hydrologic change (Harris & Mix, 1999;Rigsby et al., 2009), dynamic topography (Shephard et al., 2010), and downstream base level change (Irion et al., 2010). The available sediment burial ages, although sparse, suggest that the upper terrace may in many locations be late Pleistocene in age (Pupim et al., 2019;Rossetti et al., 2014), implying average incision rates since terrace abandonment of up to ∼1 mm/yr. Such rates are geologically rapid, particularly for a lowland region. Intermediate terraces, where found, are typically unpaired and consist of a single level between the upper surface and the active floodplain. Unraveling the influence of the many factors shaping the rivers of the Amazon region has the potential to reveal the origins and history of this unique landscape. To assess the role of periodic climate change in shaping alluvial rivers, and evaluate whether Quaternary glacial cycles could have created the high bluffs and deeply entrenched rivers of the Amazon basin, we use a combination of theoretical analysis, empirical measurements, and topographic analysis to characterize the response of the Amazon river system to cyclical climate change on Milankovitch timescales. We examine how different rivers within the basin responded to Pleistocene climate cycles and associated changes in precipitation, river discharge, and sediment flux, and we compare our model to observed river terraces and sediment burial ages. Our examples span multiple tectonic and geologic contexts, including active orogeny, foreland basin, and cratonic basin.

Model of Alluvial Profile Evolution
We model the evolution of an alluvial river's elevation profile with a non-linear advection-diffusion equation, which we derive from a power-law sediment transport expression using an approach similar to that of previous studies (Begin et al., 1981;Métivier & Gaudemer, 1999;Paola, 2000;Paola et al., 1992;Whipple & Tucker, 2002). We begin with an equation relating total sediment flux to river slope and average water discharge, where Q s is the volumetric sediment flux averaged over the channel cross-section [L 3 /T], K is a transport coefficient [L 3(1−m) T 1−m ], Q w is the water discharge [L 3 /T], and S is the channel slope [-]. m and n are exponents, and are both expected to fall between 1 and 2 (Howard, 1994;Whipple & Tucker, 2002). Similarly to the stream-power law for bedrock river incision (Lague, 2014), a transport law of this form is motivated by a shear-stress based formulation for sediment transport such as that of Engelund and Hansen (1967).
The transport coefficient K includes information about sediment transport by flowing water, physical sediment properties, hydraulic flow resistance, and the relationships between discharge and channel width and depth. As a result, this coefficient should not vary with changing climate or sediment flux. We parameterize the river discharge as a linear function of drainage area (A [L 2 ]), related by the runoff coefficient r [L/T] (Tucker & Bras, 1998): We next invoke an Exner-type volume conservation equation in one dimension (downstream channel distance), where the vertical coordinate z refers to the elevation of the channel [L] and the horizontal coordinate x is the downstream channel distance [L]. ε is the volumetric solids fraction of the bed material [-], Ω is the sinuosity of the river channel relative to the floodplain valley [-], and W is the floodplain width [L]. This model assumes that the channel and floodplain are tightly coupled on the timescales of interest, such that aggradation and degradation of the channel and floodplain are equal. This assumption is justified in the case of Amazonian rivers by the high rate of sediment exchange between the channel and floodplain (Aalto et al., 2003;Dunne et al., 1998). We are therefore modeling the evolution of the combined channel-floodplain system. Combining Equations 1-3 and replacing the channel slope with the negative downstream derivative of channel elevation gives For simplicity, we assume that the runoff and transport coefficients r and K are spatially uniform. Expanding the derivatives in Equation 4 then gives Equation 5 is a non-linear advection-diffusion equation. The first term in parentheses is a non-linear advection term representing changes in sediment transport capacity due to downstream changes in drainage area. The second term in parentheses is a non-linear diffusion term representing changes in sediment transport capacity due to downstream changes in slope. In natural systems, these two effects are of comparable magnitude, and for an equilibrium (graded) profile are necessarily equal and opposite. The effective advection celerity u and diffusivity D are given by We substitute our sediment transport law (Equation 1) to yield simplified forms: The diffusivity given by Equation 9 takes a similar form to that of past studies, which used D = Q s /WS (Castelltort & Van Den Driessche, 2003;Métivier & Gaudemer, 1999). The ratio / s E Q S depends on the discharge: higher discharge means a lower slope for the same sediment flux, which acts to increase the diffusivity. Equation 9 includes three additional terms. The inclusion of the sinuosity Ω more accurately accounts for floodplain geometry; the ratio / Ω E W gives the floodplain area per unit stream length. The packing fraction ε more accurately relates changes in volumetric sediment flux to vertical aggradation or erosion of the floodplain. The slope exponent E n accounts for the non-linearity of diffusion.
The parameters u and D can be used to estimate characteristic response times for the profile evolution of an alluvial river. We assess the relative importance of the diffusive and advective terms by examining the dimensionless Péclet number over an alluvial length scale L: If we substitute a Hack's Law power-law relationship (Hack, 1957) for drainage area as a function of down- where h is the exponent relating drainage area to channel length. We take as our length scale L the length of the alluvial portion of the river. Typical values of m and n range from one to two (Tucker & Bras, 1998;Whipple & Tucker, 2002), with estimated values for the ratio m/n between two-thirds (Howard, 1994) and one (Tucker & Bras, 1998). The Hack exponent h is typically ∼1.7-2. At the river mouth, the ratio L/x will always be less than or equal to one, since the alluvial length can be at most the total length of the river from the drainage divide. Thus, the Péclet number will be close to or less than one. Rivers with Pe ∼ 1, which occurs when L ≈ x and most of the river length is alluvial, have similar advective and diffusive timescales. Rivers with L < x and Pe < 1 have a diffusion timescale that is shorter than the advection timescale, and the system is diffusion-dominated. In either case, the diffusion timescale provides an estimate of the system response time. This is useful because the advection celerity (Equation 8) is less straightforward to quantify for natural systems due to the spatial derivative of drainage area, which is difficult to approximate for rivers with discrete tributaries. We therefore focus on the diffusivity and diffusion timescale as quantitative metrics to characterize and compare the response of alluvial rivers, which we evaluate for natural systems.
With the exception of the slope exponent n, the diffusivity (Equation 9) depends only on physical and empirically measurable quantities: the bed porosity, channel sinuosity, floodplain width, channel slope, and sediment flux. We describe the data and methods used to measure these quantities in Section 2.2. Due to uncertainties in estimating long-term sediment flux from short-term measurements, we do not use the individually measured sediment fluxes to compute the diffusivity. Rather, we regress observed sediment flux against channel slope and water discharge for a number of gauging stations to obtain a transport relationship, which we use to estimate the transport coefficient K in Equation 7. To calculate the diffusion timescale L 2 /D, we use the same length scale as above, the alluvial length of the river. This gives a similar expression to the reaction time of Métivier and Gaudemer (1999), with additional terms of sinuosity (to fully account for the relationship between river length and floodplain area), packing density, and the slope exponent (to account for non-linearity of diffusion). We also use the diffusivities to compute a propagation length δ (analogous to the skin depth in thermal diffusion) indicating the distance an environmental signal will propagate along a river channel for a fixed periodicity of forcing, using the relation    / E DP , where P is the period, which we take to be 10 5 years to simulate recent glacial cycles.

Data Sources
We calculated diffusivities for 30 gauging stations along alluvial rivers in the Amazon using Equation 7. From these diffusivities we calculated diffusion timescales and propagation lengths. We estimated the transport coefficient K from mean annual sediment flux measurements at a network of gauging stations measuring total suspended load over years to decades (Filizola & Guyot, 2009) (Figure 2), converting mass fluxes to volumetric fluxes assuming a sediment density of 2.65 g/cm 3 . While aggradation or incision depend on the total bed material load, including both bedload and suspended load, we expect that most of the sediment in these large, low-gradient rivers is transported as suspended load. We divide these rivers into three tectonic categories: Andes, foreland, and craton. Andes rivers receive sediment directly from the Andes; foreland rivers are those located within the foreland basin but without a direct Andean connection; craton rivers are those draining the Brazilian Shield or Guiana Highlands. We estimated values of K independently for each of the three tectonic settings by regressing the sediment flux measurements against the product of the discharge and slope raised to a power of 1.5, which is in the middle of the published range of values (Tucker & Bras, 1998;Whipple & Tucker, 2002) (Figure 4), with the exponent fixed to facilitate direct comparison of transport coefficients among tectonic settings. We find that the transport coefficients for Andean and foreland rivers are similar, while that of cratonic rivers is about 10 times lower ( Figure 4). This may reflect differences in sediment grain size, bed-load proportion, or channel and flow properties among these different settings, or it may indicate that cratonic rivers are not fully transport-limited and thus not well represented by our model.
We calculated channel slopes using MERIT Hydro flow routing data and MERIT DEM elevation data, at 3-s (∼90 m) resolution (Yamazaki et al., 2017(Yamazaki et al., , 2019, by performing a linear regression on the channel profiles in the vicinity of the gauging stations. We manually measured floodplain width and channel sinuosity from the same datasets. While sinuosity may have varied through time due to changes in discharge and sediment flux, most of the rivers we measured have sinuosities between 1.1 and 1.5, which suggests that potential past variations were small relative to other uncertainties. We assume that the width of the modern floodplain is representative of the width of the floodplain as the river was lowering/incising because this width sets the volume of sediment evacuated by the river to create the steep bluffs bounding the modern floodplain. We used 0.65 for the sediment packing fraction (Dunne et al., 1998). For the length scale L, we used the length of the alluvial portion of each river, as determined from topographic and geologic map data. We use the confluence of the tributary with the mainstem as the downstream boundary.

Numerical Simulations
Although estimated diffusion timescales provide a useful order-of-magnitude assessment of rivers' responsiveness to climate cycles, they only address one of the terms in our advection-diffusion model and do not account for temporal variation in D or u from changing environmental conditions. We tested the validity of this scaling analysis and explored the transient response of rivers with numerical solutions. We solved Equation 5 using a centered-space, forward-time finite difference method and a Hack's law approximation for drainage area (Hack, 1957). We calibrated the model parameters to two example Amazon tributaries, the Madeira and Guapore Rivers, representing the Andean and cratonic settings respectively (Table S1). Both of these rivers feature clear alluvial features such as meanders and oxbows within floodplains. We use m = n = 1.5 and a Hack exponent h = 1.8. We used the transport coefficients for Andean and cratonic rivers (Figure 4), respectively, and discharge from the Porto Velho and Pimenteiras gauging stations (Filizola & Guyot, 2009). For the initial condition, the upstream boundary was assigned discharge, drainage area, and slope from the gauging station, with sediment flux given by the transport coefficient and transport law. For the initial profile we used a graded (equilibrium) profile for the modern value of the runoff, at which the transport capacity in Equation 1 is equal to the upstream supply at all points in the model domain. We imposed the upstream sediment supply as a slope boundary condition, and constant base level at the downstream boundary. Although these rivers were likely subject to base-level change from ice age sea-level effects, our theoretical model is generally agnostic to the particular mechanism of perturbation.
We present two sets of simulations of oscillatory climate forcing. In the first, we varied the runoff coefficient sinusoidally by a factor of two with a 100 kyr period, similar to precipitation changes in the Amazon basin during recent glacial cycles (Figure 1), while keeping sediment supply constant. In the second, we varied the sediment supply by a factor of two under a constant runoff. Quaternary climate cycles are generally characterized by asymmetric "sawtooth" patterns, with gradual glaciation and rapid deglaciation, in contrast to the symmetric sinusoid in our simulations (Lisiecki & Raymo, 2005). Nevertheless, the 100 kyr periodicity generally reflects the time between successive maxima or minima in climatic parameters, and thus the simulated response time should characterize the degree to which the system can adjust to a new equilibrium before the next cycle. These simulations are not meant to be quantitative reconstructions of the history of the Madeira and Guapore systems, but rather attempts to understand how rivers with similar characteristics might evolve under a cyclical climate such as Earth has recently experienced. We seek Figure 4. Estimating the sediment transport coefficient. Power-law regression of mean annual sediment flux against the discharge-slope product raised to the 3/2 power. The intercept gives the transport coefficient K. We perform a separate regression for each of the three tectonic settings.
to understand a river's ability to respond on relevant timescales to any environmental change that causes profile disequilibrium.

River Response Times and Propagation Lengths
The estimated diffusivities span four orders of magnitude, from 1.6 × 10 5 m 2 /yr to 8.5 × 10 8 m 2 /yr (Figure 5a). The diffusivities show a strong dependence on tectonic setting, with the Andean rivers at the high end of the range (geometric mean 3.8 × 10 8 m 2 /yr), the cratonic rivers at the low end (1.2 × 10 6 m 2 /yr), and the foreland rivers between (1.2 × 10 7 m 2 /yr). This is mainly due to the large variations in discharge and sediment flux, although this effect is partially counteracted by the floodplain width in the denominator of Equation 2, which scales with sediment flux and drainage area. The diffusivity should generally increase downstream with increasing drainage area, discharge, and sediment flux, although the rivers we sample are generally sufficiently large and the gauging stations sufficiently far downstream that our estimates are broadly representative of the lower reaches of the rivers; we do not expect that the diffusivity increases substantially downstream of the gauging stations in our study reaches.
When we convert the diffusivities to diffusive response times (Figure 5b), the range is reduced to three orders of magnitude, with the shortest response times belonging to the Andean rivers, with a minimum of 13 kyr and a geometric mean of 36 kyr, and the longest response timescales belonging to the cratonic rivers, with a geometric mean of 340 kyr and a maximum of 9 Myr. Foreland rivers have intermediate response times with a geometric mean of 160 kyr. The reduction in range relative to diffusivities is due to the compensating effect of river length-the Andean rivers with the largest discharge, sediment flux, and diffusivity are also substantially longer than most of the foreland and cratonic rivers. Nevertheless, the short response times of the Andean rivers show that the high discharge and sediment flux outweigh their great length, even in the case of the mainstem Amazon, which has nearly 5,000 km of alluvial length. The propagation lengths for a 100 kyr period (Figure 5c) have geometric means of 4,000, 820, and 230 km for Andean, foreland, and cratonic rivers respectively. While these are likely order-of-magnitude estimates due to uncertainties in the parameters, the response times for Andean rivers (Figure 5b) are all shorter than the 100 kyr period of recent glacial cycles, and most are substantially shorter, by a factor of 2 to 6. Similarly, the propagation lengths for Andean rivers (Figure 5c) all exceed 1,000 km, indicating that environmental perturbations with a 100,000 year period can elicit a response over continental scales. Comparing the propagation length with the actual alluvial river length is analogous to comparing the response timescale with the periodicity of forcing: rivers with propagation lengths longer than their alluvial lengths have response timescales shorter than the forcing period, and vice versa.
The response times for cratonic and foreland rivers are more varied and span a wide range on both sides of this glacial timescale. Foreland rivers have response times with geometric and arithmetic means of 160 and 340 kyr and propagation lengths with geometric and arithmetic means of 820 and 1,050 km. Cratonic rivers have response times with geometric and arithmetic means of 310 and 700 kyr and propagation lengths with geometric and arithmetic means of 230 and 290 km. As we discuss in more detail in Section 4, these estimated response times and propagation lengths suggest that the Andean rivers and many foreland rivers have repeatedly aggraded and incised in response to Quaternary glacial cycles.

Numerical Results
The simulation results are consistent with our analysis of diffusive response times. The two river systems we used to calibrate the numerical simulations have estimated response times of 14 kyr (Madeira) and 1.3 Myr (Guapore), a difference that reflects their vastly different sediment loads and discharges (Table S1). Correspondingly, in the Madeira simulation, the modeled river profile rapidly adjusts to the changes in runoff or sediment flux and is always near its equilibrium profile throughout the simulation (Figure 6b). This is consistent with the Madeira River's estimated response time, which is substantially shorter than the periodicity of the modeled forcing. Analogously, our estimated propagation length of 5,000 km is longer than the 3,000 km alluvial length of the Madeira River.
The Guapore simulation behaves very differently (Figure 6c). Since the response time is much longer than the 100 kyr period of forcing, the river profile is unable to adjust to each cycle (green curves in Figure 6c). The distance of profile adjustment from one cycle is comparable to the ∼200 km propagation length we estimate. Instead, over many cycles the river approaches a profile governed by the long-term mean runoff (purple curves in Figure 6c), with the individual cycles damped by the long response time. The diverging behavior of these two simulations is a direct consequence of their different response times relative to the periodicity of the imposed climate forcing; the diffusive response time accurately captures the behavior of the advection-diffusion system we model. This dichotomous behavior occurs similarly for a periodic sediment supply with the same period ( Figure 6). A similar relationship between channel adjustment and periodicity of climatic forcing has been shown for bedrock rivers within a numerical landscape evolution model, where the response time depends instead on the erodibility (Godard et al., 2013).

Response of Amazonian Rivers to Cyclical Climate Change
The high diffusivities, short response times, and long propagation lengths for the Andean rivers of the Amazon basin are the result of their high discharge and sediment flux. These characteristics enhance Andean rivers' ability to respond to environmental change by aggradation and incision, which depend on a river's capacity to deposit or erode alluvial sediment. Planform changes in alluvial rivers also depend on the entrainment and deposition of sediment; the river response times estimated here may also be useful for understanding changes in meandering in response to climate change, and the meander migration rate of Amazonian rivers similarly has been found to vary as a positive function of sediment flux (Constantine et al., 2014). Our estimated response times of Andean-draining rivers are consistent with a floodplain recycling time of a few thousand years estimated from sediment measurements along the Amazon mainstem (Mertes et al., 1996), which further supports the applicability of an advection-diffusion model for these rivers, and likely for the similar foreland rivers as well.
The dependence of alluvial planform morphology on sediment entrainment and deposition also implies that floodplain width is not necessarily constant with time. Terrace abandonment can narrow floodplains, and lateral bluff erosion can widen them. In the case of the most recent incision event, the modern floodplain width is representative of the sediment volume eroded and evacuated since the switch from aggradation to incision and is thus the appropriate width to use when assessing the response to climate changes in the past 30 kyr. Floodplain width may have been different during earlier glacial cycles, which could have affected the diffusivity, but subsequent aggradation has buried these ancient floodplains, leaving no straightforward way to estimate their widths.
Although we model the effect of changes in runoff and sediment supply separately, they likely have both changed simultaneously in response to climate cycles and associated changes in vegetation, and may not have been in phase (Perron, 2017). The response of river systems to climate cycles depends on the net effect of changes in discharge and sediment supply on the imbalance between sediment supply and transport capacity. However, as long as these changes are governed by the same global climate cycles and exhibit a similar periodicity, a river profile's ability to adjust to any combination of external changes in discharge and  (c, d, g, and h). (c, d, g, h): Simulated river profiles using parameters for the Madeira River (c and g) and the Guaporé River (d and h). The gray line in each plot is the initial condition. For each colored profile, the solid line shows the actual computed river profile, and the darker dashed line shows the equilibrium graded profile for the instantaneous values of the runoff and sediment flux. Line colors correspond to the times marked by colored bars in a, b, e, f. sediment supply depends on the timescale of those changes relative to the river's own response timescale. In the Amazon region, the correlation between dry conditions and terrace deposition ages (Figure 1) (Pupim et al., 2019;Wang et al., 2017) demonstrates that the net effect of climate-induced changes in runoff and sediment supply in this region was a reduction in transport capacity relative to supply during dry periods, causing aggradation, and excess transport capacity during wet periods, causing incision and terrace abandonment. This suggests that changes in runoff were the dominant control on fluvial terrace evolution in the region.
We also note that, although the diffusivity and thus response timescale vary with changes in runoff and sediment supply due to non-linearities in Equation 5 and in the underlying physical processes, these changes are likely modest. A two-fold change in runoff, with all other parameters constant, changes the diffusivity and diffusion time by a factor of / 2 m n E , which is approximately a factor of 2 since m/n ∼ 1. A two-fold change in sediment supply, all other parameters constant, changes the diffusivity and diffusion time by a factor of  1 1/ 2 n E , which is less than a factor of 2 since the exponent is less than one. As a result, the relationship between a river's response time and Quaternary climate periodicity is unlikely to be substantially altered.
Although a river with a fast response time can adjust its profile in response to changes in either runoff or sediment supply, changes to these two forcings have opposite effects on sediment output for a river system that evolves according to the advection-diffusion model. In the case of runoff variations with constant sediment input, a river with a slow response time has a variable sediment output whereas a river with fast response time maintains a constant output equal to the input. Conversely, in the case of sediment input variations with constant runoff, a river with a slow response time maintains a constant sediment output whereas a river with a fast response time has variable output corresponding to the variable input. This is because a river's ability or inability to adjust its longitudinal profile and slope governs sediment transport capacity under constant discharge. These differences in sediment output between different climatic forcings have important implications for the transmission and preservation of climatic signals to the sedimentary record, since the two end-member forcings are processed differently in the alluvial river transfer system (Malatesta et al., 2018;Romans et al., 2016;Simpson & Castelltort, 2012;Straub et al., 2020). Comparison of the onshore terrace record with high-resolution offshore sediment accumulation data could thus discriminate between changes in runoff and sediment supply as drivers of geomorphic change, as has been shown experimentally by Tofelde et al. (2019)

Quaternary Landscape Evolution of the Amazon Basin
We extend our analysis of river response times ( Figure 5) to infer the recent landscape history of the Amazon region, since speleothem records indicate broadly consistent trends of wetting and drying across the basin (Cheng et al., 2013;Wang et al., 2017). The sediment-rich Andean-sourced rivers of the western and central Amazon have the capacity to quickly respond to changes in climate, and as a result have probably repeatedly aggraded during dry glacial periods and incised into these deposits during wet interglacials. This inference is supported by the short response time (∼30-60 kyr) of the mainstem Amazon, despite its extraordinary alluvial length (4,900 km), and by similarly short response times for other Andean tributaries. The deposits that today make up the terra firme regions are likely a palimpsest of past glacial episodes, formed and reworked during many successive aggradational phases. This may explain some of the anomalously old ages measured from OSL dating of upper terrace deposits in close proximity to terraces with young MIS 2-3 ages (Cremon et al., 2015): the older terraces could have formed during previous aggradational episodes without subsequent disturbance.
The spatial coherence of the upper surface is inconsistent with autogenic terrace formation, and the rapid rates of incision required by the young depositional terrace ages are too fast to have resulted from far-field uplift from tectonics or dynamic topography, although these processes may also be occurring on longer timescales. The deposition of these terrace deposits when sea level was substantially lower than present further excludes downstream sea-level change as the primary driver of terrace abandonment and incision. The discontinuous intermediate terraces, found between the modern floodplain and the upper surface, may result from autogenic dynamics of lateral channel migration during incision following the increase in runoff (Hajek & Straub, 2017;Limaye & Lamb, 2016). We note that the age of the uppermost aggradational surface does not directly indicate the periodicity of the environmental change that caused its formation, since river terraces only record one part of the response cycle (terrace abandonment and incision). The depositional ages indicate the time since the most recent peak in the response cycle (i.e., when aggradation transitions to incision), but not the time until the next equivalent peak.
The similar pattern of terraces in many foreland rivers that do not directly receive Andean sediment (Figures 2 and 3) suggests that many of these rivers behave similarly to a degree. Although their diffusivities are lower due to smaller discharges and drainage area, this is partially compensated by their shorter lengths than Andean rivers. As a result, they mostly have response times between 100 and 200 kyr. We note that a response time comparable to or even slightly longer than the periodicity of forcing still permits a detectable response to cyclical change, even if the system does not have sufficient time to fully reach a new equilibrium with each episode of wetting or drying. This is demonstrated by the propagation lengths of many hundreds of kilometers for foreland rivers for a 100 kyr periodic forcing, indicating that an observable response would be felt through much of the river system. In numerical experiments not shown in Figure 6, we observed that a river adjusts approximately halfway to a new equilibrium profile at a time of one-half the response time following a step change in discharge, gradually decaying exponentially towards its new equilibrium as time goes on.
The cyclical aggradation and incision of Andes-draining and foreland rivers has created a highly dynamic physical environment, with seasonally flooded wetlands counterintuitively dominating the region during dry glacial phases due to aggradation, and the present condition of extensive terra firme uplands following incision and localization of the floodplain during wetter interglacials. The low-sediment flux cratonic rivers of the eastern Amazon may have evolved more slowly, buffered from climate changes on glacial timescales. This may explain why many cratonic rivers in the region do not feature major terraces, indicating that they have likely not experienced significant aggradation and incision during Quaternary climatic oscillations. Such rivers often have terraces near their downstream ends driven by base-level changes from the dynamic large rivers they flow into, but these terraces diminish upstream, reflecting the short propagation length. The difference in terrace occurrence between foreland and cratonic rivers is somewhat at odds with the overlapping distributions of estimated diffusion times for cratonic and foreland rivers ( Figure 5). The systematically lower sediment flux of cratonic rivers (Figure 4) may indicate that these rivers are not fully transport-limited, even though we selected gauging stations that appeared qualitatively alluvial from topographic data and satellite imagery. If this is the case, then the transport-limited model we use in the study may not be a good description of their behavior. Partly bedrock cratonic rivers would likely imply an even slower fluvial response to environmental change than predicted by the alluvial model, which would be consistent with the less abundant and less pronounced terraces in the craton. The differing fluvial responses among different tributary basins may have contributed to regional drainage reorganization, although this effect was likely limited by the spatial correlation of response times-since response time is largely determined by tectonic setting, rivers of comparable responsivities are spatially clustered.

Amazon in a Global Context
Contrary to our findings, a past study (Métivier & Gaudemer, 1999) found evidence for buffered behavior on glacial timescales for the large Himalayan rivers of the Indo-Gangetic Plain, using a similar model. A global analysis (Castelltort & Van Den Driessche, 2003) argued that large river systems are typically buffered on these timescales. Several important differences explain the high reactivity of Amazonian rivers relative to the Himalayan rivers and many other large rivers worldwide. Rivers of the Amazon carry comparable sediment loads to those of the Himalaya, yet have floodplains an order of magnitude narrower and slopes an order of magnitude shallower than the braided Himalayan rivers. Both of these factors increase the effective diffusivity and decrease the response time of the Amazon rivers (Equation 9), making them unusually responsive for a river system of their magnitude. These differences are at least partly due to the exceptional discharge of the Amazon, owing to its location in the humid tropics and its high rainfall and runoff throughout the basin. The diffusivity depends on the ratio Q s /S, which depends in turn on the discharge-the large discharge of the Amazon allows the river to carry its sediment load at very low slopes (∼10 −5 ). The single-threaded rivers of the Amazon thus have a set of characteristics-narrow floodplains, high discharge and sediment loads-that predisposes them to adjust rapidly to environmental changes and fosters the dynamic landscape of the western and central Amazon. We note that although Castelltort and Van Den Driessche (2003) used a similar mathematical framework to estimate response times for a set of global rivers, their analysis differed substantially from ours by using the channel width instead of the floodplain width to estimate the diffusivity, and by using the entire river relief and length including bedrock and headwater reaches. This implicitly assumes that the buffering ability comes only from the channel itself, without aggradation or degradation of the adjacent floodplain, which underestimates the response times of rivers with active floodplains. Additionally, their analysis did not differentiate between alluvial and bedrock river reaches; the latter is not well-represented by a diffusion law such as we use here. The inclusion of the entire length and relief of the river overestimates the response time of a river. The combination of these assumptions likely overestimates the response time; in the case of the Amazon, the authors compute a response time of 104,000 years.

Implications for Amazonian Biodiversity and Human Settlement
The Amazon region is one of the largest and most species-rich biodiversity hotspots on Earth (Hoorn et al., 2010), and a conspicuous global exception to the correlation of high biodiversity with high-relief mountain ranges (Antonelli et al., 2018). One proposed explanation for the exceptional terrestrial biodiversity of the Amazon is the repeated expansion and contraction of the rainforest habitat due to ecological response to climate changes, with isolation in forest-fragment refugia during glacial episodes driving speciation and diversification (Haffer, 1969). We speculate that the physical changes to the landscape caused by fluvial aggradation and incision and the shifting distribution of várzea and terra firme habitats may have played a similar role. During dry aggradational episodes, seasonally flooded wetlands were dominant and upland regions restricted. Conversely, during wet incisional episodes, upland regions became extensive. The onset of Pleistocene glaciation and the intensification of glacial cycles at the Mid-Pleistocene Transition may have thus been drivers of speciation in the region through repeated habitat gain and loss. Our analysis suggests that the shift from 40 to 100 kyr periodicity would have caused a larger portion of the Amazon basin to begin to respond to these climatic changes. Moreover, the amplitude of glaciations increased with this transition, which likely increased the variability of tropical rainfall as well. Indeed, phylogenetic evidence indicates that these climate transitions correspond to elevated origination rates for both vertebrate and invertebrate lineages in the Amazon basin (Garzón-Orduña et al., 2014;Ribas et al., 2012).
The widespread terraces created by the most recent incision event also set the stage for human settlement of the Amazon. Precolonial settlements, as well as many modern settlements, are typically located along the bluffs created by these terraces-on upland terra firme, above and adjacent to rivers and floodplains (Denevan, 1996). The presence of extensive uplands that do not flood, unusual for a large lowland river system, likely contributed to the emergence of the Amazon as a culturally influential region (Heckenberger, 2013) and a center of crop domestication (Clement et al., 2015).

Conclusions
To investigate whether the large alluvial rivers of the Amazon basin have responded to Quaternary glacial cycles, we estimated diffusive response timescales for a range of Amazon tributaries using a simple theoretical model and empirical measurements. The timescales we calculate, which we support with numerical simulations, imply that Andean and foreland rivers have repeatedly aggraded and incised in response to Quaternary glacial cycles, which agrees with the observed wide range of terrace ages in the Amazon lowlands, including many recent dates (25-50 ka). The dynamic behavior of such a large river system, which differs from many other large rivers worldwide, stems from the tectonic and climatic setting of the basin, which give it abundant sediment and water discharge, as well as the narrow floodplains of Amazon rivers, which give it the power to respond rapidly to environmental change on glacial timescales and create a highly dynamic landscape under a cyclical climate.

Data Availability Statement
Supplementary data are given in the Table S1 and archived at https://doi.org/10.5281/zenodo.5600823.