Abstract

The Ebro Basin constitutes the central part of the southern foreland of the Pyrenees. It was endorheic during the Cenozoic and accumulated sediments. By the end of the Miocene, erosion and river incision reconnected the basin to the Mediterranean Sea, establishing a post-opening drainage network. Those rivers left terraces that we study in this work. We first synthesize previous works on river terraces that are widely dispersed in the basin. We provide new age constraints, up to 3 Ma, obtained thanks to cosmogenic nuclides using both profile and burial methods. We derive a unified fluvial terrace chronology and a homogenized map of the highest terraces over the entire Ebro Basin. The dated terraces labeled A, B, C, D, and E are dated to 2.8 ± 0.7 Ma, 1.15 ± 0.15 Ma, 850 ± 70 ka, 650 ± 130 ka, and 400 ± 120 ka, respectively. The chronology proposed here is similar to other sequences of river terraces dated in the Iberian Peninsula, around the Pyrenees, and elsewhere in Europe. The oldest terraces (A, B, C) are extensive, indicating they form a mobile fluvial network while from D to present, the network was stable and entrenched in 100 to 200 m-deep valleys. The transition from mobile to fixed fluvial network is likely to have occurred during the Middle Pleistocene Transition (MPT, between 0.7 and 1.3 Ma), when long-period/high-intensity climate fluctuations were established in Europe. We estimate that between 2.8–1.15 Ma and present, the incision rates have tripled.

Résumé

Le bassin de l’Ebre constitue le centre du piémont sud-pyrénéen. Il a été endorhéique et a accumulé une épaisseur importante de sédiments au cours du Cénozoique. Vers la fin du Miocène, l’érosion et l’incision fluviale ont reconnecté le bassin à la mer Méditerranée, en établissant un nouveau réseau de drainage. Les rivières constituant ce réseau ont laissé des terrasses que nous étudions dans le cadre de ce travail. Nous faisons d’abord une synthèse des travaux antérieurs sur les terrasses fluviales, qui sont largement dispersées dans le bassin. Ensuite, nous fournissons de nouvelles contraintes d’âge, jusqu’à 3 Ma, obtenues grâce aux nucléides cosmogéniques utilisant à la fois des méthodes de profil d’exposition et d’enfouissement. Enfin, nous en déduisons une chronologie unifiée des terrasses fluviales et une carte homogénéisée des terrasses les plus hautes sur l’ensemble du bassin de l’Ebre. Les terrasses, que nous notons A, B, C, D et E, sont datées respectivement de 2,8 ± 0,7 Ma, 1,15 ± 0,15 Ma, 850 ± 70 ka, 650 ± 130 ka et 400 ± 120 ka. La chronologie proposée ici est similaire à d’autres séquences de terrasses fluviales datées dans la péninsule ibérique, autour des Pyrénées, et ailleurs en Europe. Les terrasses les plus anciennes (A, B, C) sont étendues, ce qui témoignent d’un réseau fluvial mobile alors que depuis la formation des terrasses D jusqu’à aujourd’hui, le réseau était stable et canalisé le long de vallées de 100 à 200 m de profondeur. Le passage d’un réseau fluvial mobile à un réseau fluvial fixe s’est probablement produit pendant la transition du Pléistocène moyen (entre 0,7 et 1,3 Ma), lorsque des fluctuations climatiques de longue période et de forte intensité ont été établies pour les rivières européennes. Nous estimons qu’entre 2,8 et 1,15 Ma et la période actuelle, les taux d’incision ont triplé.

Introduction

Intracontinental endorheic basins (or closed basins) occupy ∼ 20% of the Earth’s land surface. They develop mostly in response to the tectonic uplift of surrounding ranges, generally under arid climatic conditions (e.g.Sobel et al., 2003; García-Castellanos et al., 2003) and are infilled by low-erodible sedimentary rocks. They represent key elements of source-to-sink systems as they preserve eroded products from surrounding catchments and large areas, as for the Tibetan Plateau and Tarim basin (Sobel et al., 2003; Han et al., 2019) or the Altiplano (Fornari et al., 2001; Sobel et al., 2003). Characterized by internally-drained depressions they often record the development of a long-lived lake (Altiplano, Titicaca Lake). The opening of an endorheic basin by the capture of internally-drained rivers and lakes has strong implications on the sediment routing and the evolution of the river landscape (Craddock et al., 2010). Despite the development of analogical and numerical modelling tools, processes by which the river network grows and propagates into a captured endorheic basin remain unclear.

The triangular-shaped Ebro sedimentary Basin in the north-eastern Iberian Peninsula is delimited by the South Pyrenean thrust belt to the north, by the Catalan Coastal Range to the south-east, and by the Iberian Range to the south-west (Fig. 1). The basin recorded a long endorheic stage during Oligocene and Miocene times (Puigdefàbregas et al., 1992; Costa et al., 2010). The re-opening of the Ebro Basin toward the Mediterranean Sea led to an incision wave that propagated within the former endorheic basin and the Catalan Coastal Range (e.g.García-Castellanos et al., 2003). The drainage divide delimiting the current Ebro catchment shows numerous evidence of geomorphic disequilibrium suggesting that it is still currently enlarging through drainage network growth, captures, and divide migration (Vacherat et al., 2018; Struth et al., 2019).

In this study, we aim to provide a better resolution of the timing of incision in the Ebro catchment and improving our understanding of the relationships between the morphologic evolution of the Ebro drainage network and the longer-term post-orogenic evolution of the Ebro sedimentary Basin. We first review constraints on Ebro river terrace identifications and ages to propose a consistent terrace chronology scheme throughout the basin. Then, we present new dates for the topmost terraces using in-situ produced cosmogenic nuclides (10Be and 26Al), hence providing constraints on the most older, i.e. initial, recorded stages of landscape development of the Ebro river network. These new data allow us to track the evolution of the Ebro river network back to the Pliocene.

Geological setting

Tectonic setting

The Pyrenean orogen resulted from the collision between the Iberian and European continental margins that occurred from the Late Cretaceous to the early Miocene (Choukroune, 1989; Roure et al., 1989; Teixell, 1990, 1996, 1998; Muñoz, 1992; Beaumont et al., 2000; Vergés et al., 2002; Mouthereau et al., 2014). The main period of orogenic growth and exhumation is dated to Eocene-Oligocene (45–30 Ma) (see Yelland, 1990; Morris et al., 1998; Sinclair et al., 2005; Gibson et al., 2007; Jolivet et al., 2007; Rushlow et al., 2013; Labaume et al., 2016). The South Pyrenees thrust sequence shows that deformation propagated from east to west and from north to south since the Late Cretaceous. This orogenic sequence is reflected, north of the Ebro Basin, in the syn-orogenic tectonic-stratigraphic relationships well preserved in the South Pyrenean Zone (SPZ), a classical thin-skin fold and thrust belt made of Mesozoic to Cenozoic sediments detached in the Triassic evaporites (Muñoz, 1992; Vergés et al., 1995, 2002; Mouthereau et al., 2014; Carola et al., 2015; Saura et al., 2016). It was followed by a post-orogenic phase of exhumation reported throughout the Pyrenees at about 10 Ma (Morris et al., 1998; Fitzgerald et al., 1999; Whitchurch et al., 2011; Fillon et al., 2013, 2020; Mouthereau et al., 2014; Bosch et al., 2016; Monod et al., 2016; Huyghe et al., 2020; see also the review by Calvet et al., 2020).

The Iberian chain, south of the Ebro Basin, forms a NW-SE trending doubly vergent intraplate mountain belt. It results from the tectonic inversion of Mesozoic rift basins comprising from east to west the Columbrets, Asturian, Basque-Cantabrian, Maestrat, south Iberian, and Cameros basins (Roca and Guimerà, 1992; Salas and Casas, 1993; Salas et al., 2001; Etheve et al., 2018). Constraints on the recent exhumation pattern of the Iberian chain mainly come from the Cameros basin (Del Rio et al., 2009; Rat et al., 2019). ZFT thermochronological constraints show a minimum age for the inversion initiation of the Cameros basin at ∼ 60 Ma, which is coincident or slightly predates the widespread increase in exhumation rates during the Eocene documented in the Pyrenees. The main exhumation phase is recorded in the Cameros basin at 35–25 Ma coevally with the Pyrenees and the Cantabrian belt (Fitzgerald et al., 1999; Fillon et al., 2016). The increase of exhumation is temporally consistent with the closure of the connection of the Ebro Basin with the Atlantic Ocean by 36 Ma, leading to the deposition of thick alluvial sediments (Costa et al., 2010). As in the Pyrenees, this exhumation was followed by a Miocene period of planation that resulted in the paleosurfaces still visible on top of the Iberian chain (Calvet et al., 2015). The deposition of late Miocene conglomerates, sealing the main Cameros thrust, was followed by erosion after ∼ 9 Ma (Rat et al., 2019). This result is in agreement with the regional drainage reorganization described throughout the Ebro Basin (García-Castellanos et al., 2003; Fillon et al., 2013, 2020).

The Catalan Coastal Ranges (CCR), located between the Ebro Basin and the Valencia Trough, formed by the tectonic inversion of Mesozoic rift basins during the Paleogene, contemporaneously with the Pyrenean orogeny (Salas et al., 2001; Gaspar-Escribano et al., 2004). From Late Oligocene onward, extension in the Valencia Trough associated with opening of the Gulf of Lion led to the development of ENE-WSW to NE-SW-striking horst and graben, triggering uplift (Lewis et al., 2000; García-Castellanos et al., 2003).

During the Late Eocene, uplift in the Western Pyrenees (Muñoz et al., 1986; Puigdefàbregas et al., 1992) closed the Ebro Basin. This resulted in a sudden increase of sedimentation rates and continentalization of the endorheic basin (Costa et al., 2010). The deposition of a large amount of coarse alluvial deposits (Coney et al., 1996; Babault et al., 2005b) from Bartonian to Late Oligocene (Beamud et al., 2003, 2011) filled the paleovalleys that developed at the front of the South-Pyrenees (Vincent, 2001). Such conglomeratic bodies are also observed against the Catalan Coastal Range and the Iberian Range. The center of the Ebro basin became a long-lived lake filled with lacustrine and sandy deposits during the Oligocene and the Miocene (Pérez-Rivarés et al., 2002, 2004; García-Castellanos et al., 2003; García-Castellanos, 2006; Larrasoaña et al., 2006; Vázquez-Urbez et al., 2013).

In the Pyrenees, in addition to the low temperature thermochronology data already mentioned, post-10 Ma uplift is suggested by geological data (Calvet et al., 2020) including recent paleoelevation estimates (e.g.Huyghe et al., 2020). This uplift is consistent with the reconstruction of the evolution of river profiles suggesting a post-30 Ma incision of the whole of Iberia (Conway-Jones et al., 2019). In the absence of significant crustal shortening, the drivers of this post-10 Ma uplift are heavily debated. They include changes related to increasing temperature in the asthenopheric mantle, regional scale mantle flow, lithospheric folding, or magmatic addition in the lithospheric mantle (Cloetingh et al., 2002; Boschi et al., 2010; Faccenna et al., 2014; Conway-Jones et al., 2019; Huyghe et al., 2020).

Timing of Ebro Basin opening to the Mediterranean Sea

The timing and factors controlling the opening of the Ebro basin toward the Mediterranean Sea have long been debated, mainly due to the limited chronological constraints. Extension in the Catalan Coastal Range (CCR) triggered rift flank uplift, which, together with dry climatic conditions, favored the endorheic basin stage (García-Castellanos et al., 2003). This extension finally resulted in the Valencia Trough opening during Miocene (Fontboté et al., 1990; Roca et al., 1990, 1999; Sàbat et al., 1995; López-Blanco, 2002). Coney et al. (1996) proposed that the Valencia Trough opening was responsible for the capture of the Ebro paleo-lake. Seismic and sedimentary analysis of the Ebro delta in the Valencia Trough have shown that the first record of siliciclastic sediments coming from the Iberian margin is Serravalian-Tortonian in age (Castellon Group, Evans and Arche, 2002). García-Castellanos et al. (2003) then proposed that the age and amount of sediments preserved in the Valencia Trough are consistent with an onset of basin opening between 13 and 8.5 Ma, further facilitated by wetter climatic conditions that increased base level in the paleo-lake. The timing of basin opening has been refined to 12–7.5 Ma by combining the flexural isostatic compensation of the eroded volume with available constraints on sediment age (García-Castellanos and Larrasoaña, 2015). This timing of basin incision is coherent with low-temperature thermochronological studies reporting AHe cooling ages from the SPZ (Fillon et al., 2013) and AHe ages from the western Axial Zone in the Pyrenees (Bosch et al., 2016; Fillon et al., 2020) documenting post-orogenic exhumation between 9 and 11 Ma. In its lower reaches, the Ebro seems to have crossed the CCR before the Messinian Salinity Crisis (MSC), as evidenced by clasts sourced from the Pyrenees (Arasa-Tuliesa and Cabrera, 2018). The MSC did not lead to the entrenchment of the Ebro at the bottom of a deep canyon, unlike what has been observed for the Nile (Chumakov, 1967) and the Rhone (Loget et al., 2005; Mocochain et al., 2009). This issue was resolved by detailed offshore seismic constraints that document canyons buried offshore in the Ebro delta filled by pre-Messinian sediments (Urgeles et al., 2011). Based on the above arguments, we consider the Ebro Basin opening toward the Mediterranean Sea occurred between 7.5 and 13 Ma (see also Calvet et al., 2020 for a synthesis).

Beyond the issues about the age of incision of the Ebro Basin and whether its capture is related to sediment overfilling or headward retreat incision (García-Castellanos et al., 2003), little is known about the evolution of a river network after the onset of the capture of Ebro basin.

Ebro Basin paleoclimates

At the beginning of our study period, during the Miocene, the climate of the Ebro basin was subtropical and dry (López Martínez et al., 1987; Barrón et al., 2010), favoring the development of endorheic lakes (García-Castellanos, 2006). During the Middle-Late Miocene and Early Pliocene, northern Iberia recorded alternations of cold-wet and hot-dry periods (i.e., Bessais and Cravatte, 1988; Barrón et al., 2010). Afterwards, the Panama Isthmus closure at ∼ 3 Ma induced the establishment of the present-day oceanic circulation and global climate (cf.Maier-Reimer et al., 1990; Haug and Tiedemann, 1998; Molnar, 2008). This led, during the Quaternary, to more humid and colder climate conditions characterized by dry glacial periods and humid interglacials, especially for the Mediterranean area (Suc and Popescu, 2005; Jiménez-Moreno et al., 2013). Finally, during the Lower-Middle Pleistocene Transition (MPT, ∼ 1.3–0.7 Ma), the dominant orbital periodicity changed from 41 ka to 100 ka (Tziperman and Gildor, 2003; Lisiecki and Raymo, 2005; Siddall et al., 2010).

Review of Ebro Basin fluvial terraces, with emphasis on the topmost terraces

In the Ebro Basin, the flights of fluvial terraces along the main rivers and tributaries are widespread and can be organized to first-order into two groups (e.g.Lucha et al., 2012). The most elevated terrace levels (Fig. 1), at > 100–150 m above the current river level, consist of widespread and high-elevation pediment surfaces covered by up to ∼ 20 m of rounded pebbles of varying lithology. These upper levels are often located at the summit of hillcrests and seem locally disconnected from the present drainage network (e.g.Bomer, 1979). On the opposite, below ∼ 120 m, most of the terrace levels have a limited extent within modern valleys and along river courses. This morphological break highlights a transition between wandering fluvial systems with efficient lateral planation and beveling capabilities (e.g.Cook et al., 2014; Bufe et al., 2016, 2017) to entrenched and more stable river courses (e.g.Bomer, 1979). The term “mobile” refers to lateral river mobility, enhancing lateral erosion, whose greatest effect is valley widening up to the beveling of large platforms. For example, the Barbastro-Balaguer anticline has been completely flatten during the setup of these high levels (Lucha et al., 2012) before being subsequently entrenched by the Cinca, Noguera-Ribargorzana, and Segre rivers. The contrasting behaviour of the rivers after this transition is well illustrated by the terraces sequences of the Cinca (Lewis et al., 2017) or Alcanadre rivers (Calle et al., 2013). It may correspond to a transition from a depositional stage under sub-arid conditions, toward a more erosive stage characterized by a more focused incision. This could correspond to the climatic change recorded during the Lower-Middle Pleistocene transition (e.g.Stange et al., 2016). The formation of terraces younger than the Middle Pleistocene transition appears to be related to global climatic cycles (e.g.Macklin et al., 2002; Santisteban and Schulte, 2007; Arboleya et al., 2008), though other parameters such as uplift may be also important in defining the number of terraces and their altitudinal distribution (Santisteban and Schulte, 2007).

Due to the marked interest in the relationship between the incision and glacial cycles, most previous studies in the Ebro and Duero basins have focused on the lowest terrace system, close to the actual course of the rivers (e.g.Lewis et al., 2009; Calle et al., 2013; Stange et al., 2013a, 2013b). In addition, the lower terraces are much easier to date with classical methods like OSL, TL, ESR, 10Be cosmonuclide, and even 14C, usually much more efficient for ages under 300 ka (e.g., Fuller et al., 1998; Lewis et al., 2009, 2017; Calle et al., 2013; Stange et al., 2013a; Duval et al., 2015, 2017; Sancho et al., 2016, 2018; Delmas et al., 2018; Soria Jáuregui et al., 2019).

The highest preserved terraces levels in the Ebro and Duero basins, supposed to be the oldest evidence of fluvial incision, have often been considered as regionally coeval based on their elevation only (Bomer, 1979; Julián Andrés, 1996; Santisteban and Schulte, 2007; Lucha et al., 2012). This relative approach is important because dates alone sometimes do not allow deciphering the various terraces. Moreover, these dates are difficult to obtain. They are mainly derived from ESR, magnetostratigraphy, and more recently cosmogenic pairs of nuclides (Fig. 1). In the following paragraphs, we review the labels and dates of these highest levels and the chronological transition toward the lower levels, i.e. when the river began to entrench. We adopt a roughly East-West progression, first for the Pyrenean (northern) reaches and then for the southern reaches.

Along the Segre river flowing from the Pyrenees, 8 levels of terraces have been documented (Stange et al., 2013b), as well as in its major tributary, the Noguera Ribagorçana. The three top-most terrace levels are covering the interfluve to the west in between the Noguera Ribagorçana and the Cinca Rivers (the “highest” levels described before; Bomer, 1979; Stange et al., 2013a) but have not been documented along the Segre river. These three top-most terraces are characterized by important cementation (Bomer, 1979; Delmas, 2019). They lie 110–160 m above the present-day Noguera Ribagorçana river course (Bomer, 1979; Lucha et al., 2012; Stange et al., 2013a). The uppermost level has been dated to be ∼ 200 ka BP by inversion of 10Be depth profiles by Stange et al. (2013a) (Fig. 1), but this date is questioned regarding the inversion technique adopted (Delmas, 2019) and the correlation with the Ebro high fluvial terraces (see below, Bomer, 1979; Lucha et al., 2012). In the upper Noguera Pallaresa, a tributary of the Segre, cave deposits were recently dated using the cosmogenic burial method to ∼ 1.1 Ma (26Al,10Be pair in Llenes cave, ∼ 100 m above the current river course, Fig. 1) (Genti, 2015).

Along the Cinca river, 10–11 terrace levels are documented (Bomer, 1979; Peña-Monné and Sancho, 1988; Santisteban and Schulte, 2007; Lewis et al., 2009, 2017). The two topmost levels are equivalent to the three topmost levels mentioned to the west of the Noguera Ribagorçana valley. They are labelled Qt1 and Qt2 by Santisteban and Schulte (2007) and are ∼ 180 m and ∼ 130 m, respectively, above the present-day river course. The terrace level immediately below (Qt3 in Santisteban and Schulte, 2007), shows a normal magnetic polarity whereas Qt1 and Qt2 polarities are reverse. In addition, based on the degree of soil weathering, Lewis et al. (2017) proposed an age of 401 ± 117 ka for Qt3. Qt5 is dated to 178 ± 21 ka (OSL, Lewis et al., 2009) and Lewis et al. (2017) estimate likely ages of 999–1070 ka and 780–999 ka for Qt1 and Qt2, respectively, on the basis of paleomagnetic data. ESR analyses were also performed on these terraces, yielding ages of ∼ 1.3 and ∼ 0.8 Ma for Qt1 and Qt2, respectively (Duval et al., 2015; Sancho et al., 2016). Remnants of Qt1 and Qt2 are located on the interfluve between the Cinca and the Alcanadre rivers (Fig. 1).

The Alcanadre river is the next main river to the west and is a tributary of the Cinca river. 9 terrace levels have been mapped along the Alcanadre river. The two topmost levels, Qt1 and Qt2 were presented in the preceding paragraph. They lie at 160–200 m and 100–190 m, respectively, above the present-day Alcanadre river, and are characterized by inverse magnetic polarity (Calle et al., 2013; also dated with ESR as discussed in the following Duval et al., 2017) (Fig. 1).

The next river to the west is the Gállego river. 10 to 12 terraces are described along its course (Benito et al., 2000; Santisteban and Schulte, 2007 and references; Lewis et al., 2009; Benito et al., 2010). The 5–6 uppermost levels are older than the Brunhes-Matuyama transition. The Brunhes-Matuyama transition is at ∼ 150 m above the current Gállego river (Santisteban and Schulte, 2007). To the south of the city of Zaragoza, the Rio Huerva displays a similar, well-developed terrace sequence up to 150 m above the river course, but not dated (Guerrero et al., 2008). The highest terraces described along the Ebro river course are those found close to Zaragoza (Benito et al., 2000) and 30 km downstream (T1 and T2 ∼ 200 m above the current Ebro course, Guerrero et al., 2013) (Fig. 1). Locally, the terrace sequence is disturbed by karstic processes; this is particularly true in the Zaragoza area (Mensua and Ibañez, 1977; Benito et al., 1998, 2000, Guerrero et al., 2008, 2013). These remarks probably explain that the important Gállego set of terraces do not appear in large-scale correlation (e.g., Santisteban and Schulte, 2007).

The last important tributary sourced from the Pyrenees is the Aragón river. It has been mostly studied in its upper reaches in the mountain, for the glacially-derived terraces (García-Ruiz et al., 2013). Along its course within the Ebro Basin, Bomer (1979) indicates 5 fluvial terraces. The topmost one lies 120–140 m above the present-day river and builds a large plateau testifying for a mobile river when it formed. This terrace is likely equivalent to the large set of terraces found at the top of the Ribagorzana and Cinca valleys (Bomer, 1979). In its lower part, it is highly deformed by diapiric movements (Casas et al., 1994). To the north of the Aragón river, along its Cidacos tributary, an extensive terrace sequence is preserved (Fig. 1). The difference in elevation of the Cidacos terraces appears to increase towards the Ebro Basin.

In the upper reaches of the Ebro river, high terrace levels surrounded by glacis are observed (Soria Jáuregui et al., 2019). In the Duero basin, close to the drainage divide with the Ebro Basin (Vacherat et al., 2018), the highest Arlanzon river terraces T3, T4, and T5, were dated, using ESR, between 1.1 Ma, and ∼ 0.6 Ma (Moreno et al., 2012) (Fig. 1).

To the south, numerous rivers with documented terraces are draining toward the Ebro, like the Najerilla river, which has a well-preserved terrace record composed of 10 levels, but with sparse chronological information (Julián Andrés, 1996) (Fig. 1). The only one southern tributary with chronological data is the Guadalope river whose terrace record comprises 13 levels with some IRSL chronological up to 250 ka BP (Fuller et al., 1996, 1998). Along the Guadalope river, the uppermost terraces, lying ∼ 56 and 81 m above the current river stream, could be equivalent to the topmost terraces along the Cinca and Gállego valleys (Santisteban and Schulte, 2007).

Methods

Towards a unified model of Ebro terrace succession

There has been a little attempt, so far, to provide a unified model of terrace succession over the Ebro drainage basin. Tentative correlations have been proposed however in some parts of the basin, mainly along the Segre and its tributaries: Alcanadre, Cinca, Noguera Rigaborçana (Bomer, 1979; Alberto et al., 1983; Stange et al., 2013b). Rare correlations between terraces along distant valleys were proposed, like the correlation between the Gállego and Alcanadre river terrace sequences (Mensua and Ibañez, 1977; Lewis et al., 2017) or the one between main river terraces on both sides of the Pyrenees (Delmas et al., 2018; improved from Cordier et al., 2017). No attempt was made to correlate with fluvial terraces further west, probably because of the distance and because the area of Zaragoza and the Gállego valley is disturbed by local deformations produced by the mobilization and dissolution of evaporite (Benito et al., 1998, 2000, Guerrero et al., 2008, 2013).

We have correlated and unified the terrace labels over the entire basin at a scale of 1:50 000 (workflow on Fig. 2). We have first considered individual geographical areas (usually a valley) before to evaluate correlations between the successive locations. At the valley scale, we used georeferenced 1:50 000 geomorphological maps (Mapa Geomorfológico de España a escala 1:50 000 published by IGME; UTM 30 or 31 reference system; info.igme.es/cartografiadigital/tematica/Geomorfologico50.aspx?language=es). These maps contain a detailed mapping of the terraces and a nomenclature in alphabetical order starting with ’a’ for the oldest terrace. The terrace layouts (polygons) were taken from the digital geological map (MAGNA 50 – Mapa Geológico de España a escala 1:50 000 (2ª Serie) with shapefiles, IGME; http://info.igme.es/cartografiadigital/geologica/Magna50.aspx) when the contour on the geomorphological map corresponded to the boundaries on the digital geological map (Fig. 2). When the difference between the two was judged to be too large, we manually retraced the contour of the terrace as drawn on the geomorphological map. In the rare case of disagreement from one map to another for a single valley, we chose to join or instead distinguish between entities. These geomorphological maps cover only ∼ 20% of the basin, but most of the main valleys east of Zaragoza. In the areas not covered, we relied on the geological map and manual contouring of flat areas from the shaded DEM (MDT25 dataset from IGN) and derived slope map (Fig. 2).

The mapping of terraces has been improved by comparison with published works on the terraces of the studied valleys (e.g.Benito et al., 2000; Stange et al., 2013b; Lewis et al., 2017), as well as broader syntheses from which we georeferenced the main synthesis maps (e.g.Alberto et al., 1983; Benito et al., 2000; Calle et al., 2013) (Fig. 2).

This work of detailed mapping focused only on the highest terraces. We did not include any terrace in the immediate vicinity of the active stream, either geographically or in terms of elevation (< 50 m difference in elevation). Based on geological and geomorphological maps and published works (e.g.Lewis et al., 2017), we have verified that a priori all the terraces older than 300–400 ka were mapped. This work resulted in more than 500 polygons, with ∼ 250 of them being associated with a terrace level identified in the published geomorphological map or identified in a publication.

The second part of the work was to correlate the sequences of terraces. For the multiple tributaries of the Segre to the east, this work has been facilitated by the number of previously published works and existing geomorphological maps, although without any real homogenization for terrace labels (see Stange et al., 2013a). Similarly, the sequence of the Gállego is well defined (Benito et al., 1998, 2000). We correlate these two areas based on the mapping by Alberto et al. (1983) of the high terraces and the one by Mensua and Ibañez (1977) of the terraces and glacis in a large area around Zaragoza. For the Aragón and Ebro valleys, the sequences of terraces are far from the Segre/Gállego area. Despite this difficulty, we established a correlation which is first based on the lower sequence, in particular the level E, a large ubiquitous level that can be found in every valley of the Ebro Basin. This correlation is then refined on the basis of differences in altitude between the terraces (e.g.Julián Andrés and Chueca Cía, 1998), degrees of preservation, and an attempt of correlation by Bomer (1979), who correlated the topmost surfaces of the Segre and the Aragón valleys based on their geomorphological position and the presence of highly cemented conglomerates.

At the end of this work, certain ambiguities remained. Ultimately, to resolve these ambiguities, we have used the ages and magnetic polarities published in the literature as well as the new age constraints presented below (Fig. 2). The latter adjustments are proved to be minor, indicating our homogenisation procedure is quite robust. A few remaining mismatches reveal a potential error limited to a difference of one level. Differences in terrace level recognition also occur between our map and one of the published ones, or between two different publications (for example between the maps from Mensua and Ibañez, 1977 and from Alberto et al., 1983). Following our terrace mapping workflow, we labeled each terrace level using capital letters, similar to the convention used on geomorphological maps and deliberately different from the labels used in publications, which are different for each river studied. The shapefile is available at https://doi.org/10.5281/zenodo.4721256.

Cosmogenic nuclide dating methods used (depth profiles and P-PINI)

In-situ cosmogenic nuclides are produced within the Earth’s surface material through nuclear reactions (mostly spallation but also muonic capture) promoted by cosmic rays. Spallation processes occur within a couple of meters below the Earth’s surface (e.g., Lal, 1991; Gosse and Phillips 2001). The concentration of in situ-produced (i.e., produced within the crystal lattice) reflects the time the mineral spent in the production zone (c.a. from surface to few upper meters below surface). They have been extensively used for dating alluvial surfaces, assuming that both the initial mineral concentration at deposition time (thereafter called inheritance) and the denudation rate of the alluvial terrace are negligible (e.g., Bierman et al., 1995; Brown et al., 1998; Regard et al., 2005). For old terraces (> 300 ka), surface reworking may be quite significant. To overcome this problem and avoid inheritance effects, it is advisable to analyze the depth profiles of cosmogenic nuclide concentration (e.g., Anderson et al., 1996; Repka et al., 1997; Ritz et al., 2003; Hidy et al., 2010), although this cannot allow to date terraces much older than 500 ka. For longer timescales, cosmogenic nuclides can be used by pairs. The analysis is very different in this case: the two cosmogenic isotopes are produced in a similar fashion in the same mineral during his life nearby the Earth’s surface, allowing the knowledge of their initial concentration ratio. Then, the mineral is buried far from the surface, where cosmogenic nuclide production is negligible. Thus, the concentration ratio evolves following the various disintegration rates of the two cosmogenic nuclides. This method is called burial method (Balco and Shuster, 2009). It may be modified to take into account the production after the mineral has been buried, which has been popularized under the name of isochron method (Balco and Rovey, 2008). This method allows dating alluvial terraces up to a couple of million years (e.g., Balco and Shuster, 2009; Jungers and Heimsath, 2016).

In this paper, we report new chronological constraints obtained with various cosmogenic nuclide techniques: exposure depth profiles, simple burial and P-PINI, which is an adaptation of the isochron burial method.

Basic systematics

The basics of cosmogenic nuclide systematics can be found in Gosse and Phillips (2001) and Dunai (2010). We sum up in the following the most important information. Cosmogenic nuclides are produced at Earth’s surface by secondary cosmic rays (neutrons muons). Cosmic rays are quickly absorbed at depth, so that the production rate exponentially decreases against the mass of overlying material: P(z) = P0eρz/Λi, where P0 is the surface production rate (atom.g−1.yr−1), ρ the density (g.cm−3), and Λ (g.cm−2) corresponds to the attenuation length of cosmic-ray particles. Production accounts for three types of secondary cosmic rays: neutrons, slow muons and fast muons (thereafter corresponding to i = 1, 2 and 3, respectively). At an alluvial surface, the concentration C in cosmogenic nuclide depends on the age of the alluvial surface (t), its erosion rate (ε) and the initial cosmogenic nuclide content, often called inheritance (H), and the decay constant (λ):
C(t)=Heλt+iP0eρϵt/Λiλ+ρϵΛi(1e(λ+ρϵΛi)t).
graphic

If t and ε are unique for a single terrace, H varies for each sample. To determine the most probable (t, ε) pair, depth profiles are often used. Numerical processing of concentration vs. depth data is not straightforward and numerical procedures have been proposed, in particular, the one we used, developed by Hidy et al. (2010). They explore (t, ε) pairs and select the most probable ones by minimizing the difference between the observed and modelled concentrations.

Input parameters in the program provided by Hidy et al. (2010) are a density of 2 to 2.4 g.cm−3, a neutron attenuation length of 160 ± 10 g.cm−2 and a minimum number of profiles tested of 100 000. Other site-dependent input parameters were fixed after a few iterations to encompass all the range of the likely values. They are the explored age range, the maximum erosion rate, an erosion threshold. We also set the maximum inheritance to a value close to the minimum measured concentration.

Different cosmogenic nuclides can be used. For quartz-rich lithologies, the most convenient one is 10Be. Interestingly, 26Al can be also used. 10Be and 26Al are produced proportionally; the radioactive decay of 26Al being is approximately twice as fast as the 10Be (half-lives of 0.705 and 1.387 Ma, respectively). Consequently, the concentrations of surface samples in 10Be and 26Al are linked (their ratio R is ∼ 6.6 if their source area denudation rate is not too slow). If buried, it is possible to record the burial time by recording the evolution of the ratio in the two nuclides, following their disintegration. The burial must be at high depth, far below the surface, where there is no more cosmogenic nuclide produced. This holds true for depth > 10 m. On the contrary, it is possible that production during burial is not negligible in alluvial sequences, thus the burial duration is a minimum estimate. It is possible to correct for this post-burial production by comparing the concentration in cosmogenic nuclide of various rocks sampled at a similar depth, following a method called isochron (e.g., Balco and Rovey, 2008; Jungers and Heimsath, 2016). P-PINI (Particle Pathway Inversion of Nuclide Inventories) is an adapted isochron method that has been recently setup to account for complex pre-burial histories (Knudsen et al., 2020). The P-PINI method applies a Monte Carlo source-to-sink model approach, in which a large range of possible exposure and erosion histories are simulated for both the source and sink areas. For the source area, this involves highly complex exposure histories and abrupt erosion events associated with the presence of ice covers. The source model thus includes three episodes of ice cover with durations ranging from 2–100 ka for the two former ones and 2–10 ka for the latter one. These ice-covered periods are characterized by a complete halt in the production of cosmogenic nuclides, and a stochastic erosion event with an erosion depth randomly selected from the range 1–10 m. Periods with no ice cover are characterized by the production of cosmogenic nuclides with a range of production rates that corresponds to the elevation range of the estimated source area. The ice-free periods are accompanied by uniform and continuous erosion with rates randomly selected from the range 1–1000 m/Ma. Each simulation of the exposure and erosion history in the source area produces a (26Al,10Be) nuclide pair that is used as input for the sink model. This model simulates a large range of possible burial ages, assuming that the sediments were not covered by ice after deposition and that they experienced continuous erosion with rates in the range 1–500 m/Ma throughout the burial period. The production rate in the sink model corresponds to the location and elevation of the study site after correction for the sampling depth. Altogether, the combined P-PINI model approach includes 1 × 107 simulations for each study site, and a method known as rejection sampling is used to accept simulations that yield 10Be and 26Al concentrations that agree with measured concentrations. P-PINI includes the production of 10Be and 26Al associated with spallation, negative muon capture, and fast muons. A full in-depth description of the P-PINI model approach can be found in Knudsen et al. (2020).

Sampling sites

Segre terraces

The highest Segre river terrace is characterized by a large, ∼ 30 km long pedimentation surface pediment striking NNW – SSE made of three distinct sub-levels, lying on the interfluve between the Noguera Ribagorzana and Cinca rivers (Bomer, 1979; Lucha et al., 2012; Stange et al., 2013a, 2013b). This set of surfaces is limited to the north, by the South Pyrenean Frontal Thrust (SPFT), close to the upstream valley of the Noguera Ribagorzana river. The highest Segre river terrace levels the Barbastro-Balaguer anticline, which outlines the SPFT (Lucha et al., 2012). The organization of terrace remnants suggests an eastward shift of the Noguera Ribagorzana during deposition.

From these three alluvial surfaces, the level Qt0a (B in our terminology, see thereafter), where TSEG2 has been collected (0.573069°E, 41.70244°N, 357 m; Fig. 3), is the highest, ∼ 160 m above the present-day Noguera Ribagorzana river. It consists in ∼ 10 m thick of gravels deposited on top of Oligocene strata, with intercalation of sandy channels, topped by a very thin soil layer (< 10 cm). Although the top surface is flat, it is slightly dipping towards the valley indicating lateral erosion on the edges. We collected 6 samples along a depth profile from 15 to 380 cm from the surface (TSEG2, Fig. 3). We choose a sampling site with a relatively flat top to minimize post-deposition erosion, on an ancient quarry site west of the locality of Rossello.

We also sampled (TSEG1) the lowest of these uppermost levels, Terrace Qt1 (D in our terminology). This terrace is situated at ∼ 115 m above the present-day active channel of the Noguera Ribagorzana river. It consists of polygenic imbricated conglomerates overlying Oligocene strata. The top of this section (from 0 to ∼ 30 cm) is cemented by a calcrete as mentioned by Bomer (1979). TSEG1 (0.551336°E, 41.83316°N, 389 m) consists in 7 samples along a depth profile from 15 to 300 cm from the surface (Fig. 3). The sampling site, south to the Alfarras locality, is a small ancient quarry with no evidence of important surface erosion.

Alcanadre terrace

Sample TALC1 (−0.034394°E, 41.908655°N, 486 m) comes from an upper terrace level that forms the divide between the Alcanadre and Cinca rivers (Terrace Qt1 of Calle et al., 2013; Duval et al., 2015). This large terrace (∼ 15 km long; B terrace level, see the following) is gently sloping southward (0.5%). Its elevation above the modern Alcanadre river decreases from 200 m to the North to 160 m to the south. At the sampling place, to the SW of Peralta de Alcofea, the top of the terrace is at 176 m above the current river. Samples come from a thick (> 10 m) succession of conglomerates with well-rounded gravels and pebbles in a sandy matrix and occasional large cross-beddings and sand lenses (Fig. 3). The upper part shows the development of thick soil (∼ 3 m) with massive calcrete. Sample TALC1 consists of quartz pebbles and sand taken at depth of 8.50 m below the top of the terrace. Paleomagnetic data of Calle et al. (2013) indicate normal magnetic polarity near the base of the succession (at depth of ∼ 10 m) but normal magnetic polarity above. ESR data of Duval et al. (2015) give ages of 1606 ± 187 and 1269 ± 157 ka for a sample at 14 m depth (sample ALC1202) and 1504 ± 229 and 1282 ± 139 ka for a sample from the upper part (depth of 4 m) of the sequence (sample ALC1201). Duval et al. (2015) consider that the most reliable age of this terrace is 1276 ± 104 ka.

Aragón terrace

The highest Aragón river terrace forms a ∼ 25 km-long pedimentation surface (T1), slightly dipping (slope ca. 0.3%) toward the Ebro river to the south. This terrace T1 surface splits into four remnant pieces from north to south and is situated ∼ 140 m above the Aragón river. Lower T2 to T7 terrace levels are deeply entrenched in the Aragón valley contrary to T1.

Two samples have been collected on this surface (A terrace level, see the following). The first sampling site (TARA0, −1.60158°E, 42.19839°N, 395 m) is located on the southernmost remnant and corresponds to a quarry, to the NNW of the Arguedas locality. It is characterized by>20 m of gravels with intercalations of sandy channels. The top (> 1 m) appears highly cemented by calcrete. We collected five samples from rounded quartz pebbles to sand for burial dating, at ∼ 20 m depth (TARA0, Fig. 4). The second sampling site (TARA2, −1.419875°E, 42.38064°N, 466 m) is located on the northernmost remnant, on an ancient quarry, NE of the locality of Carcastillo. The top is highly cemented by calcrete, up to ∼ 1 m depth. We collected 5 pebbles of sandstone and quartz at ∼ 9 m depth for burial dating (TARA2, Fig. 4).

Ebro terrace

The highest preserved terrace record of the Ebro river (T1), located between Logroño to the west and Calahorra to the east, is largely affected by erosion and shows important alteration and soil development. Terrace T1 is situated at ∼ 200 m above the active Ebro channel (likely B terrace level, see the following). We collected 5 pebbles of quartz (TEBR2, −2.207328°E, 42.40442°N, 504 m, Fig. 4) at ∼ 10 m depth for burial dating.

Najerilla terrace

The Najerilla river is a tributary of the Ebro whose headwaters are located in the Iberian Chain. The highest terrace record in the vicinity of the Najerilla river corresponds to a staircase of four distinct terraces, characterized by a slightly NW dipping (slope ca. ∼ 0.6%). The lowest level (T4) is the most extensive (∼ 10 km large), the highest (T1) only show a limited remnant (∼ 700 m large).

Terrace T1 is situated at ∼ 245 m above the Najerilla active channel (a priori corresponding to A or B terrace level, see the following). Field observations indicate important surface erosion on the top of the terrace and soil development. We collected 5 samples of quartz and sandstone pebbles at ∼ 7.5 m depth for burial dating (TOJA1, −2.775386°E, 42.461922°N, 701 m).

Laboratory analyses and age calculations

All the samples were processed at CEREGE following standard procedures (cf.Braucher et al., 2015). They were measured at ASTER AMS (Tab. 1; Arnold et al., 2010).

Results

Mapping and homogenization of terrace sequences resulted in the morphological maps presented in Figures 5 and 6. The mapping is comprehensive for the highest levels from A to D. The terraces corresponding to level E have been mapped where they represent an important reference in the establishment of the sequence (e.g. along the Gállego or Aragón rivers, Figs. 5 and 6).

We recognize the highest, oldest terrace level, labelled “A”, in the upper Ebro and Aragón valleys (Figs. 5 and 6). This level A constitutes probably the uppermost terrace remnants, sparsely remaining along the Ebro valley, lying more than 200 m above the current river (Fig. 7a). This level is deformed near the confluence between the Aragón and Ebro rivers, suggesting it is much older than the subsequent ones. “B”, “C”, and “D” levels have been defined following the three levels identified on the interfluve between the Noguera Rigaborcana and the Cinca rivers (Alberto et al., 1983; Stange et al., 2013a; Lewis et al., 2017) (Figs. 5 and 6). Remnants of last level “D” are not parallel to the current river network and lie generally ∼ 100 m over the current rivers (Fig. 7). “E” is a usually extensive level, usually lying 40–60 m above current river courses (Figs. 5 and 6); it corresponds to Qt4 of Lewis et al. (2017).

Dating results

For TSEG1, eight samples have been measured (Tab. 1 and Fig. 8). Two 10Be Concentrations cannot be measured (Tab. 1). The 26Al concentration of the point at 50 cm depth is probably erroneous (it is 3 to 6 times less than the neighboring samples, Tab. 1) and has been discarded. 10Be and 26Al profiles lead to best age estimates of 693 ka (most probable age range 523–1095 ka), and 422 ka (most probable age range 400–971 ka), respectively (Tab. 2). The profiles indicate that erosion rates are most probably less than 1 m/Ma, and certainly less than 2.1 m/Ma (Tab. 2). The overlap of 10Be and 26Al results defines a common range of 523–971 ka for TSEG1 (Figs. 5 and 6).

Six TSEG2 samples have been measured (Fig. 8). 10Be and 26Al profiles lead to best age estimates of 532 ka (most probable age range 489–1961 ka), and, 435 ka (most probable age range 374–1104 ka), respectively (Tab. 2). Erosion rates are most probably less than 1.8 m/Ma (Tab. 2). The overlap of 10Be and 26Al results defines a common range of 489–1104 ka for TSEG2 (Figs. 5 and 6).

For TARA0, P-PINI indicates a best terrace age of 4731 ± 1447 ka (Tab. 3). For TARA2, P-PINI computed 2215 ± 778 ka (Tab. 3). For TEBR2, TOJA1, and TALC1, P-PINI indicates ages of, respectively, 1034 ± 275 ka, 957 ± 161 ka, and 5016 ± 1702 ka (Tab. 3 and Figs. 5 and 6). The large uncertainty for TALC1 and TARA2 comes from a low (insufficient) number of samples and that the concentrations are too similar to be well constrained by the P-PINI or Isochron method (Balco and Rovey, 2008).

Unified terrace chronology

Previously, we have presented our relative terrace chronology, leading to terrace levels labelled in alphabetical order starting with A. We now analyze the absolute chronology.

Previous works mainly focused in the eastern part of the basin. In the Segre catchment, study sites from B and C levels provided inverse magnetic polarity while the magnetic polarity is normal for sites from D (Fig. 9) (Calle et al., 2013; Lewis et al., 2017). ESR dates proposed are 1276 ± 104 ka and 817 ± 68 ka (Duval et al., 2015), for terraces mapped as B and C, respectively (Fig. 9). Our P-PINI burial age for B at TALC1 is 5016 ± 1702 ka appears erroneous maybe because it is poorly constrained, probably because of too little data. Our dates from cosmogenic depth profiles are 489–1104 ka for B (TSEG2) and 523–971 ka for D (TSEG1, Figs. 6 and 9). Besides, Lewis et al. (2017) proposed an age of 401 ± 117 ka for the next terrace level (E, Fig. 9).

Towards the west, around the Aragón valley, even if it is tempting to correlate the topmost terrace of Bomer (1979) with the B level, the topmost terrace is likely to be much older, as emphasized by the old cosmogenic nuclide ages (> 2 Ma) found for TARA0, TARA2 upstream from the confluence from the Aragón and Ebro rivers (Figs. 5, 6 and 9). This antiquity is further confirmed by the deformation of this terrace level, whereas the lower levels, including the well-developed E level in the area, are not deformed (see Casas et al., 1994 and the corresponding 1:50 000 geological map).

The A level was not previously dated. It corresponds to two sites from the present study (TARA0 and TARA2). The P-PINI dates for these two sites are quite different but both agree on a fairly old age. None of the dates are very reliable due to the low number of data. Since we cannot choose between the two, we adopt the mean age of ∼ 2.8 ± 0.7 Ma (Tab. 4). We note, however, that this age represents well the samples TARA0-B, TARA0-D, TARA2-C and TARA2-D (Supplementary Material, Figs. S10 and S13). The B level was deposited with inverse magnetic polarity (Calle et al., 2013; Lewis et al., 2017), and dated at ∼ 1.3 Ma (Duval et al., 2017). TSEG2, sampled from this level, yields an age of 489–1104 ka, by a cosmogenic profile method (Fig. 9). Also, TOJA1 (960 ± 160 ka), in the upper reaches of the Ebro main valley, could correspond to this level. We had some doubt regarding which terrace level, A or B, TEBR2 belongs to. Its P-PINI age of 1.03 ± 0.27 Ma clearly falls in the range of the level B (Tab. 4). Consequently, and considering B is older than C (see below), we regard an age of 1.15 ± 0.15 Ma to be likely for this level B (Tab. 4). The age of C, the next level, is constrained by an inverse magnetic polarity (Lewis et al., 2017), and one ESR estimate of ∼ 0.8 Ma (Duval et al., 2017). Its likely age is 850 ± 70 ka after our analysis (Fig. 9 and Tab. 4). The terrace level D mostly displays a normal magnetic polarity (Calle et al., 2013; Lewis et al., 2017). We estimate an age of 523 to 971 ka (TSEG1) from our cosmogenic depth profile that is much older than the estimate by Stange et al. (2013a). Together with the C and E terrace dates, which must be older and younger, respectively, we propose an age of ∼ 650 ± 130 ka for this D level (Fig. 9 and Tab. 4). Finally, terrace E has been dated to be 401 ± 117 ka (Lewis et al., 2017) (Fig. 9). We summarize the absolute age constraints on Table 4: we estimate ages for terraces A, B, C, D, E to be 2.8 ± 0.7 Ma, 1.15 ± 0.15 Ma, 850 ± 70 ka, 650 ± 130 ka and 400 ± 120 ka, respectively (Fig. 9 and Tab. 4).

Discussion

The youngest sediments of the endorheic stage preserved in the Ebro Basin are late Miocene (see Section 2.2) lacustrine and fluvial deposits (Pérez-Rivarés et al., 2002; García-Castellanos et al., 2003; Pérez-Rivarés et al., 2004; García-Castellanos, 2006; Larrasoaña et al., 2006; Vázquez-Urbez et al., 2013; Calvet et al., 2020). They are observed 500 m above the Ebro river course, and have possibly been covered by younger sediments now eroded (García-Castellanos and Larrasoaña, 2015). The subsequent erosion of the Ebro basin was limited and affected mostly the Miocene deposits that are currently preserved in the form of alluvial deposits but also pediments. Plio-Pleistocene alluvial remnants are mostly terraces found along (and parallel to) the main rivers in valleys within the basin. Some of these remnants are lying on interfluve crests between the main valleys. Among them, good examples are represented by terraces between the Noguera Ribagorçana and the Cinca valleys (levels B, C, and D). These uppermost terraces argue for a mobile fluvial network while the lowermost ones reflect river entrenchment without significant lateral mobility. This change in fluvial behavior seems to have occurred at the time of deposition of terrace level D, shortly after the Brunhes-Matuyama transition at 780 ka (Fig. 9). This agrees with the change in fluvial style at the Mid-Pleistocene Transition proposed by Sancho et al. (2016).

Terrace chronology

The proposed chronology for terraces A, B, C, D, E is 2.8 ± 0.7 Ma, 1.15 ± 0.15 Ma, 850 ± 70 ka, 650 ± 130 ka and 400 ± 120 ka, respectively (Fig. 9). These dates are the result of a compilation of ages obtained with various methods (Tab. 4): ESR, paleomagnetism, cosmogenic nuclide depth profiles, and cosmogenic nuclide burial ages from the isotope pair (10Be, 26Al) (in caves and in terraces through P-PINI inversion). The number and quality of data for each level of terrace are variable. For example, paleomagnetism only indicates if the terrace level is younger or older than the Brunhes-Matuyama transition (780 ka). Also, after our own experience, cosmogenic nuclide depth profiles are not accurate when the ages exceed 500 ka, and can be significantly altered by terrace post-abandonment surface evolution. For example, the ages computed from depth profiles published in Stange et al. (2013a) are probably underestimated (see also Delmas et al., 2018; Delmas, 2019). Similarly, our own depth profiles, which do not give an unequivocal result, probably tend to underestimate the ages, possibly because of a recent acceleration of erosion. The 2.8 ± 0.7 Ma for the level A is not very good because the two estimates on the uppermost Aragón terrace largely differ (TARA0 at 4730 ± 1450 ka and TARA2: 2210 ± 780 ka). The level B benefits from a large number of dating attempts. Its P-PINI cosmogenic burial dates agree with ESR (except for TALC1) and an inverse magnetic polarity. The age suggested by our depth profile data (TSEG-2) is a bit younger than the average age of ∼ 1150 ka calculated for level B but is older than the age proposed by Stange et al. (2013a). Level C delivered only two but overlapping pieces of information, making it reliable: an inverse magnetic polarity (Lewis et al., 2009) and an ESR age of 817 ± 68 ka (Duval et al., 2017). On level D, we produced a depth-profile cosmogenic age, which, despite the aforementioned caveats, and with the constraint of normal magnetic polarity, allows to restrict the range of possible ages to 650 ± 130 ka.

In the nearby Arlanzon valley in the upper part of the Duero basin (location shown in Fig. 1), Moreno et al. (2012) used the ESR technique to date the uppermost terraces: 1140 ± 130 ka for the T3 terrace, 780 ± 120 ka and 939 ± 100 ka for the inverse polarity T4 terrace, 700 ± 100 ka, 600 ± 100 ka, and 700 ± 70 ka for T5, and 400 ± 90 ka, and 370 ± 70 ka for T8. These ages agree well with our own estimates: the Arlanzon terraces T8, T5, T4, and T3 correlate with the Ebro Basin terraces E, D, C and B respectively (cf.Fig. 1). Similarly, on the northern side of the Pyrenees, the uppermost Têt river terrace (T5) has provided two ESR dates at 1099 ± 179 ka and 1133 ± 159 ka (Delmas et al., 2018), again contemporaneous to B. Along the Têt river, age of T4 is still debated while the T3 terrace level below has been dated with ESR technique to be ∼ 400 ka (440 ± 39 ka, and 374 ± 47 ka)(Delmas et al., 2018). The Têt terraces T5 and T3 should then correspond to our B and E levels, respectively, as already noted by Delmas et al. (2018). In the Pyrenees, some karst dates fall within the age range we propose in the Ebro basin. The Llenes cave, located ∼ 100 m above the present course of the river, and dated by 26Al/10Be burial dating to 1120 ± 160 ka (Genti, 2015), could certainly be correlated with level B (see Calvet et al., 2015 and refs included for the correlation between cave systems and terraces). More doubtful is Cave 5 of Agosto which has a similar age (1240 ± 220 ka Genti, 2015), but at a much higher elevation (800 m above the present river). Finally, it appears that the chronology of the Ebro basin shares lots of similarities with adjacent areas. We conclude that the distribution of terrace ages is a regional signal, not specific to the Ebro Basin.

We observed that rivers begin to entrench around terrace levels B to D: it marks the transition from large terraces now lying on interfluves to sequences of young terraces developing parallel to the current river courses. On the basis of physical models (Bufe et al., 2016; Baynes et al., 2020) as well as natural systems (Bufe et al., 2017), it turns out that large terraces are produced by high channel mobility, which in turn depends on water and sediment fluxes, themselves a function of climate. As already proposed by Sancho et al. (2016), it is tempting to relate this entrenchment to climate as it has occurred during or shortly after the Middle Pleistocene Transition (MPT, 1.3 to 0.7 Ma ago). Indeed, the MPT is due to dominant orbital periodicity change from 41 ka to 100 ka (Tziperman and Gildor, 2003; Lisiecki and Raymo, 2005; Siddall et al., 2010). This led to an increase of climatic contrasts, triggering intense glacier fluctuations in the surrounding mountain ranges and controlling sediment discharge and river incisions, similar to what we observe in the Ebro Basin (Stange et al., 2013a, 2013b; Lewis et al., 2017). Our observations further suggest that before the MPT, the upper Pliocene/lower Pleistocene climate favored a mobile river network. This idea is, on the one hand, challenged by Santisteban and Schulte (2007) who claimed for a diachronous incision beginning along the Iberian rivers. On the other hand, the MPT increase in incision rates is attested all over Europe (Gibbard and Lewin, 2009). We infer that the river network has been fixed when the climate underwent larger and longer glacial/interglacial oscillations.

Incision: rates and significance

The incision appears to have accelerated for the last million years. In the center of the basin, around Zaragoza, the youngest endorheic sediments are currently lying at an elevation of ∼ 700 m (Pérez-Rivarés et al., 2002, 2004; Vázquez-Urbez et al., 2013), ∼ 500 m above the nearby Ebro river. A and B terraces lying 200–220 m and 140 m above Ebro river, respectively. The total amount of incision is ∼ 300 m between the Ebro Basin opening (13–7.5 Ma) and the topmost terraces A. This is a minimum, given the possibility these strata have been buried under a couple of hundreds of meters of younger sediments before being exhumed (García-Castellanos and Larrasoaña, 2015). Thus, in the Zaragoza area, the incision was ∼ 40 m/Ma before 2.8 Ma then ∼ 40 m/Ma for the period 2.8–1.15 Ma and ∼ 120 m/Ma (x3) since then (Fig. 10). Along the Aragón valley, A and B are 130 m and 90 m above the current river, respectively: the incision rate varies from ∼ 25 m/Ma for the period 2.8–1.15 Ma to ∼ 80 m/Ma (x3, for uncertainties see Table S1, Supplementary Material) for the last 1.15 Ma (Fig. 10). We tested the robustness of these observations taking into account various ages and terrace elevation datasets presented in Figure S1 (Supplementary Material). It indicates that most of the Ebro basin rivers entrenched at a rate of 110 to 220 m/Ma. As proposed by Sancho et al. (2016), such a rate cannot have lasted for many millions of years, in agreement with a recent acceleration of incision. This analysis also confirms that incision acceleration may have occurred around 1 Ma ago, as shown by the two areas displaying old terraces: the lower Gállego and the Aragón valleys (Fig. S1, Supplementary Material). In the Duero basin, Rodríguez-Rodríguez et al. (2020) documented an increase in catchment-wide denudation rate by factor of two or three (x2 to x3) since 1 to 2 Ma. However, the authors relate this acceleration to expansion of the Duero drainage network driven by long-term base-level lowering and autogenic processes rather than to a climatic signal.

We note that the incision calculated in the Zaragoza area of 500 m is possibly underestimated if the latest pre-opening sediments have been buried before exhumation (García-Castellanos and Larrasoaña, 2015). Thus the pre-2.8 Ma period incision is (slightly?) lower than in the southern flank of the Pyrenean Axial Zone for the last 9 Ma deduced from low-temperature thermochronology (700–1600 m, equivalent to an exhumation rate of 130–500 m/Ma, Fillon et al., 2013). This could be partly explained by the larger amount of sediments at the foot of the axial zone, where paleovalleys were filled by up to 2000 m of conglomeratic sediments now eroded (e.g.Vincent, 2001; Babault et al., 2005a; Fillon et al., 2013). It also indicates that Neogene erosion has been stronger at the foot of the Axial Zone than in the centre of the basin, which is consistent with a significant uplift of the Pyrenees relatively to the basin since 10 Ma (e.g.Conway-Jones et al., 2019; Calvet et al., 2020).

Most upper terraces developed at the outlet of Pyrenean streams when they enter into the Ebro sedimentary basin. Extensive terrace sequences developed similarly along rivers flowing from the Iberian Chain to the south. Rivers flowing mainly within the Ebro Basin, in contrast, show less extensive terrace development, as along the Arba river. In the NW of the basin, a comparison between the Cidacos and Arga rivers is enlightening: despite of small length, the Cidacos river left large terraces, while the longer nearby Arga river shows sparse terraces with little lateral extent. Maybe, as the Arga is larger, it eroded more the ancient terraces and the Cidacos terrace remnants could be conceived more like pediments? Nevertheless, it suggests that large terrace sequences are found if (i) they form large/thick surfaces, and (ii) later river mobility was not sufficient to remove them.

Terrace distribution also indicates that older terraces developed preferentially out of the basin center and from east to west (Santisteban and Schulte, 2007). In the eastern part of the basin, the absence of the topmost level A, together with the location of the lower levels on top of the interfluves suggests that a large reshaping occurred at the Plio-Pleistocene transition. This reshaping may have produced a pediment landscape such as in the current Arba catchment (Centre-north of the basin). There, large pediments are found, in a depressed position relative to the Aragón terraces to the north-west. This is consistent with observations from other rivers of Europe for which slight incision documented during the Latest Pliocene/Early Pleistocene led to the development of wide and shallow valleys (Gibbard and Lewin, 2009). The east-west difference is also marked by a decrease of incision to the west, which could be explained by upstream incision wave migration, as recently documented in the Duero (Rodríguez-Rodríguez et al., 2020). This incision wave is currently affecting the drainage divide between the Duero and Ebro catchments (Vacherat et al., 2018; Struth et al., 2019). Hence, the large and old terraces found upstream could correspond to the smooth landscape characterizing the area before Late Pliocene/Pleistocene incision. The geomorphic processes active at this time may have been diffuse, explaining the large scatter between the ages found (cf. burial dates of TARA0, TARA2). This upstream incision wave migration is potentially enhanced by an erosional isostatic rebound more pronounced in the center and eastern part of the basin than in the western part (García-Castellanos and Larrasoaña, 2015). This isostatically induced rebound is due to basin unfilling but can also be caused by moderate to important denudation of the Pyrenees and Central Iberian Massif in Plio-Pleistocene times (García-Castellanos and Larrasoaña, 2015; Stange et al., 2016). In the center-east part of the basin, this uplift does not necessarily differ from the center of the basin to the Pyrenees foothills, explaining why terrace levels are parallel such as for the Segre and its tributaries (Fig. 7b). This observation seems to be valid also for the Gállego river (Fig. 7c), despite the karstic deformation occurring in its lower reaches (Benito et al., 1998, 2000, Guerrero et al., 2008, 2013). To the west, the Aragón river terraces show a clear decreasing incision trend downstream (Fig. 7e), while the trend is opposite along the upper Ebro River (e.g. the Oja/Najerilla system near the Ebro-Duero drainage divide). This trend is however not well resolved and would need additional works if we are to discuss its correlation with the post-Miocene tilt shown by lacustrine strata covering the Ebro and Duero basins (García-Castellanos and Larrasoaña, 2015). At a shorter time scale, incision rates significantly increased as already noted by Sancho et al. (2016). Post MPT increased erosion is likely controlled by erosion waves characterized by knickpoints migrating upstream, possibly originated at the coast as a consequence of eustatic fluctuations (e.g.Antoine et al., 2000; Crosby and Whipple, 2006; Bridgland and Westaway, 2008; Loget and Van Den Driessche, 2009). While migrating upstream, the amplitude of erosion waves could be attenuated with higher incision rates downstream as recently documented on the Duero (Rodríguez-Rodríguez et al., 2020). Hence, the more pronounced incision in the east of the basin may be driven by a position closer to the Mediterranean Sea, possibly enhanced by a Quaternary flexural erosional uplift of the entire basin up to the Catalan Coastal Range (Lewis et al., 2000; Stange et al., 2016).

Conclusion

A unified fluvial terrace chronology and mapping have been produced for the Ebro Basin. It is based on a review of ages from the literature and new cosmogenic nuclide ages obtained with the profile and P-PINI methods. From the topmost terrace downward, the terraces are labelled in alphabetical order. Terraces A, B, C, D, and E are dated to 2.8 ± 0.7 Ma, 1.15 ± 0.15 Ma, 850 ± 70 ka, 650 ± 130 ka and 400 ± 120 ka, respectively, in agreement with other dated terraces in the Iberian Peninsula and around the Pyrenees. The oldest ones (A, B, and C) are extensive, testifying for wandering fluvial systems while from D to present, the network has been stabilized and is currently entrenched in 100 to 200 m-deep valleys. We therefore hypothesize that a transition occurred from mobile to localized fluvial network. We propose this transition is likely related to the Middle Pleistocene Transition (MPT, between 0.7 and 1.3 Ma), when long-period/high-intensity climate fluctuations were established in Europe. We estimate that between 2.8–1.15 Ma and present, the incision rates have tripled, from 25–40 m/Ma to 80–120 m/Ma. After the MPT, the terrace disposition also shows more important incision and basin unfilling in the east possibly enhanced by flexural rebound.

Acknowledgements

V. Zavala and J. Rat participated to the field work. We thank M. Delmas, R. Braucher, and S. Carretier for fruitful discussions. This work benefited from support and funding by OROGEN TOTAL-BRGM-CNRS project. We thank two reviewers, M. Calvet and D. García-Castellanos, and Editor A. Teixell for their careful reading and constructive comments.

Cite this article as: Regard V, Vacherat A, Bonnet S, Mouthereau F, Nørgaard J, Knudsen MF. 2021. Late Pliocene-Pleistocene incision in the Ebro Basin (North Spain), BSGF - Earth Sciences Bulletin 192: 30.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.