Topographic, lithospheric and lithologic controls on the transient landscape evolution after the opening of internally-drained basins. Modelling the North Iberian Neogene drainage

– The opening of internally-drained (endorheic) sedimentary basins often leads to a major drainage change, re-excavation of the basin sedimentary in ﬁ ll, and transient landscape. The timing of such basin openings can be dated only in exceptional cases in which the youngest sedimentary in ﬁ ll remains preserved. For this reason, the processes and timing involved in their transient landscape evolution are poorly known. We explore the role of erodibility, basin geometry and ﬂ exural isostasy during the capture of internally-drained basins by means of numerical modelling techniques constrained by recent terrace cosmogenic dating and geomorphological analysis, addressing the issue as to why the Duero and Ebro rivers, draining two Cenozoic sedimentary basins in N Iberia with similar geographical dimensions and drainage histories, have undergone a markedly different erosion evolution leading to distinctly different present morphology. To evaluate how these intrinsic parameters affect the transient landscape evolution, we design a synthetic scenario inspired by those basins. The results show that, once a basin becomes externally drained, its drainage integration and erosion rates are strongly dependent on (1) the basin elevation above the base level; (2) the width of the topographic barrier, (3) its erodibility; and (4) the rigidity of the lithosphere. The results show that transient landscape evolution can last for tens of millions of years even in absence of tectonic activity and changes in base level or climate. Basins isolated by wide and resistant barriers such as the Duero Basin may undergo a multi-million-year time lag between drainage opening and basin-wide incision. In the case of the Duero Basin, this delay may explain the paradoxical time lag between the last lacustrine bulk sedimentation dated at 9.6 Ma and the onset of widespread basin incision variously estimated at 3.7 to 1 Ma.

Internally-drained (endorheic) basins have no connection to the oceanic base level, evaporating the collected water runoff in central lake systems.Accordingly, they trap all the products of erosion within the continent, often leading to the formation of low relief at high elevation (Han et al., 2019).Their transit to exorheism is driven by either capture, overfilling, or overspilling through a basin outlet located at a topographic sill, the lowest point along the basin's drainage divide (García-Castellanos, 2006).This causes a sudden lowering of the geomorphological base level and an increase in the stream power near the spillway, where the water so far evaporated in endorheic lakes meets the steep slope at the divide with the river basin that takes the role of basin outlet.As the lake's outlet is incised, the lake's level is reduced and, due to its shallow waters, its area becomes rapidly reduced.This is independent on whether the opening is dominated by basin overfilling or by fluvial capture (piracy).The spillway and the just extended drainage basin then tend to adjust to the new base level (Howard et al., 1994) resulting in an upstream propagation of an erosion wave (e.g., Loget and Van Den Driessche, 2009;Struth et al., 2019;Rodríguez-Rodríguez et al., 2020a, 2020b).This retrogressive incision wave erodes the elevated low relief formed during the endorheic conditions.The time-scale of the landscape response to these perturbations is poorly understood (e.g.Whipple, 2001;Whittaker and Boulton, 2012).Previous works analyzed the effect of base level changes (e.g.Bishop, 1995;Burbank and Anderson, 2001), erodibility (e.g.DiBiase and Whipple, 2011), climate (e.g.Schumm, 1979), or elastic thickness (e.g.Chang and Lijun, 2019) in topographic evolution.
Here, we make use of a numerical model integrating all these mechanisms applied to a scenario inspired by the endorheic-exorheic transition of the Duero (Douro in Portuguese) and Ebro drainage basins, covering most of northern Iberia.Remarkably, these two basins have a similar catchment area (∼98,000 km 2 and ∼86,000 km 2 , respectively) and geological history but a completely different present topography.During the Late Miocene, both the Duero and Ebro basins became overfilled with sediment, initiating their incision by a new drainage network.This is based on the mammal fauna (Santisteban et al., 1997) and on magnetostratigraphic dating of the top of the Páramo Formation of the northeastern most Duero Basin (9.7-9.6 Ma according to Krijgsman et al., 1996; ∼9.1 Ma according to Beamud et al., 2006).The same technique combined with flexural isostatic modeling at the Ebro Basin yields an age of opening of 7.5 to 12 Ma (García-Castellanos and Larrasoaña, 2015;García-Castellanos et al., 2003).However, the onset of drainage entrenchment through the Iberian Massif and the Catalan Coastal Range (CCR) barriers (Fig. 1A) seems much slower or delayed by several million years (Cunha et al., 2019a) and their present-day topography and drainage networks are very distinct, with a more eroded, concave-up, and equilibrated Ebro River profile relative to the convex-up Duero River profile (Fig. 1B).The paradox presented by the Duero Basin is that the end of lacustrine sedimentation (9.6 Ma; Krijgsman et al., 1996) in the central areas largely precedes the estimated onset of fluvial incision of the basin (3.7 to 1 Ma; e.g., Silva et al., 2017;Cunha et al., 2019a;Rodríguez-Rodríguez et al., 2020b) and the last minor deposits in the margins of the sedimentary basin (Quaternary alluvial fans and depositional glacis or rañas; e.g., de Vicente et al., 2018;Cunha et al., 2019b).These latest sedimentary infill is considered to be roughly coeval with the onset of basin incision and terrace formation.Recently, however, based on terrace cosmogenic dating, Rodríguez-Rodríguez et al. (2020b) have argued for an earlier onset of basin incision.The time gap between lacustrine sedimentation and the onset of exorheism and the excavation of the sedimentary basin is intriguing because, by definition, internal drainage is a trap for the erosion products of the drained areas and therefore lake sedimentation cannot end sooner than exorheism.
To investigate this seeming conflict, the nature and timing of this transient landscape and the key mechanisms controlling the basin drainage evolution, we use numerical modelling of isostasy and river incision constrained with published geological data, knickpoint analysis, and erosion and sedimentation rates and volumes (e.g.Vacherat et al., 2018;Struth et al., 2019;Cunha et al., 2019a) (Tab.1).We evaluate specifically the effects of geographical (barrier width, basin elevation), lithospheric (elastic thickness; T e ) and lithological (barrier erodibility) parameters on the evolution of the topography under the situation of a high-altitude spill point.

Modelling methodology
We model the topographic evolution using the code TISC, which calculates lithospheric flexural isostasy (elastic thin plate) in response to surface mass transport.This transport is in turn calculated via a fluvial erosion-limited incision and stream-power sediment transport (García-Castellanos, 2002;García-Castellanos and Jiménez-Munt, 2015).The model therefore disregards a significant role of other sources of relative topographic changes such as eustasy of dynamic topography, but isostasy seems to capture the most of the vertical motions in the region (García-Castellanos and Larrasoaña, 2015).The initial model setup is an idealization of the Duero and Ebro basin geology (Figs.1-3A).The input parameters are summarized in Table 1.We adopt a precipitation rate (actually the runoff) linearly increasing with elevation to roughly capture the present-day runoff distribution in the two drainage basins, neglecting past changes in precipitation.The values adopted (Tab. 1) correspond to the average runoff in the region and a conservative evaporation rate that allows for the lakes to cover an area similar to the lacustrine sediment outcrops.In TISC, water flows from cell to cell in a rectangular grid following the steepest descent among eight neighboring cells, forming lakes that fill topographic minima and finding their outlet (if after evaporation the lake keeps extra water, therefore exorheic) or the extent of the endorheic lake (the lake's area is reduced until finding a water balance between its inflows and the lake surface evaporation).The erosion model follows a hybrid detachment-and transportlimited approach in which incision rates are a power law of the basal shear stress of the river flow: where z is the elevation, t is time, k e is the erodability (García-Castellanos and O'Connor, 2018), t c is the shear stress at the base of the water flow, and t c is the critical shear stress.
For simplicity, here we adopt t c = 0 and a = 1.5.This leads to a stream power law for erosion with the form (see García-Castellanos and Jiménez-Munt, 2015): where r is the density of water, g the gravitational acceleration, n the roughness coefficient, k w the proportionality of channel width to discharge.The sediment load q acquired by the river (i.e., the erosion rate in Eq. ( 2)) is limited by a transport capacity q eq , adopted proportional to stream power: Specifically, TISC multiplies the erosion rate in Equation ( 2) by (q eq -q)/q eq in order to cancel erosion and trigger aggradation when the river exceeds its transport capacity.When the river transports a sediment load q larger than q eq then dz/dt becomes positive, implying fluvial aggradation in areas of low channel slope, or delta sedimentation in lakes or in the ocean, where K cap is assumed null.
As an innovation relative to the modeling methodology in García-Castellanos and Jiménez-Munt (2015), we have incorporated here the flow of underground water along the sedimentary units.For this, we solve the Darcy's law in finite difference assuming an isotropic porous medium and using a permeability of 10 À14 m 2 for the sedimentary cover.Remarkably, underground flow multiplies by a factor 2 to 5 the migration velocity of the drainage divide between both sedimentary basins, allowing to reproduce the migration of tens of kilometers observed at the Ebro-Duero divide (Mike s, 2009;Vacherat et al., 2018).

Model setup
The model starts at t = À10 Myr with both basins set very close to overspilling (based on Krijgsman et al., 1996, for the Duero andGarcía-Castellanos et al., 2003, for the Ebro), and ends at 0 Myr (present).We choose our preferred set of parameters (hereafter named reference model) by fitting the knickpoint wave position and retreat rate (Struth et al., 2019;Rodríguez-Rodríguez et al., 2020b), delta sediment volume and basin erosion rates (see Tab. 2).We set the preserved sedimentary infill elevation in the reference model at 800 m above sea level, a compromise between the 535-750 m derived for the Ebro Basin (García-Castellanos and Larrasoaña, 2015) and the 750-850 m at the Duero Basin (Antón et al., 2012, Antón et al., 2014).
This reference scenario adopts an equivalent elastic thickness for the lithosphere of T e = 15 km (representative for the values obtained by Ruiz et al., 2006, Kaban et al., 2018;Gaspar-Escribano et al., 2004), a width of 175 and 25 km for the Duero and Ebro barriers, respectively, with an erodibility of 2.6 10 À8 m yr À1 Pa À1.5 for both bedrock barriers and 4.0e À7 m yr À1 Pa À1.5 for the sediment infill.The former erodibility values, needed to fit the measured erosion rates in the Duero Arribes outlet area, are the lowest ever derived from long-term river erosion, to our knowledge, consistent with the solid, unaltered granite lithology forming most of the highest reaches of the outlet of the Duero sedimentary basin (García-Castellanos and O'Connor, 2018).We calculate mean erosion and sedimentation rates for the basin and the delta area.
Starting from this reference model, we test the effect of varying the lithospheric elastic thickness (T e ), the barrier width and the barrier erodibility.

Modeling results
The results from the reference setup (Figs. 2 and 3; see animation in https://youtu.be/nF0obRQWL60)show that, as the basin's drainage opens finding an outlet towards the ocean, incision gradually propagates into the continent first causing the extinction of the lakes, then the reintegration of the drainage, and then the incision of the sedimentary basins.For large values of elastic thickness T e , the lithosphere behaves too rigidly to induce distinct isostatic vertical motions in response to erosion (Fig. 4A).However, for values lower than 27.5 km (encompassing most T e estimations in Iberia) the isostatic rebound in response to erosion of the barrier and the subsidence related to offshore sedimentation significant cause vertical motions (Fig. 3E).These, in the case of the wider barrier of the western (Duero) basin, further delays in the integration and incision of the lake in the new drainage network.Wider barriers (Fig. 4B) and higher erodibility (Fig. 4C) both extend the time needed to fully capture and extinguish the lake and to entrench and incise the new drainage network.The earliest penetration of the incision wave into the sedimentary basin is thus obtained for higher effective elastic thickness, lower barrier width and higher erodibility.
In spite of the similar initial setting and overspill adopted, the times of maximum basin erosion and delta sedimentation rates (Fig. 5) differ by 7 Myr between the left and right basins of the reference model.The thin layer of sediment prescribed on the Duero barrier and the corresponding spillway allows a fast capture of the basin (a rapid end of endorheism) while the isostatic rebound imposes a slow incision propagation through the less erodible barrier bedrock.Narrow topographic barriers and smaller Te values imply that erosion and sedimentation in the margins induce more isostatic vertical uplift at the edge of the continental basin.The large sediment erosion in the Ebro Basin (unloading of ca.40,000 km 3 , in agreement with estimations in García-Castellanos et al., 2003) results in a similar sediment accumulation in the delta (22% of the sediment is lost across the boundaries of the model) and its corresponding flexural subsidence (Fig. 3).This surface mass transport causes the isostatic rebound of the basin area by up to 600 m (Fig. 3D).This rebound is enhanced by the lesser contribution of the forebulge uplift in response to the sediment loading in the delta (Fig. 3D, E).Interestingly, this mechanism provides an explanation for the exhumation of older sediments in the eastern part of the Ebro basin (see the age of the outcropping sedimentary infill in Fig. 1B; also Lewis et al., 2000).Wider barriers inhibit the forebulge isostatic effect by separating the sediment source from the sink, as apparent in the Duero side in Figures 3D, E and 6.
The barrier width and erodibility also determine the erosion pattern and headwards propagation.On one hand, a wider and less erodible barriers delay the erosive wave retreat towards the basin.The difference in barrier width (175 km for the Duero barrier against 25 km for the Ebro barrier), results in a time lag between the erosion peak at both basins of about 7 Myr (see Fig. 4A).After reaching that peak, erosion rate decreases towards the steady-state in both basins.The results show with two different knickpoint retreat velocities for the bedrock and the cover at 0.02 to 0.03 m yr À1 and 0.10 to 0.20m yr À1 , respectively, close to what is derived from knickpoint analysis by Struth et al. (2019) and cosmogenic dating by Rodríguez-Rodríguez et al. (2020b).This has been tuned by finding the best-fitting ratio of basin-sediment erodability relative to bedrock (a factor 20; Tab. 1).The time needed for  the faster cover-knickpoint wave to reach the central divide (x = 0) is 1.4 Myr in the case of the Ebro and between 7.4-6.0Myr for the Duero.The sooner arrival of the Ebro Basin is due to the narrower topographic barrier, which also explains that both the cover and bedrock knickpoints propagate undistinguishable at the same speed.Importantly, the adopted erodability values are also compatible with the average incision rates obtained for both basins, ensuring that the surface process model is capturing the first order response of both vertical incision and headwards erosion retreat.What ensures that this reference model is eroding at reasonable rates is the present-day (t = 0) average erosion rates of 21.The 4-6 Myr delay (depending mainly on the erodibility of bedrock) between the onset of exorheism and the widespread incision of the basin predicted for the basin on the left side (Fig. 2) is a consequence of the wider hard-rock barrier separating this basin from the ocean.This mechanism also leads to the coeval incision of the center of the basin while small sedimentation rates are still occurring at the basin margins (best visible in the video).This suggests an explanation for the deposition of the Quaternary raña glacis (e.g.Alonso-Zarza et al., 2002) long after the end of bulk sedimentation in most of the Duero Basin.

Discussion on the timing of transient landscape evolution
Internally-drained basins form often amid uplifting orogens that create topographic lows while blocking them  from the income of humid air (e.g.Carroll et al., 2010;Han et al., 2019).The dry conditions allow cancelling the outlet of these basins and prolonging their life (García-Castellanos, 2006).Our modeling results suggest that basins without tectonically active boundaries (or long after the ending of tectonic deformation) can also have their opening delayed by low values of lithosphere rigidity, implying a more local flexural response to the unloading of the basin divides (e.g.Heller et al., 2011).Eroding the boundaries of a closed basin in a low lithospheric rigidity region (normally associated with a young tectono-thermal age of the lithosphere) leads to flexural uplift of the basin flanks, and then to a delay in erosion and reintegration of the internal basin.This mechanism was first reported in the Ebro Basin (García-Castellanos et al., 2003) and the Great Divide Basin (Wyoming, USA) by Heller et al. (2011) and our results show that it is also consistent with the evolution of the Duero Basin and partially explains its different evolution relative to the neighboring Ebro Basin.
The timescales of transient landscape evolution are hard to determine because conditions of paleoclimate and tectonics are often not quantified, and because the identification of a steadystate topography in the field is highly speculative (Willett and Brandon, 2002;Goren, 2016).According to Beeson et al. (2017), the time required to reach the steady-state after a tectonic or climatic perturbation is related to the timescale of river steepness adjustment at individual basins, to balance erosion and rock uplift, defining a transient landscape time.Previous studies (e.g.Paola et al., 1992;Beaumont et al., 2000;Whipple, 2001) suggest that each geomorphic system has an intrinsic response time (or equilibrium time) to approach a balance between the relief-generating mechanisms (e.g., tectonics or isostasy) and the gradual levelling of the surface carried out by surface transport processes.This relates to the time-scales needed after a climatic perturbation, a change in drainage connectivity, or a change in tectonic deformation rates for the average elevation to asymptotically reach stability (Castelltort and Van Den Driessche, 2003) and for sedimentation/erosion rates to remain constant.This response time is often estimated in the order of millions of years (Pazzaglia, 2003, Whipple, 2001, Whipple and Tucker, 1999;García-Castellanos et al., 2003).Based on an analytical formulation, Whipple (2001) estimated a range of 0.25 to 2.5 Myr to adjust to a change in rock uplift rate or base level fall, depending on the non-linearity of the incision rule and the magnitude and type of perturbation.Whittaker and Boulton (2012) calculated a fluvial response time of 1-3 Myr for drainage areas between 10 and 70 km 2 in Hatay Graben (Turkey) and the Apennines.Our results indicate that, even under a hypothetically simultaneous drainage opening with a similar base level fall of 800 m at ∼10 Ma for the two basins, the Ebro-like side reaches the stabilization of erosion after ∼9 Myr (time needed to recover 50% of the erosion rate peak; Fig. 5A) while the Duero-like barrier leads to a much slower and delayed network incision (see stabilization in terms of sediment volume in Fig. 5C).This highlights that the different initial topographic and lithological configuration may explain the differential timing of incision.To see the sensitivity to each parameter, and to allow the portability of our results to other scenarios, we also tested the effect of the drainage area and the initial basin elevation.Doubling the drainage area or precipitation rate leads to a 2.0 and 3.2 Myr earlier entrenchment of the left (Duero) basin in the model, respectively, due to the larger amount of water discharge collected by the main river.Figure 7 shows the result of doubling the base level fall from 800 m (initial elevation of the sedimentary basins in the reference model) to 1600 m, obtaining a nearly double value for the erosion/sediment rate and for the isostatic rebound (Fig. 7).Rates of erosion and sedimentation reach a maximum value peak before adapting to the new base level (Figs.3B and 4A) and then decrease towards a new steady state.Topographic adjustment to the new base level is performed through the propagation of two knickpoint families, a slow one for the bedrock and a fast one along the sedimentary cover.The wider topographic barrier separating the Duero Basin from the ocean implies a time of transitory landscape (time needed to reach a concave-up river profile over the entire basin) of about 30 Myr, 2 to 3 times longer than in the Ebro Basin.
These results support a timescale of landscape transition significantly longer than obtained by the aforementioned studies, up to 30 Myr in our Duero Basin scenario, due to the hard lithology of Variscan granite forming its barrier.In absence of tectonic activity, the transient landscape evolution during the opening of internally-drained basins is mainly controlled by the width and erodibility of the topographic barrier, the size and elevation of the basin above the base level, and the rigidity of the lithosphere.Because erodibility varies by more than eight orders of magnitude (García-Castellanos and O'Connor, 2018), the bedrock lithology along the spillway can exert far more influence on the topographic evolution than the geographical configuration (basin area and elevation).The isostatic vertical motions due to the mass transport during the basin capture are less important, but low values of lithospheric rigidity (typical of regions characterized by young tectonothermal ages such as those affected by the Alpine orogeny) imply a higher amplitude of flexural rebound due to erosion, uplifting the basin margins and delaying the arrival of fluvial incision to the basin, sometimes by millions of years (Fig. 4A), as documented in our study region.

Implications for the drainage evolution of North Iberia
The contrasting preservation of the top of the Duero Basin sedimentary infill versus the advanced erosion of the Ebro Basin can be explained by the wider and harder granite barrier of the former, and does not necessarily imply a much older transition to exorheism of the latter.Our modeling results show that the two basins, in spite of their similar area, base level fall, and age of their sedimentary overfilling, have undergone very different transient landscape evolution due to the wider and less erodible barrier that separates the Duero Basin from the Atlantic Ocean.The time lag between the end of the main Cenozoic sedimentation including lacustrine deposits in the Duero Basin (9.6 Ma; Krijgsman et al., 1996) and the onset of basin incision (possibly starting as early as 5 Ma but surely accelerating after 2.5 Ma; Cunha et al., 2019a;Rodríguez-Rodríguez et al., 2020b) can be explained by the slow propagation of the main knickpoint along the wide, resistant granitic barrier of the Duero Basin (Arribes gorge area).As exemplified by the reference model, a hydrological overspill of the Duero towards the Atlantic Ocean around 9.6 Ma may have terminated the shallow, playa-lake, endorheic system while the granite lithology along the outlet could then prevent the basin incision for several million years.Only when the main knickpoint fully traverses the granitic barrier and reaches the Duero sedimentary basin upstream from Miranda do Douro and Zamora does the incision fully propagate into the center of the basin.This implies that, under conditions of resistant lithology of a wide and flat topographic barrier bounding the endorheic basin, exorheism does not immediate lead to basin incision as often assumed (García-Castellanos, 2006;Cunha et al., 2019a), but in contrast a multi-million year delay can be expected.In other words, while the bulk of incision may have started only recently in geological time-scales (3.7 to 1 Ma), this does not exclude that exorheism may have taken place much earlier, explaining the lack of widespread sedimentary infill in the Duero Basin after 9.6 Ma (except for marginal glacis).
The reference model, constrained by measured river incision rates in both basins, also provides a process-based explanation for the two families of knickpoints (bedrock and sediment cover) retreating at rates of ∼0.02 and ∼0.20 m/yr described by Struth et al. (2019) and Rodríguez-Rodríguez et al. (2020b).In this scenario, the marginal depositional piedmont or glacis (also known as rañas) and fluvial fans next to the ranges can be regarded as the result of pediment formation at the foot of the peri-basinal ranges, in areas where the erosional wave had not yet reached during the Pliocene-Quaternary transition.The Duero Basin, in contrast with other dated basin openings/captures such as the Ebro Basin, provides a geological setting where basin overfilling and drainage overspilling and opening ca.9.6 Ma does not rapidly translate into widespread basin erosion and drainage entrenchment.The Duero basin, due to its wide western topographic barrier made of very hard granitic lithology, underwent a rather delayed acceleration of incision rates only after 2.5 Ma (Rodríguez-Rodríguez, 2019b) or after 1.8 Ma (Cunha et al., 2019b).This scenario can explain the long time elapsed since the end of the main sedimentation in the center of the basin at 9.6 Ma (Krijgsman et al., 1996) until the younger development of incision and terrace formation coeval with pediment or glacis formation in the periphery of the sedimentary basin.It is also consistent with the late transition of the terrace source materials, from recycled basin materials to a bedrock (mountain range) origin, about 1-2 Ma (Rodríguez-Rodríguez et al., 2020b).

Fig. 1 .
Fig. 1. (A) Current topography of the Ebro and Duero basins in the northern part of Iberia.Black dashed line represents the current position of the bedrock incision knickpoint wave in the Duero Basin (Struth et al., 2019).(B) Swath profile of both basins indicating the maximum (red line), mean (green line) and minimum elevation (blue line; river profile).Z max indicates the maximum infill elevation.CCR: Catalan Coastal Range.The x-axis indicates the age of the bedrock and sediment outcropping along the profile.

Fig. 2 .
Fig. 2. Four stages of the evolution of the reference model.Each stage shows the calculated elevation and drainage network (left) and the incision rate distribution (right).The horizontal red lines locate the sections in Figure 3.A video animation of this model can be accessed at https://youtu.be/nF0obRQWL60.

Fig. 3 .
Fig. 3. Reference model results.Panels (A) and (B) show the initial setup at t = À10 Myr.(C) Topography obtained at 0 Myr (present).(D) crosssection showing isostatic subsidence and rebound (dashed line and arrows).(E) Cumulative uplift (dashed) and subsidence (bold lines) and sediment thickness.(F) Sedimentation (À)/erosion (þ) rates at t = 0. River profiles with knickpoint location for the Duero (G) and Ebro (H) river networks.The points displayed in these boxes indicate the location of bedrock (yellow, B_KP) and cover (red, C_KP) knickpoints (compare with Struth et al., 2019).Arrows indicate the distance traveled by the knickpoints from the initial location (yellow circle at x = 50 km).
2 m/Myr (Duero) and 70.1 m/Myr (Ebro), close to the 26 m/Myr measured by Schaller et al. (2016) in the Esla R. tributary of the Duero R. In contrast, the vertical incision along the main river in the reference model (Fig. 3F) is up to 144 m/Myr, in agreement with the 122-250 m/Ma estimated from cosmogenic dating (Rodríguez-Rodríguez et al., 2020b) for the central Duero River.As for the Ebro Basin, the erosion rates in Figure 3F satisfactorily reproduce the 82-127 m/Myr range compiled by Regard et al. (2021).These satisfactory results support the appropriateness of the analytical expression in Equation (3) (details in García-Castellanos and Jiménez-Munt, 2015), and particularly the exponents m = 0.56 and n = 1.22 derived for a = 1.5.

Fig. 4 .
Fig. 4. Model sensitivity test.Age of lake extinction and basin incision onset for varying elastic thickness (A), barrier width (B), and erodibility (C).The grey band indicates the value used for the reference model.t = 0 Myr indicates present day.

Fig. 5 .
Fig. 5. Model results for the Ebro (black lines) and Duero (red lines) basins from À10 to þ20 Myr.(A) Mean erosion rate Er of the two sedimentary basins, showing a trend currently decreasing in the Ebro and increasing in the Duero.(B) Mean sedimentation rates Sr in the delta areas.The inset shows the onset of sedimentation at the Duero delta and synchronous initiation of the incision of the basin.(C) Volume of sediment in the deltas, showing a current stabilization in the Ebro delta contrasting with the increasing volume of the Duero delta (diamonds indicate published data; Table DR-2).

Fig. 6 .
Fig. 6.Distribution of the cumulative sedimentation and erosion (black line) and the isostatic rebound (red dashed line) through the transect displayed in red in Figure 3.

Fig. 7 .
Fig. 7. Results from an identical setup as in Figure 3A but for double elevation of the basins (double base level fall: 1600 m).(A) Topography; (B) cross section; (C) isostatic rebound and subsidence in response to surface mass transport (contours) and sediment thickness; sedimentation/ erosion rates (D) and river profiles with knickpoint location for the Duero (E) and Ebro (F) river networks.Note the more pronounced erosion, integration, and uplift of both basins, relative to the results for the reference model in Figure 3.

Table 1 .
Input parameters for the reference setup of the TISC surface-processes model.

Table 2 .
Results from the reference model and comparison with published observations and results.