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

The study of fluvial network rearrangements provides a key to understand past and future landscape evolution. River profiles respond to tectonic or climatic perturbations on water supply, slope or base level changes (e.g.Whipple and Tucker, 1999). However, the large drainage-area changes triggered by fluvial captures (e.g., García-Castellanos et al., 2003; Willett et al., 2014; Viaplana-Muzas et al., 2018) impose even more abrupt changes in the incision and transport regime. An extreme scenario for such large drainage area changes is provided by the opening of internally-drained basins (endorheic-exorheic transition) by either river piracy, basin overfilling, basin overspilling or a combination of the three mechanisms (García-Castellanos, 2006).

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² and ∼86,000 km², 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.

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⁻¹⁴ m² 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š, 2009; Vacherat et al., 2018).

The model starts at t = −10 Myr with both basins set very close to overspilling (based on Krijgsman et al., 1996, for the Duero and Garcí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⁻⁸ m yr⁻¹ Pa^−1.5 for both bedrock barriers and 4.0e⁻⁷ m yr⁻¹ 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.

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³, 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⁻¹ and 0.10 to 0.20 m yr⁻¹, 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.0 Myr 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.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.

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.

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 steady-state 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² 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.

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).

This work is supported by MITE(CGL2014-59516), SUBTETIS (PIE-CSIC-201830E039) and ALPIMED (PIE-CSIC-201530E082) projects, funded by the Spanish Government. L.R.-R. was funded by the Marie Curie-Clarín COFUND program financed jointly by the 7th WP of the European Union-Marie Curie Actions and Gobierno del Principado de Asturias (Reference ACA-17-19). The work by A.J.-D. was supported by a Juan de la Cierva-Formación postdoctoral contract (ref. FJCI-2016-28878) from the Spanish Ministry of Science, Innovation and Universities. Pedro Cunha, Laurent Husson, and Ciro Cerrone are kindly acknowledged for their detailed and constructive reviews of an earlier version of this manuscript.