First Chronological Constraints for the High Terraces of the Upper Ebro Catchment

: The Cenozoic sedimentary basins in the Iberian Peninsula show a change from long-term basin inﬁll to incision, a transition that indicates a period of major drainage reorganization that culminated in the throughﬂow of the networks to the Atlantic and Mediterranean oceans. Both the cause of the transition from aggradation to degradation and the linkages to tectonic, climatic, and geomorphic events hinge on the chronology of the ﬂuvial network incision and excavation of the basin’s sedimentary ﬁlls. In this paper, we describe the ﬁrst chronologic data on the highest ﬂuvial terraces of the upper area of the Ebro River, one of the largest ﬂuvial systems in the Iberian Peninsula, to determine the onset of incision and excavation in the basin. For this purpose, we combine electron spin resonance (ESR) and paleomagnetism methods to date strath terraces found at 140, 90, and 85 m above the current river level. Our results show ages of ca. 1.2 and 1.5 Ma for the uppermost river terraces in the upper Ebro catchment, constraining the minimum age of the entrenchment of the upper Ebro River.


Introduction
The sedimentary records of the Tertiary basins in the Iberian Peninsula show a change from long-term basin infill to degradation, a transition that marks a period of major drainage reorganization culminating in the throughflow of the main river networks to the Atlantic Ocean and Mediterranean Sea (e.g., [1][2][3][4][5][6]).The timing of the change from endorheic to exorheic conditions and whether this was linked to tectonic, climatic, and geomorphic events hinge on the chronology of the fluvial network incision and excavation of the basin's fill.Although river incisions of most of the basins in Iberia appear to have initiated in the Quaternary, the precise timing of their beginning is not fully understood (e.g., [4,7]), except for certain rivers, such as the Tagus [8] Fluvial terrace staircases provide direct records of both base-level and river incisions, and are known to be a result of the combination of climatic changes and tectonic uplift (e.g., [9][10][11][12]).Therefore, providing direct chronometric ages for such fluvial terraces will elucidate the switch from an endorheic to exorheic context in the Iberian fluvial basins.
One of the Iberian basins that experienced such evolution is the Ebro Basin.This basin became the South Pyrenean Foreland Basin, bounded by the Catalan Coastal Ranges to the east the Iberian Range to the south, and the Pyrenees to the north.The formation and development of the present-day Ebro River fluvial network was mainly controlled by the tectonic evolution of the Pyrenean orogenic belt during Tertiary times.The basin started as an endorheic system, isolated from the ocean for several million years, until it opened up to the Mediterranean (e.g., [13]).The timing and process(es) leading to the major drainage reorganization that eventually led to the basin's opening are currently not well defined (e.g., [14][15][16][17][18][19][20]).
Magnetostratigraphic data suggest the Middle Miocene period as the age of the youngest preserved endorheic sediments in the Ebro Basin (from 12-11 Ma ago) ( [18,[21][22][23]).A conglomeratic unit near the Iberian Range (labelled unit T8 by [24]), which post-dates the Middle Miocene lacustrine carbonates, may reflect exorheic drainage conditions [22], but whether they alternatively correspond to an internal alluvial system cannot be ruled out.Additionally, clastic deposits, or a mantled pediment, of Pyrenean provenance that truncate the culminating lacustrine limestone in the center of the Ebro catchment [25] may also (although not necessarily) represent early exorheic drainage, although reliable numerical chronologies are still missing.Overall, there is no solid evidence of post-Middle Miocene endorheic sedimentation in the Ebro foreland basin, a fact that suggests the lower age limit for the beginning of the exorheism.Many fewer studies have been devoted to the timing and formation of the fluvial terraces formed by the progressive downcutting of the fluvial network that followed the aperture of the basin to the Mediterranean Sea (e.g., [26][27][28][29][30][31][32][33][34][35][36][37]).Even less data exist for the terraces of the trunk river [38][39][40][41].
The longest terrace sequence of the Ebro River is formed by 11 terrace levels [42,43], whereas the most complete terrace systems in the Ebro catchment are made up of 12 levels associated with the Gállego [27] and Huerva [44] rivers.So far, the oldest numerical age constraint available for a fluvial terrace in the Ebro catchment comes from the Alcanadre River, which provides an age of ca.1.3 Ma [45].In that river, terraces T1 and T2, at ca. +180 and +120 m, respectively, above the current river level, produce dominant reverse magnetization directions, indicating a Lower Pleistocene age.ESR dating of optically bleached quartz grains from those two terraces provided ages of 1276 ± 104 ka and 817 ± 68 ka for T1 and T2, respectively [46], which is consistent with the palaeomagnetic results.In the Gállego River, the highest terrace is at +150-175 m. 26 analyzed the terrace sequence and suggested that the Matuyama-Brunhes reversal (0.78 Ma) occurred between terraces at +85-90 m and +60-70 m near Gurrea.Finally, [38] presented further paleomagnetic data for a terrace level located at +210 m above the Ebro River, supporting a Lower Pleistocene age (pre-Matuyama-Brunhes boundary); however, it is worth mentioning that these last two studies did not yield any new numerical age results.The corollary is that a time gap exists between the youngest clear and dated evidence for endorheism in the Ebro foreland basin sector, of Middle Miocene age, and the oldest fluvial terraces, dated at about 1.3 Ma, at about 180 m above the present-day river.
In the mid-Ebro River, a number of localities show gravelly deposits that are higher than the studied terraces, and therefore presumably older.For example, the locality known as Loma Negra (Zaragoza) shows a glacis level at about 400 m above the present-day river [25], but has an unknown chronology.Additionally, a number of fluvial terraces east of Alfaro are found at about 200 above the Ebro River [47], which presumably may be older than the targeted deposits.
Data are scarce in the upper Ebro catchment.According to the existing geological maps, the youngest evidence for endorheism in the upper sector is given by the lacustrine limestones of Upper Miocene [48].However, no ages are yet available for either the oldest fluvial terraces or the marginal alluvial fans in that area.The main goal of this paper is therefore to provide new data on the timing of formation of the highest river terraces preserved and generated by the Ebro River in the upper catchment [43].This will allow further constraint of the timing of the highest preserved fluvial aggradation phases across the catchment following the findings of the studies mentioned above.Electron spin resonance (ESR) and paleomagnetism, when possible, were used to determine the chronology of the two fluvial terraces in the upper Ebro River.

Regional Setting
The Ebro River, located in northeastern Spain, is one of the longest rivers and drains the second largest fluvial catchment in the Iberian Peninsula.The Ebro River finds its source in the southern foothills of the Basque-Cantabrian Mountains and, running eastward, flows for over 970 km into the Mediterranean Sea, forming the Ebro delta.The formation and development of the present-day Ebro River fluvial network was mainly controlled by the tectonic evolution of the surrounding Pyrenean orogenic belt during Cenozoic times.The study area is located in the upper Ebro catchment, structurally separated from the Ebro foreland basin by the Sierra Cantabria-Montes Obarenes thrust, the continuation of the Southern Pyrenean thrust (Figure 1) [49].Here, the morphostructure corresponds to a folded ejective relief, made up of a succession of wide synclines and narrow anticlines (Figure 1).The synclines constitute concordant negative reliefs (e.g., Villarcayo and Miranda-Treviño basins) [50], and the anticlines form positive reliefs (e.g., Sierra Cantabria-Montes Obarenes, Valderejo-Sobrón anticline).The river terraces that we studied are located specifically in the Miranda-Treviño basin sector (MT Basin), a piggyback basin that developed at the back of the Sierra Cantabria-Obarenes Mountains thrust (Figure 1).In this area, the trace of the river is largely controlled by the local geologic structure, following synclines and forming canyons across anticlines.The mountainous edges are mainly made of Cretaceous marine deposits (limestones, marls, and sandstones), whereas the MT Basin is filled by Paleogene and Neogene continental deposits (conglomerates, sandstones, and lacustrine limestones).The southern limit of the MT Basin is the Sierra de Cantabria-Montes Obarenes thrust, which places the Mesozoic cover over the Ebro foreland basin toward the south.The trace of this thrust fault and the associated mountain range delineates the southern margin of the MT Basin, and hence of the upper Ebro catchment.
The Ebro River flows through a canyon incised in Mesozoic rocks just north of Haro (Figure 1), leaving the MT Basin and entering the Ebro foreland basin.On the hanging wall block of the thrust where the MT Basin is found, an endorheic drainage system developed during at least the Late Miocene, as suggested by a unit of lacustrine limestones [48].A conglomeratic unit with limestone and quartz cobbles, unconformably overlying Oligocene sediments west of Miranda de Ebro, have been assigned to the Pliocene [48], although their age and significance have not yet been determined.The fact that they are tilted reveals some tectonic activity, which may correspond to the latest stages of the Sierra Cantabria-Obarenes Mountains thrust.In any event, sometime after that, but definitely after the Late Miocene, the Ebro River developed a sequence of 11 strath terrace levels over the Miocene sedimentary bedrock in the Miranda sector of the catchment [41,51,52].Their distribution is rather heterogeneous; they are more abundant toward the Miranda basin, and the lowermost terraces are generally better preserved (Figure 1).The majority of river terraces are mainly composed of well-sorted, well-rounded limestone pebbles and cobbles, often imbricated, and embedded in a sandy matrix [41].The river terraces that we studied are located in the upper Ebro catchment [41,43], in the upper area of the catchment (northern Spain), specifically in the Miranda basin (Figures 1 and 2).

Material and Methods
The terrace sequence of the upper Ebro catchment was redefined and correlated using aerial photographs (scale 1:18,000), orthophotos (scales 1:50,000, 1:25,000, 1:10,000, and 1:5000), and geological maps (IGME, scale 1:50,000).This information was implemented by field data.Elevation data were obtained from 5 m cell size LiDAR-derived DEMs published by the Instituto Geográfico Nacional (IGN).Fluvial facies were described in suitable outcrops and further sampled for ESR dating and paleomagnetic studies.

Sampling
We sampled the highest identified fluvial terraces, T1 and T2, located at 140, 90, and 85 m above the present-day river and near the town of Miranda de Ebro (Figure 2, Table 1, Figures S1-S3).A total of 5 sediment samples were collected for ESR dating purposes (see Supplementary Material SM1).Sampling was performed by hammering PVC tubes into layers of fine to medium sand.Two samples (MIR1601 and MIR1602) were taken from T1 deposits (+140 m.a.p.r.l.) at the so-called Monte Miranda-Yarritu locality.Three samples were collected from the T2 deposits at two localities, Monte Miranda-Arinorsa (+90 m; samples MIR1603 and MIR1604) and Aeródromo Miranda (+85 m; MIR1605).In situ measurements of the gamma dose rate were performed at each sampling spot.Additional sediment samples were collected for laboratory analyses (radioelement and water content).With regard to samples for paleomagnetism, oriented specimens were taken at each locality sampled for ESR (see Figures S1-S3).We used standard cylindrical plastic containers (~10 cc) pushed onto the sediments where they were cohesive enough.Otherwise, unconsolidated sediments were hardened with a solution (1:1) of sodium silicate before or just after orienting and removed from the outcrop to obtain ~10 cc cubes.All samples were oriented in situ with a standard compass and clinometer.Specimens were stored in a mu-metal cylinder in the laboratory to prevent further viscous magnetization.Sample preparation and ESR dose evaluation were performed at CENIEH (Burgos, Spain), using the same analytical procedure and experimental conditions as the previous work of [45]) on the Alcanadre deposits.Full details are provided in Supplementary Material.Quartz grains were dated via the multiple aliquots additive dose (MAAD) method and following the multiple center (MC) approach initially defined by [53].The ESR signals of the aluminum (Al) and titanium (Ti) centers for each sample were acquired in separate spectra using specifically optimized parameters.Each aliquot (14 in total: 1 natural, 1 optically bleached, and 12 gamma-irradiated) was measured 3 times after ~120 • rotation in the cavity for Al and Ti signals in order to consider angular dependence of the signal due to sample heterogeneity.To evaluate the repeatability of the D E values, measurements were repeated two to three times over several days.
The ESR intensity of the Al signal was extracted from peak-to-peak amplitude measurements between the top of the first peak (g = 2.0185) and the bottom of the 16th peak (g = 1.9928) [54].The ESR intensity of the Ti centers was measured following options A, C, D, and E sensu [46].Final dose response curves (DRCs) were obtained by pooling the mean ESR intensities and associated errors (1 standard deviation) derived from each repeated measurement, using either an exponential + linear function (EXP + LIN) for the Al center or the so-called Ti-2 function, as in [46].A complete description of the analytical procedure can be found in Supplementary Material SM2.Fitting results are given in Tables 2 and 3.

Dose Rate Evaluation
The total dose rate was obtained from a combination of in situ and laboratory measurements.In situ external gamma dose rates were obtained with a NaI or LaBr probe connected to a Canberra InSpector-1000 multichannel analyzer and calculated by the threshold technique [55].In addition, the corresponding radioelement (U, Th, K) concentrations for each sample were determined by ICP-MS analysis of dry raw sediment following a four-acid digest preparation procedure (Table 4).External alpha and beta dose rate components were obtained using the dose rate conversion factors from [56] to derive the concentration values, while the gamma dose rate was taken from the in situ measurement.Dose rates were calculated assuming a mean grain size of 150 µm and a removed thickness by HF etching of 20 µm.The internal dose rate was assumed to be 30 ± 10 uGy/a, based on work by [57] and using an alpha efficiency of 0.07 ± 0.01 [58].Values were corrected with beta and alpha attenuation values for spherical grains [59,60] and water attenuation formulae from [61].Current water content was obtained in the laboratory by drying the sediment at 50 • C in an oven for 3 weeks.A water content of 15% (% wet weight) in the sediment was assumed for the age calculation.The cosmic dose rate was calculated using formulae from [62], with depth, altitude, and latitude corrections [63].

ESR Age Calculation
ESR age was calculated using non-commercial Scilab based software, with error calculations based on Monte Carlo simulation, with the following sources of uncertainty: concentration, depth, water content, gamma dose rate, beta dose attenuation, and D E value.ESR ages are given as 1σ (Table 5).Recently, [64] showed that this software typi-cally provides dose rate values that are very close (within a few %) to those obtained by DRAC [65].In accordance with the MC approach, ESR age results were calculated for the Al, Ti (option D), and Ti-H (option C) centers (Table 5).Indeed, recent works showed that option D can provide overall accurate results for Early to Late Middle Pleistocene fluvial deposits (e.g., [66][67][68][69][70][71]), while option C (Ti-H) is usually more appropriate over a younger time range, from the Late Middle to Late Pleistocene [58,72].However, because we are aware of the existing variability within the community regarding how the ESR intensity of Ti centers is measured [73], we also provide the D E estimates of the Ti center derived from options A and E (Supplementary Table S1).If necessary, the reader can derive the corresponding ESR age estimates by using the total dose rate values given in Table 5.  1 and 2) and gamma source dose rate (2.3%).

Paleomagnetism
All measurements were made at the Archeomagnetism Laboratory of the Geochronology Facilities at CENIEH (Burgos, Spain).Natural remanent magnetization (NRM) and progressive demagnetization were measured in a 3-axis SQUID magnetometer (2G Enterprises) housed in a ~9 m 3 shielded space of Helmholtz coils (residual field <3000 nT).The cryogenic magnetometer has a built-in system of three degaussing coils for alternating field demagnetization to 170 mT.Thermal demagnetization was performed on a TD48 furnace (ASC) placed in the Helmholtz cage.The direction of the remanence components was determined by principal component analysis line-fitted after visual inspection of the demagnetization using Zijderveld-type plots.Characteristic remanent magnetization (ChRM) directions were computed using standard software and plotted on a lower hemisphere using equal area stereographic projection (data summary in Table 6).

Geomorphological and Sedimentological Characteristics
Terrace T1 is a topographically dominant feature of the landscape, located over the town of Miranda de Ebro (Figure 1).It is the only remnant of the culminating fluvial phase in the entire upper Ebro catchment.T1 elevation fluctuates between 600 and 594 m above sea level (+140 m above the Ebro River).The studied outcrop (Monte Miranda-Yarritu) shows a maximum visible thickness of approximately 32 m.The fluvial sediments linked to the outcrop are quite homogeneous, and are mainly made up of horizontally and crossstratified well-rounded pebbles and cobbles organized into metric planar beds, indicating deposition in channels and gravel bars [74].The coarse units are typically topped by massive sandy units, interpreted as the result of deposition during declining flooding events [75].The recurrent vertical stacking of coarse horizontally stratified units is overlaid by fine beds, pointing toward an extended period of sediment deposition.
The presence of terrace T2 is also minimal.This level is composed of only two remnants, elevated between +90 and +85 m above the Ebro River, and both located within the Miranda basin.The remnant closer to terrace T1, known as Monte Miranda-Arinorsa, is located at an elevation of 540 m above sea level and shows an elongated shape from south to north, reaching an approximate length of 600 m (Figures 1 and 2).The outcrop displays a maximum visible thickness of 37 m.From a sedimentological point of view, the fluvial sediments are composed of horizontally and cross-stratified well-rounded gravels organized into metric to decimetric beds bounded by concave-up and tabular erosive surfaces.These sediments are interpreted to indicate deposition in channels and gravel bars [74].Finer units are also present, especially close to the surface, pointing toward the final stages of fluvial aggradation during the formation of the terrace.The second remnant of T2 is a small terrace flight located at the end of a cul-de-sac valley, topographically bounded to the north and separated from the main valley by a Cretaceous limestone and dolomite crest (Figures 1 and 2).This terrace level occupies an area of 0.03 km 2 and is located +85 m over the Ebro River.The main outcrop has a visible thickness of 13 m, and is completely made of fine sediments.The bottom half of the outcrop is composed of trough cross-bedded sands representing deposition of sand dunes at the channel bottom.The top half of the sequence consists of various units of laminated fine-grained sands, silts, and mud that show desiccation cracks and bioturbation features, indicating deposition from suspension in an overbank environment [74].It is interpreted that this area functioned as a sediment trap during major flooding events and was subsequently abandoned by the progressive entrenchment of the fluvial valley.

ESR Data
ESR DRCs and fitting results obtained for the Al center are displayed in Figure 3 and Table 2, respectively.Bleaching coefficients vary within a very narrow range (49.1-50.6%),suggesting very similar bleaching conditions for the five samples.Repeatability of the ESR intensity is overall very good, with variation that does not exceed 2.5%.This results in excellent repeatability (variation of <6%) for four out of five samples, the only exception being MIR1605 (15.4%).Finally, three out of five samples show excellent goodness of fit, with adjusted r 2 value >0.99, while the other two (MIR1602 and MIR1603) show lower values, between 0.98 and 0.99 (Table 2).These proxies highlight the overall robustness of the ESR dataset.The calculated D E values cluster into two main groups, one with values in the range of 1600-1850 Gy (MIR1601, MIR1602, MIR1603) and one with lower values between 1100 and 1400 Gy (MIR1604, MIR1605).
ESR data obtained from options D and C (sensu [46]) are displayed in Table 4 (fitting results) and Figures 4 and 5. Fitting results derived from options A and E (sensu [46] are given in Supplementary Table S1 for comparison but will not be further discussed in this work.The variability of the ESR intensities obtained from the repeated measurements of the Ti center (option D) remains limited (<2%) for most of the samples, with the exception of MIR1601 (4.6%).This results in an overall slightly higher variability of the D E estimates compared with that achieved for the Al center (7.0 vs. 5.6% on average): repeatability is good (<10%) for four out of five samples, while MIR1605 shows a higher value (15.8%).Goodness of fit is excellent for most of the samples, with adjusted r 2 values of >0.99, the only exception being, again, MIR1605.In summary, all proxies indicate good robustness of the ESR data collected for all samples, except for MIR1605, which shows relatively worse D E repeatability and goodness of fit.D E values vary within a relatively narrow range (890-1060 Gy) for samples MIR1601 to MIR1604, while MIR1605 shows a lower value of about 740 Gy.
ESR DRCs and fitting results obtained for the Ti-H center (option C) are displayed in Figure 5 and Table 3, respectively.ESR data show relatively high variability of the ESR intensities over repeated measurements, resulting in regular to poor repeatability of the D E estimates (variation between 7.1 and 38.1%).Additionally, they provide overall poor goodness of fit, with only two samples showing an adjusted r 2 value above 0.95 (MIR1602 and MIR1604; Table 3).Together, these proxies suggest the limited reliability of the ESR data derived from the Ti-H center for these five samples, and especially for MIR1605, which displays very poor D E repeatability.Interestingly, D E estimates obtained for the five samples vary within a relatively narrow range between 520 and 630 Gy.
A comparison of the D E values derived from the Al and Ti centers for each sample is shown in Figure 6.All samples show the same D E pattern, with the Al center systematically providing the highest D E values and the Ti-H center (option C) the smallest ones.This may be interpreted at first glance as a typical bleaching pattern, and may simply reflect the differences in the bleaching kinetics of each center.Based on the principles of the MC approach, the D E pattern indicates that the Al signal was most likely not fully reset during sediment transport.Consequently, the Al center provides only a maximum possible estimate of the true burial dose.Although we might initially assume that the differences observed between Ti D E values derived from options C and D are also the result of an incomplete signal reset of the latter, the data should be interpreted with extreme caution.First, the robustness of the ESR data collected for the Ti-H option C was questioned earlier.Second, previous studies showed the limited reliability of the center to provide D E estimates above 300 Gy (see [73] and references therein).One of the very few exceptions was the Cuesta de la Bajada site, where the Ti-H signal provided an apparently accurate D E result of around 600 Gy.However, this is likely a unique situation resulting from very peculiar quartz samples showing exceptionally high ESR intensities of the Ti-H signal (see [76]).Third, the fact that all D E values were consistent within an error around 500-600 Gy might suggest signal saturation.Taken together, these issues suggest that the Ti-H center most likely underestimates the true burial dose and should therefore be considered as providing a minimum dose constraint in the present study.However, we cannot reasonably exclude that the relatively low environmental dose rate (between 450 and 1100 µGy/a) might instead make the Ti-H signal suitable for dose evaluation in the present study (see [73] and references therein).This question will be further discussed in the next section.Finally, Ti center option D yields D E estimates of between 740 and 1060 Gy, which is the typical dose range in which this signal produces accurate dose results (see [73]).Consequently, it most likely provides the closest estimate to the true burial dose rate.

Dose Rate Considerations
ICP analyses of the raw sediment samples show relatively low radioelement concentrations, with 0.61-0.75ppm of U, 1.67-3.11ppm of Th, and 0.24-0.87% of K (Table 4).These resulted in laboratory gamma dose rates ranging from about 190 to 373 µGy/a.In comparison, in situ gamma dose rates were between 22% lower (MIR1602) and 14% higher (MIR1604) than the laboratory-derived values.These values are nevertheless consistent within error for all but one sample (MIR1602), suggesting that the unmediated environment around the sampling spots is overall relatively homogeneous, which is consistent with field observations (Supplementary Figures S1-S3).Regarding MIR1602, the difference observed between laboratory and in situ values may be the result of either a relatively heterogeneous environment in the vicinity of the sampling spot, a slight disequilibrium in the U-238 decay chain, or possibly a combination of both.
Total dose rate values (Table 5) show an overall relatively low dose rate environment: they do not exceed 900 µGy/a, except for one sample (MIR1601, ~1100 µGy/a).This is, to our knowledge, very uncommon in the Iberian Peninsula.Our previous ESR dating works performed on different Spanish fluvial systems systematically showed significantly higher dose rate environments, e.g., 1850-3300 µGy/a for various terrace deposits associated with the Arlanzón River (Duero Basin) located a few tens of km farther South of Miranda de Ebro [77], 950-1750 µGy/a for the highest terraces of the Alcanadre River (Ebro basin) [5], 1150-2250 µGy/a for the T4 terrace of the Alfambra River valley (Alfambra-Teruel depression) [78], 3250-3950 µGy/a for T4 terrace deposits associated with the Miño River (Galicia) (69), 3250-7250 µGy/a for various deposits of the Jarama and Manzanares Rivers (Madrid) [79], and 1700-4500 µGy/a for different fluvial terrace systems in the eastern Cantabrian margin [80].

Paleomagnetism
Only the specimens from site MI04 showed stable behavior upon alternating field demagnetization (Figure 7).Most of the remaining specimens showed spurious or uninterpretable magnetization directions.In one sample there is a clear indication of a strong overprint, which cannot be successfully removed by alternating fields (specimen MI01-6, Figure 6).This specimen reveals decayed NRM up to 40 mT, above which level the magnetization directions are scattered.Nevertheless, the linear regression of the remanence vectors upon demagnetization bypass the origin in the Zijderveld diagram, revealing a great circle on the stereographic projection.This circle passes by a southward and negative direction, suggesting the presence of a higher coercivity component (most likely hematite) of reverse magnetization.Only locality (MI04 (Aeródromo de Miranda) produced reliable, interpretable directions (Figure 8).Progressive AF demagnetization often reveals two components of remanence.A soft, low-coercivity component is eliminated at a maximum field of 10 mT, after which the demagnetization trend decays linearly toward the origin of the coordinates in the Zijderveld diagram.Some specimens (e.g., MI04-9, Figure 7) show a slight curvy demagnetization trend, which may reflect gyroremanent magnetization produced during the AF demagnetization in the magnetometer.All but one characteristic remanent magnetization (ChRM) direction were southward and negative (Figure 8), consistent with reverse polarity.Mean inclination was 67 • (absolute value), indistinguishable from the expected geomagnetic axial dipole.

Age Determination
All samples show the same pattern for the ESR age results (Table 5): the Al center yields the oldest estimates by far, between 1.7 and 2.4 Ma, while the Ti center (option D) provides ages that are consistent with the second part of the Early Pleistocene (between 0.8 and 1.5 Ma), and finally the Ti-H center gives the youngest results, with late Early to Early-Middle Pleistocene age estimates (between 0.5 and 1.1 Ma).Following the principle of the MC approach, ESR age results derived from the Al center should be regarded as maximum possible age constraints.In other words, the true burial age should be similar to or younger than the Al ESR estimates.In contrast, the Ti center (option D) is expected to provide the closest estimate to the true burial age, given the age and dose range considered in the present study.Previous works focused on either the second part of the Early Pleistocene (e.g., [5,66,67]) or on D E estimates ranging between 700 and 1100 Gy [69], showing that Ti (option D) provides accurate ESR ages.In particular, the closest analogue to these deposits in terms of environmental dose rate (800-1350 µGy/a) and chronology (1.1-1.5 Ma) is probably the fluvial system of the Moulouya River [66,67].There, the late Early Pleistocene ESR ages derived from the Ti center (option D) were consistent with the independent age control provided by the luminescence and paleomagnetism methods.In contrast, for the same reasons, Ti-H is initially expected to yield underestimated age results, thus providing a minimum age constraint for the deposits.
The two samples collected from the T1 terrace deposits at +140 m, MIR1601 and MIR1602, yielded very scattered Ti ESR ages of 856 ± 83 and 1531 ± 155 ka, respectively.It is somewhat surprising that these samples, positioned in the highest (and therefore the oldest) terrace deposits, produced such an apparently young age, and we suspect that the age of MIR1601 is likely underestimated.Since the D E estimates are very similar for both samples, the cause of such discrepancy is likely to be related to the dose rate evaluation.MIR1601 returned the highest dose rate value of the dataset by far.Since the dosimetric data collected from the laboratory and in situ are consistent, we have no reason to question their reliability.The only possible explanation would be that the sediment experienced some remobilization of radioelements that may have locally impacted the radioactivity of the surrounding environment.Paleomagnetic data offer some perspective in this sense.The site under consideration, MI01, gives very poor results, as only three characteristic directions can be computed.All directions are downwards and north seeking, consistent with a normal polarity field.Both observations, poorly defined remanence directions and the appearance of north-seeking directions, might suggest a recent remagnetization of the sediments after the last major reversal (772 ka).To sum up, although we do not have empirical evidence to support the hypothesis of remobilization of radioelements, both ESR and paleomagnetic data seem to support the occurrence of post-depositional processes that would affect both the paleomagnetic signature and the underestimated ESR age.Consequently, the age of the formation of the T1 terrace is most likely better constrained by the ESR age of 1531 ± 155 ka obtained for MIR1602.
The two samples collected in the +90 m deposits, MIR1603 and MIR1604, yield very consistent Ti ESR ages of 1263 ± 123 and 1207 ± 106 ka, respectively.A weighted mean of 1231 ± 80 may be considered as the best age estimate for the T2 terrace.In contrast, sample MIR1605, collected from another outcrop, provides a much older Ti ESR age of 1579 ± 198 ka.Such ESR age is consistent with the reverse polarity revealed by the paleomagnetic analysis.The ESR result, nevertheless, seems to be strongly overestimated, because it is stratigraphically inconsistent with the other two obtained for T2.This may be partly due to the limited robustness of the ESR data collected for this sample, which creates significant uncertainty around the D E result (Table 3 and Figure 5).Moreover, sample MIR1605 displays the lowest dose rate of the dataset, and we cannot exclude that it may also play a role in the overestimation, since the ability of both Al and Ti (option D) signals to yield accurate estimates of small D E values at a maximum of a few hundreds of Gy may be reasonably questioned (e.g., [81]).
In contrast, the Ti-H center was found to be of particular interest in low-dose-rate environments (<1000 µGy/a) (see an overview in [73] and references therein) and in detecting small D E values (e.g., [71]).Therefore, we cannot reasonably rule out that it provides a reliable estimate of the burial dose for samples MIR1602 to MIR1605.However, the thermal lifetime of the signal is unknown, and if it has been demonstrated to be useful for Late Middle Pleistocene samples, its potential usefulness for older time ranges is simply unclear.In that regard, the Ti-H results obtained for the +140 m and +90 m deposits (samples 1601 to 1604) are stratigraphically inconsistent, and also in disagreement with the paleomagnetic results obtained for T2, although those are based on only two samples.Consequently, they should only be considered as minimum possible age estimates.In contrast, the Ti-H age of 1120 ± 163 ka (MIR1605) is consistent with the weighted mean age of 1231 ± 80 ka obtained for T2.This indirectly suggests that the Ti-H result obtained for this sample may be considered as a finite age rather than a minimum estimate.However, this age should be seen as purely indicative given the poor quality of the ESR data obtained for this sample.We consider that the formation age of terrace T2 remains better constrained by samples MIR1603 and MIR1604 and the resulting weighted mean age of 1231 ± 80 ka.In this regard, paleomagnetic data for such terrace at site MI04 (i.e., site MIR1605l; see Table 6) provides strong evidence for a reverse polarity chron, fully compatible with the proposed Matuyama age (>0.772Ma) for terrace T2.
Our new results, as well as published reports on terraces in the lower Ebro, suggest that we can reasonably discard the idea that this unexpectedly young chronology resulted from a systematic age underestimation by the ESR method due to signal saturation or limited lifetime.Indeed, these samples (i) are found in a relatively low-dose-rate environment, and (ii) do not show exceptionally high D E values (<1000 Gy for Ti, option D).Moreover, previous studies have shown that ESR could provide accurate age constraints for samples of at least 2.0 Ma (e.g., [70]).

Lower Pleistocene Terraces
Both T1 and T2 terraces date back to the mid-Early Pleistocene (1531 ± 155 and 1231 ± 80 ka, respectively).The only existing numerical data on terrace formation within the Ebro catchment comes from the Alcanadre tributary.There, Qt1 terrace (+200-160 m) was formed during the Middle to Early Pleistocene (1276 ± 104 ka) [82], coeval with the T2 of our study.As for other fluvial chronological data elsewhere in Iberia, 77 reported an ESR age of 1.14 ± 0.13 Ma for terrace T3 (+78-70 m) in the Arlanzón River, located in the Duero catchment.This age should, however, be regarded as a maximum age constraint for the deposits, as only the Al signal was measured.6 also suggested an Early Pleistocene age for the highest terrace levels in the Duero River based on 10 Be- 26 Al cosmogenic isotope depth profiles.Available datasets support the idea that both the Ebro and Duero catchments experienced fluvial aggradation phases during the Early Pleistocene.
The profound change in base level after the opening of the upper Ebro Basin toward the Ebro foreland basin through the MT Basin sometime after the Upper Miocene, played a major role in the partial excavation of the upper catchment basin and the beginning of the strath terrace system formation [41,51,52].In that sense, it has been suggested that the response time for a given fluvial system to recover equilibrium conditions after experiencing a major perturbation fluctuates between 1 and 5 Ma (e.g., [83,84]).The opening of a basin, a major perturbation, results in a new base level for the fluvial system that, in order to adjust to the new level, generates upstream-migrating knickpoints or a knickzone (e.g., [1,85,86]).As far as the Ebro catchment, it has been demonstrated that the Ebro River and some of its most important tributaries already adapted their courses to the Mediterranean base level and that the upper course of the river is also adapted to the local base level found in the foreland basin (e.g., [36,43,87].Therefore, there must be other factors to explain the development of fluvial terraces during the Quaternary in the Ebro fluvial system. Fluvial terrace staircases are known to be the result of alterations introduced within fluvial systems by the combination of climatic changes, tectonic uplift, and base-level changes (e.g., [10][11][12]88]).Previous studies in the upper Ebro catchment interpreted that fluvial activity was controlled by climatic factors which, in turn, were superimposed on a continuous valley entrenchment related to a long-term regional crustal uplift [41].The rate of this erosional dynamic can be calculated by using the height of the terraces and their ages and can be interpreted as a proxy for uplift rates.Estimated incision rates are 0.09 m/ka for T1 and 0.07 m/ka for T2.These values are very similar to those estimated for the Arlanzón River (0.06 m/ka) [77] in the Duero catchment and the Tagus River (0.1-0.07 m/ka; Cunha et al., 2008), but much lower than those found in the Alcanadre River, Ebro catchment (0.16-0.12 m/ka) [82].A number of uplift mechanisms affecting the Iberian Peninsula [89][90][91][92] may explain such contrast in uplift rates, but the scarcity of data does not allow definitive conclusions to be drawn about the uplift mechanisms affecting these Cenozoic basins.In the MT Basin in particular, estimated incision rates could be related to the crustal shortening happening since the Miocene in northern Iberia [93] and/or erosional isostatic adjustments taking place in the Ebro foreland basin [18].
Most studies link the formation of fluvial terraces to glacial-interglacial-driven changes in watershed hydrology and sediment flux (e.g., [94] and references therein).Under such premise, glacial-interglacial cycles are ultimately driving all or most of the incision/aggradation cycles.The sedimentological features of the resulting deposits in the MT Basin are very different, though; those deposited during cold stages are made up of coarse particles, whereas interglacial deposits seem to be composed of fine particles [41].In the upper Ebro catchment, terraces formed during those cold stages were correlated with higher water discharge and higher sediment delivery linked with glacial advances in the headwaters and periglacial conditions in the bounding reliefs.In turn, interglacial conditions favored an overbank deposition thanks to the lower water discharge and sediment delivery related to the expansion of vegetation cover and hillslope stabilization [41,43,[95][96][97][98]).In terms of paleoclimate, the timing of formation of T1 and T2 was around or at the onset of the so-called Early-Middle Pleistocene climate transition (EMPT).
Several records (e.g., [99]) reveal a significant change in climate dynamics starting at 1.4 Ma, when eccentricity-dominated (125 kyr) cycles become dominant [100].It has been reported that fluvial systems, at a global scale, experienced an increase in their incision rates, causing a deep entrenchment of their fluvial valleys around the EMPT [101].

Final Remarks and Conclusions
Our ESR study on quartz grains from fluvial sediments in the upper Ebro catchment within the MT Basin produces the first numerical age constraint for the highest terraces.Our best estimates for the ages of terraces T1 (+140 m) and T2 (+90-85 m) are 1531 ± 155 ka and 1231 ± 80 ka, respectively.Paleomagnetic results for terrace T2 furnish a clear reverse polarity in one site, fully consistent with the ESR results.ESR age underestimation in terrace T2 is observed at two localities, and is presumably related to post-depositional processes after the formation of the terrace.High porosity and permeability could definitely have favored such a process, which also could explain the presence of secondary magnetization directions in the deposit.Field observations show precipitates of Mn-oxides in some gravelly lenses within the terrace, which reveal that fluid circulation occurred at some point after gravel deposition.
At present, few chronometric age constraints exist for the oldest, highest Ebro fluvial terraces.Much farther east of the study area, in the central Ebro catchment, terraces T1 and T2 along the Alcanadre tributary, at ca. +180 and +120 m above the current river level, provided ESR ages of 1276 ± 104 ka and 817 ± 68 ka, respectively [5], so slightly younger yet not statistically different from our new ages for terrace T2 in the upper catchment.From a methodological point of view, these results are directly comparable with ours, as both are based on the exact same analytical procedure.The magnetic polarity in those terraces in the Alcanadre is reverse, showing altogether a Lower Pleistocene age, as in the Miranda area.The results described here, in the waterhead of the Ebro River, show ages almost coeval with the results in the mid-Ebro River, so fluvial terraces, from the geochronological point of view, could be coetaneous, as the ages are not statistically different.Be that as it may, the overall results indicate that the formation of the "modern" Ebro River and concomitant development of fluvial terraces were relatively recent.
Although it is beyond the scope of this paper, it worth pointing out that both the MT Basin and the Ebro foreland basin show evidence of internal drainage and lacustrine sedimentation until the Middle Miocene, so regardless of the timing of the present-day drainage network and the age of the connection between the upper and central Ebro River, the opening of the interior basin must have started at or after the mid-Miocene.In this regard, mantled piedmont truncating lacustrine carbonates in various locations in the Ebro Basin (e.g., [4,25]) indicate a major drainage change.Although such veneer of alluvium does not necessarily imply exorheic conditions, the deposits certainly predate the proper entrenchment into the bedrock.Therefore, the origin of the mantled piedmont could date back to the Pliocene or even Late Miocene.In this sense, our results are compatible with the conclusions by [14], in that the present-day Ebro River was not connected to the Mediterranean Sea before the Pliocene.The existence of a "proto Ebro", as envisioned by [102], back in the Pliocene or even earlier is still compatible with an Early Pleistocene catchment of the center and upper Ebro Basin, as shown in this study.The Castellón Group, a well-known progradation detrital unit (Middle-Late Miocene) in the Valencia Trough ahead of the present-day Ebro delta [103], certainly represents a significant input of sediments from the continent.Yet, the formation of such a clastic unit could have resulted from a catchment restricted to the Catalan coastal ranges [104], and therefore not necessarily be an indication of the opening (= connectivity) of the Ebro Basin to the Mediterranean (e.g., [14]).
This study presents chronological results obtained from the highest terraces of the Ebro River in its upper catchment within the MT Basin.The basin, according to the available geological data, was not connected to the main trunk of the Ebro River before the Late Miocene.The analysis of both studied terraces indicates that their formation was related to the onset of a period characterized by pronounced paleoclimatic changes, when climatic oscillations changed from 41 to quasi-100 ka cycles and the duration of interglacials experienced a relative decrease (the EMPT) [100].A general entrenchment of the river valleys took place around the onset of the EMPT.Incision rates ranging between 0.09 and 0.07 m/ka have been estimated, which are similar to values obtained in other Iberian fluvial systems but substantially lower than those calculated for some Pyrenean tributaries of the Ebro River.
At the scale of the Iberian Peninsula, although it is beyond the scope of this paper, it is worth noting that recent ESR data on fluvial terraces of the Tajo River also provided much older ages than anticipated (>1.8 Ma) [8].Overall, and besides the fact that Tajo and Ebro drain to the Atlantic and Mediterranean, respectively, it seems that, as in the Ebro case, the Atlantic draining networks might have acted as an exorheic system much earlier than anticipated [105][106][107][108].
Finally, our results confirm the utility of combining multiple, independent dating techniques to chronologically constrain fluvial deposits, in line with previous works by [66,68].As observed in the present study, this combination is essential to reliably date Lower Pleistocene sediment, as it is the only way to overcome some limitations that may exist with the use of a single method on a given sample.Future work will be focused on enlarging the number of samples, in particular from the higher terraces of the Ebro River, for ESR and paleomagnetic purposes in order to reinforce and complete the chronostratigraphic framework that is being established.

Figure 2 .
Figure 2. Map of terrace levels developed in Miranda Basin.Bottom left inset shows reconstructed river profile for each terrace level based on terrace surface points extracted from a 5 m resolution DEM (color range is the same as that used for terraces on the map; see legend).Below are field images of terrace deposits: A, Monte Miranda-Yarritu; B, Monte Miranda-Arinorsa; C Aeródromo.White capital letters refer to water laid deposits: channels (CH); gravel bars (GB); sandy bedforms (SB); floodplain fines (FF).

Figure 3 .
Figure 3. ESR dose response curves (DRCs) of the Al center obtained for five samples dated in present work.All ESR intensities derived from repeated ESR measurements were pooled for fitting.Vertical error bars represent experimental errors associated with ESR intensities, corresponding to 1 standard deviation.

Figure 4 .
Figure 4. ESR dose response curves (DRCs) of Ti centers (option D) obtained for five samples dated in present work.All ESR intensities derived from repeated ESR measurements were pooled for fitting.Vertical error bars represent experimental errors associated with ESR intensities, corresponding to 1 standard deviation.

Figure 5 .
Figure 5. ESR dose response curves (DRCs) of Ti-H center (option C) obtained for five samples dated in present work.All ESR intensities derived from repeated ESR measurements were pooled for fitting.Vertical error bars represent experimental errors associated with ESR intensities, corresponding to 1 standard deviation.

Figure 6 .
Figure 6.Multiple center approach.Comparison of D E results derived from measurements of Al and Ti centers.

Figure 7 .
Figure 7. Progressive demagnetization diagrams of representative samples.Each data point represents the NRM end vector for individual demagnetization steps projected onto the horizontal (solid symbols) and vertical (open symbols) plane on the Zijderveld diagram.All samples were demagnetized with alternating field (steps shown in milliTesla).Natural remanent magnetization (NRM) shows initial magnetic intensity previous to AF procedure.Sample MIR01-6 shows a linear trend that bypass coordinate origin during AF demagnetization.On stereographic projection, a great circle is observed, which suggests the presence of a higher coercivity reverse magnetization component unresolved with AF.Specimens MI04-8 and MI04-6 show progressive decay after 10 mT, where a viscous component is erased.

Figure 8 .
Figure 8. Lower hemisphere stereographic projection of paleomagnetic directions of site MI04.Stereonet on left shows characteristic remanent magnetization (ChRM) direction, high coercivity components.Stereonet on right shows viscous, low coercivity component.Open (closed) symbols represent vectors projected onto upper (lower) hemisphere.

Table 1 .
Sample location and details of studied fluvial terraces.

Table 2 .
ESR data derived from measurement of Al center.Bleaching coefficient is expressed as relative difference between ESR intensities of natural and bleached aliquots.Repeatability of ESR intensities is assessed through variability (1 relative standard deviation) of mean ESR intensities obtained after each day of measurement.Similarly, repeatability of D E values corresponds to variability (1 relative standard deviation) of D E values calculated for each day of measurement.

Table 3 .
ESR data derived from measurement of Ti centers.Repeatability of ESR intensities is assessed through variability (1 relative standard deviation) of mean ESR intensities obtained after each day of measurement.Similarly, repeatability of D E values corresponds to variability (1 relative standard deviation) of D E values calculated for each day of measurement.Option D (Mixture of Ti-Li and Ti-H) Option C (Pure Ti-H)

Table 4 .
[57]oelement concentrations measured from raw sediment by ICP-MS analysis.Laboratory gamma dose rate was calculated using dose rate conversion factors from[57], measured water content, and radioelement concentrations.In situ gamma dose rates are based on threshold approach.Probe employed indicated in parentheses.

Table 5 .
Detail of ESR age estimates and dose rate components.Errors are 1 sigma.Note that D E errors are a combination of fitting errors (from Tables

Table 6 .
Paleomagnetic data.Declination/inclination: mean characteristic remanent magnetization direction; A95: radius of cone of 95% confidence; kappa: Fisher dispersion parameter around the mean; N: number of samples used in statistics; HC/LC: high/low coercivity component of magnetization.Sites MI02 and MI03 did not provide any results.