Plio-Quaternary uplift of the Iberian Chain (central–eastern Spain) from landscape evolution experiments and river profile modeling

Geomorphology 246 (2015) 48–67

Contents lists available at ScienceDirect

Geomorphology journal homepage: www.elsevier.com/locate/geomorph

Plio-Quaternary uplift of the Iberian Chain (central–eastern Spain) from landscape evolution experiments and river profile modeling E. Giachetta ⁎, P. Molin, V.N. Scotti, C. Faccenna Dipartimento di Scienze, Università degli Studi Roma Tre, Largo S.L. Murialdo, 1-00146 Rome, Italy

a r t i c l e

i n f o

Article history: Received 15 May 2014 Received in revised form 6 May 2015 Accepted 2 June 2015 Available online 6 June 2015 Keywords: Landscape evolution Iberian Chain Uplift Rivers Numerical modeling

a b s t r a c t The Iberian Chain (central–eastern Spain) is characterized by a dome-shaped topography with poorly incised and relatively flat landscape in its interior at an elevation of ~1300 m. A recent regional tectonic event started during or after the Late Pliocene, inducing deep valley incision on the flanks of the Iberian Chain. The rate and timing of this deformation, as well as the magnitude of the uplift-induced erosion processes, still lack a quantitative validation. We used 3D (SIGNUM) and 1-D (river profile) numerical models to investigate the landscape evolution of the Iberian Chain from the Late Pliocene to the present and to study the response of topography and rivers to tectonically induced base level lowering. Our model results are supported by geological, geomorphological and geochronological data. The results of the SIGNUM experiments show that the Iberian Chain topography has experienced regional uplift since ca. 3 Ma with rates between 0.25 and 0.55 mm y−1. In response to this uplift, rivers have incised at an average long-term incision rate of about 0.22 mm y−1 and have progressively captured drainage areas from the chain interior. To further validate the results of landscape evolution modeling, we use a 1-D river incision model to explore the incision history of several rivers draining the Iberian Chain. The results constrain the range of values of bedrock erodibility from 5.8 × 10−5 to 5.0 × 10−4 m0.4 y−1 and predict long-term incision rates similar to those from the SIGNUM experiments. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The topography of a mountain chain is the result of the complex interactions between tectonics (uplift) and erosion (Willett, 1999; Hovius, 2000; Willett et al., 2001). Tectonics generates topography, and surface processes redistribute rock by erosion, transport and sedimentation, resulting in loading and unloading of the crust and mantle lithosphere (e.g., Willett, 1999; Roe et al., 2006; Stolar et al., 2006; Burov and Toussaint, 2007; Cloetingh et al., 2007). A multidisciplinary approach is required for studying topography because of the complexity of the interactions among tectonics, erosion and climate that govern the development of topography in a given setting. Nevertheless, the timing and rate of uplift are usually difficult to quantify from landscape topography and from geological and geochronological data. For this reason many numerical models of landscape evolution have been developed in the last decades that facilitate our understanding of topographic responses to changing tectonic boundary conditions (Tucker and Slingerland, 1994; Braun and Sambridge, 1997; Densmore et al., 1998; Tucker et al., 2001; Carretier and Lucazeau, 2005; Hancock et al., 2011; Refice et al., 2012).

⁎ Corresponding author. E-mail address: [email protected] (E. Giachetta).

http://dx.doi.org/10.1016/j.geomorph.2015.06.005 0169-555X/© 2015 Elsevier B.V. All rights reserved.

As an example of how numerical models can provide a more complete understanding of the tectonics of a mountain chain, we combine a landscape evolution model and a river longitudinal profile model to investigate the uplift history of the Iberian Chain (central–eastern Spain). The Iberian Chain is an ideal location for this type of investigation because the timing and rate of the relief generation are still debated (Gutiérrez et al., 2008; Casas-Sainz and De Vicente, 2009; Scotti et al., 2014 and references therein). The Iberian Chain is an intraplate thrust-belt originally formed during the Late Cretaceous to the Middle Miocene (Álvaro et al., 1979) and has a dome-shaped topography characterized by a poorly incised relict landscape in its interior. This old landscape, standing at a mean elevation of ~ 1300 m, includes one or more nearly flat Miocene–Middle Pliocene erosion surfaces (Birot, 1959; Solé Sabarís, 1979; Peña et al., 1984; Simón, 1984; Gracia-Prieto et al., 1988; Gutiérrez Elorza and Gracia, 1997). In contrast, the flanks of the Iberian Chain are deeply incised by rivers that have been capturing the drainage network of the chain interior as these rivers begin to incise the relict landscape in response to a base level change (Scotti et al., 2014 and references therein). It is generally accepted that a regional tectonic event uplifted the relict landscape and induced rivers to incise the Iberian Chain (Simón, 1984, 1989; Muñoz-Martín and De Vicente, 1998; Gutiérrez et al., 2008; Casas-Sainz and De Vicente, 2009). Geological and geomorphological evidence indicates that this uplift started at or after the Late Pliocene (Scotti et al., 2014 and

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

references therein); however, the rate and timing of this deformation, as well as the magnitude of the uplift-induced erosion processes, still lack a quantitative validation. Here we explore how a numerical landscape evolution model can be used to estimate time and space variation in uplift rate and how it can improve the knowledge on the recent deformation history of an area like the Iberian Chain, where quantitative constraints are poor or missing. At the scale of a mountain range, tectonic processes operate on time scales equivalent to 106 y, whereas the response time of landscape to tectonic or climatic perturbations ranges from 104 to 106 y (Whipple, 2001; Wegmann et al., 2007). Therefore, the spatial scale of the Iberian Chain (~300 km2) and the temporal scale at which the processes of interest (uplift and river incision) operate necessitate the use of numerical experiments to simulate the development of topography and hydrography following regional uplift. The ultimate goal of our experiments is to model the Plio-Quaternary evolution of the Iberian Chain, characterizing the rates and pattern of uplift and erosion that have produced the main features of the present topography. For this purpose we have designed a two-step experiment. In the first step we use a landscape evolution model called SIGNUM (Refice et al., 2012), which has been successfully used to investigate the feedbacks between uplift and erosion in active tectonic settings and the reorganization of drainage network at the orogen scale (e.g. Capolongo et al., 2011; Giachetta et al., 2014). SIGNUM is suitable to simulate the transient response of a synthetic topography representing the Iberian Chain to a tectonic perturbation over 106 year time scale. We ran several numerical models starting with different initial topographies and applied different magnitudes and patterns of uplift rate until we obtained the synthetic topography that best fits the present Iberian Chain. In the second step, we numerically modeled the present day river longitudinal profiles extracted from digital topography to test the time duration of the landscape evolution models, the pattern of uplift, and the stream power parameters that predict the best-fit topography. To test the reliability of the river profile modeling, we also performed a sensitivity analysis to compare the stream power law parameters predicting the best-fit river profiles with those obtained by previous works. The second step of our experiments focuses on the stream longitudinal profiles because among geomorphological systems, rivers are the most sensitive to variations in boundary conditions including tectonics, climate, base level, and lithology. 2. Geological and geomorphological settings The present-day topography of the Iberian Chain, an intraplate double vergent thrust belt, is characterized by low relief surfaces standing at an average elevation of ~ 1300 m (Fig. 1a). These features have been interpreted as a unique one, the so-called Main Planation Surface of the Iberian Chain (Birot, 1959; Solé Sabarís, 1979; Peña et al., 1984; Simón, 1984), or as two distinct nested surfaces (S2 and S3; Gracia-Prieto et al., 1988; Gutiérrez Elorza and Gracia, 1997). The ages of these flat or gently undulated surfaces are from the Miocene to the Late Pliocene according to their morphological correlation with the lacustrine and fluvial deposits of intermontane basins (Gracia-Prieto et al., 1988; Casas-Sainz and De Vicente, 2009) (Fig. 1b). An older erosion surface, the Intra-Miocene Erosion Surface (Peña et al., 1984), or S1 (Gutiérrez Elorza and Gracia, 1997) developed on Mesozoic rocks and is observed on the higher peaks at an elevation of 1700–2000 m. Small ranges, like the Castilian and Aragonian branches and the Maestrazgo Mts. (Fig. 1b), stand a few hundreds of meters above the Main Planation Surface and are interpreted as relics of ancient ridges (Peña et al., 1984; Gutiérrez Elorza and Gracia, 1997). The formation of the high-standing planation surfaces has been linked to two major plate tectonics events occurred in the geological history of the Iberian plate: the collision between Africa and Europe and the opening of the Valencia Through in the western Mediterranean (Casas-Sainz and De Vicente, 2009; Scotti et al., 2014).

49

In the Late Cretaceous, the collision of Iberia with Africa and Europe generated the main mountain ranges at the eastern plate boundary: the Pyrenees (northern Spain), the Iberian Chain (central–eastern Spain) and the Central System (central Spain) (Álvaro et al., 1979; Muñoz, 1992; Teixell, 1998; Casas-Sainz and Faccenna, 2001; Vergés et al., 2002; De Vicente et al., 2007). The Iberian Chain results from the inversion of normal faults that initiated during Mesozoic extension of the Iberian Basin (Álvaro et al., 1979; Guimerà et al., 2004) and were subsequently inverted to accommodate Cenozoic intraplate shortening (Casas-Sainz and Faccenna, 2001). Conversely, the Castilian (RodríguezPascua and De Vicente, 1998) and Aragonian (Ferreiro et al., 1991; Calvo Hernández, 1993; Cortés-Gracia and Casas-Sainz, 1996; Casas-Sainz et al., 1998) branches show a strike-slip component along NW–SE structures with basement-induced positive flower structures (Wilcox et al., 1973; Scotti et al., 2014). The western margin of the Iberian Chain ends against the Central System, a double-vergence intraplate belt, elongated NE– SW and resulting from thick-skinned Tertiary compression involving metamorphic–granitic Variscan basement (De Vicente et al., 2007). The linkage between these ranges is affected by several NW–SE dextral strike-slip faults. The formation of the Pyrenees, Central System and Iberian Chain produced topography that allowed for the development of large compressive basins containing shallow perennial lakes (Calvo et al., 1993; Wright et al., 1997; Casas-Sainz and De Vicente, 2009): the Duero Basin to the NW, the Ebro Basin to the N, the Guadiana Basin to the S, and the Madrid (High Tagus) Basin to the SW (Fig. 1b). The filling of these large endorheic basins started in the Eocene, generating thick continuous succession characterized by lacustrine and alluvial facies (Calvo et al., 1993; Calvo, 2004). Contemporaneously, the uplift of the Central System separated the Duero and Madrid (High Tagus) basins (Calvo, 2004). At local scale, small-size internally drained basins (hereafter referred as intermontane basins) developed in the Iberian Chain interior and were filled by similar but discontinuous continental deposits (Calvo, 2004). The origin of these basins is controversial, being related to either compressional episodes or extensional structures that can be correlated with NE–SW extensional structure along the Mediterranean coast and related to the Neogene opening of the Valencia Trough (Guimerà and Alvaro, 1990; Roca and Guimerà, 1992; Vegas, 1992; Cortés-Gracia and Casas-Sainz, 2000). From the Late Oligocene to the Early Miocene, the large-size basins were characterized by a very discontinuous sedimentation pattern, indicating a varied tectonic activity (Calvo et al., 1993; Calvo, 2004; Casas-Sainz and De Vicente, 2009). Extensional faulting affected the Mediterranean side of the Iberian Chain (the Maestrazgo Mts.) from the Late Miocene to Pleistocene times (Simón et al., 2013). The intermontane basins remained distinct lacustrine areas until the Late Miocene, when sedimentation became widespread, overlapping the sedimentary thresholds (Armenteros et al., 1989; Alonso-Zarza and Calvo, 2000). This phenomenon is recorded by the lacustrine carbonates of the Páramo Formation (Late Miocene–Early Pliocene), the last record of endorheism before the progressive capture of the internally drained basins since the Pliocene (Gutiérrez et al., 2008). A general asymmetric dome-like uplift of the Iberian Chain started in the Late Pliocene (Gutiérrez Elorza and Gracia, 1997; Scotti et al., 2014 and references therein). In the Pliocene–Quaternary, extensional tectonics generated new depressions (Gallocanta, Munabrega, and Jiloca depressions), while others such as the Calatayud and Teruel depressions were reactivated (Cortés-Gracia and Casas-Sainz, 2000; Gracia et al., 2003; Rubio and Simon, 2007; Gutiérrez et al., 2008). The transition from internally drained basins to an exorheic drainage network, mainly driven by capture of the Ebro, Tagus, Jalon and Turia rivers, led to the incision of the ancient landforms and sedimentation of coarser deposits in the marginal areas of the Iberian Chain (Martín-Serrano, 1991; Gutiérrez et al., 1996; Gutiérrez Elorza and Gracia, 1997; Casas-Sainz and Cortés-Gracia, 2002). The uplift drove

50

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

a

4° W

3° W

2° W

1° W

0°

42° N

III’ EB Ca

DB

d yu ta la B.

41° N

TA1 TA2

Jiloca B.

I II 40° N

II’

Te rue B. l

IC

MA1 MA2

MB

Mij

are

sR

.

Elevation m

2500 2000

I’

1500 1000 500 1

III GB 39° N

0

50

100

150

200

Kilometers

b

4°W

3°W

2°W

1°W

0° Legend

CAMEROS MTS.

42°N

Normal faults Thrust faults

DUERO BASIN EBRO BASIN

Paramo Fm. Quaternary

AR

Pliocene

AG NI

O

Miocene Paleogene

AN

CA

ST

M TE

L

Cretaceous

B.

41°N

Jurassic

ILI

AN

S SY

A TR

Triassic

B.

MAESRATZGO MTS.

N CE

Paleozoic Proterozoic

MADRID (TAGUS) BASIN

40°N

MEDITERR. SEA GUADIANA BASIN

39°N

PREBETICS

0

25

50

100

150

200 Kilometers

Fig. 1. Geological and geomorphological settings of the Iberian Chain. a) SRTM DEM of the Iberian Chain showing the location of the studied rivers and the sampling location (circles) of calcareous tufas laying on fluvial strath terraces (after Scotti et al., 2014). The location of three swath-profiles is indicated by white dashed rectangles. IC: Iberian Chain; MB: Madrid Basin; GB: Guadiana Basin; EB: Ebro Basin. b) Geological map of the study area draped on a shaded topographic relief map derived from the SRTM DEM. Description of the outcropping lithologies — Proterozoic: gneiss, schist, marble, volcanic rocks; Paleozoic: quartzite, slate, conglomerate, sandstone, limestone, evaporite, plutonic and volcanic rocks; Triassic: conglomerate, sandstone, limestone, clay; Jurassic: limestone (locally fractured and karstified), dolostone, marl, conglomerate, sandstone; Cretaceous: limestone (locally fractured and karstified), dolostone, marl, sandstone; Paleogene: conglomerate, sandstone, clay; Miocene: conglomerate, sandstone, clay; Pliocene: conglomerate, sandstone, clay; Quaternary: gravel, sand and silt; Páramo Formation (Late Miocene–Early Pliocene): lacustrine limestone.

the headward incision of rivers that deeply incised the Iberian Chain flanks, but could not reach the interior apart from the Tagus, Jalon and Turia rivers, and that led to the progressive capture of the intermontane basins (Gutiérrez et al., 2008). The post-capture configuration of the Iberian Chain drainage was largely controlled by neotectonic

activity, revealed by normal faults affecting Plio-Quaternary deposits (Gutiérrez et al., 2008). The present-day river network dissects the planation surface but few portions of the Iberian Chain interior are still endorheic like the Gallocanta half-graben or part of the Jiloca depression. At present, Upper Pliocene–Quaternary deposits are confined

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

Elevation (m)

mostly in internally drained or recently captured intermontane basins (Rubio and Simon, 2007; Gutiérrez et al., 2008; Lafuente et al., 2011a; and references therein). The planation surfaces as well as the endorheic basins appear not to be affected by the significant base level drop due to the Messinian Salinity Crisis (~ 1500 m below the present sea level; Clauzon et al., 1996). Indeed, this base-level lowering, which lasted for around 600 ka, produced a knickpoint along the Ebro River currently located at 10 km from the present coastline and completely buried by Pliocene deposits (Babault et al., 2006; Loget and Van Den Driessche, 2009 and references therein). In the Iberian Chain the present seismicity is scarce (Casas-Sainz and De Vicente, 2009). It is low to moderate and mostly concentrated in the intermontane basins (Gutiérrez et al., 2009, 2012; Lafuente et al., 2011b; Simón et al., 2012) and along the western margin of the Valencia Through (Perea et al., 2012; Simón et al., 2013). Recently, Scotti et al. (2014) provided important evidence for the timing of late Cenozoic uplift of the Iberian Chain and the spatial pattern of deformation by quantitative geomorphological analysis coupled with geological and geochronological data. The topographic setting of the Iberian Chain is described by three swath profiles (Fig. 2) and a map of local relief (Fig. 3) (Scotti et al., 2014) extracted from a 90 m SRTM DEM, whose resolution is suitable for geomorphological analyses at a regional scale. The swath profiles (Fig. 2) show the dome-shaped topography of the Iberian Chain, characterized by a widespread low-relief surface (planation surface) standing at ~ 1300 m in its interior. Remnants of ancient ranges rise from the planation surface up to 1600 m of elevation. The maximum values of local relief (dashed line in Fig. 2) are mostly around 200 m even in the intermontane basins. The consequent uniform pattern throughout the chain suggests that rivers have been incising equally throughout the range in response to a main tectonic event (Scotti et al., 2014). The maximum amplitude of cyclic sea level fluctuations since the Pliocene reached ~ 100 m (e.g. Waelbroeck et al., 2002). Sea level drawdown could have contributed to enhance river incision but these changes in base level are too small to account for the total amount of incision observed in the Iberian Chain (Scotti et al., 2014).

2000

I

In order to describe the spatial variation in river incision throughout the study area, a map of local relief has been made by arithmetically subtracting two surfaces describing the general configuration of valley bottoms and peak elevations (Fig. 3; after Scotti et al., 2014). The resulting map describes the Iberian Chain as a massif poorly incised by rivers since the values of local relief are low (less than 200 m), more similarly to the Madrid and Duero basins than to all the other mountain belts of the Iberian Peninsula. The lower values are found in the chain interior, but they increase slightly on the flanks. The highest values of local relief are in the eastern sector of the Iberian Chain, along the Mediterranean coast, where active or recent extensional faults are located (Simón et al., 2012). Scotti et al. (2014) analyzed several rivers that drain the Iberian Chain arranged in a radial pattern as a consequence of the dome-like uplift. They extracted river longitudinal profiles (Fig. 4), quantifying the relative river geometry by the steepness index ksn, a measure of channel gradient, and the concavity index θ, the rate of change of channel gradient with drainage area (Hack, 1957; Flint, 1974). The results indicate that most of the rivers draining the Iberian Chain are far from an equilibrium state and probably still adjusting to the regional uplift (Scotti et al., 2014). The concavity index θ has a low average value (0.27). ksn progressively increases from SW to NE as well as from W to E, reaching the highest values in the NE sector of the Iberian Chain (Maestrazgo Mts.; Fig. 1), indicating a local higher uplift rate that recently raised the river heads (Scotti et al., 2014). The highest elevation of the range presents remnants of planation surfaces at high altitude (~2020 m a.s.l.). Scotti et al. (2014) calculated the asymmetry factor AF of 14 catchments, which reflects the plan-view shape of a drainage basin. The dome-like uplift of the massif induced the tilting of the Iberian Chain flanks so that the channels oriented transverse to its axis shift toward the peripheral areas (for example, the Tajo, Jalón, and Jiloca rivers), and those slightly transversal (the Cabriel, Turia, and Guadalope rivers) tend to rotate toward an orientation parallel to the maximum slope of the chain flanks, similarly to the pattern of symmetric basins (the Tajuña, Júcar, Magro, Mijares, Martín, Huerva rivers). The asymmetry of the regional uplift is also recorded by the presentday geometry of the Páramo Formation: it has gentle dips from the I’

IC

WNW

51

ESE

1000

Local relief Mean elevation Minimum elevation Maximum elevation

0

100

0

200

Distance (km)

Elevation (m)

II 2000

E

MB

1000

0

II’

IC

W

200

100

0

300

Distance (km) Elevation (m)

III 2000

NNE

GB

1000

0

III’

IC

SSW

0

EB

100

200

300

Distance (km) Fig. 2. Swath profiles across the Iberian Chain showing the trend of maximum, mean and minimum heights (modified from Scotti et al., 2014). Local relief is indicated by the dashed curve. IC: Iberian Chain; MB: Madrid Basin; GB: Guadiana Basin; EB: Ebro Basin.

52

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

LEGEND Local Relief (m)

Fig. 3. Map of local relief including the along-channel variation in ksn. Modified from Scotti et al., 2014).

top of the Iberian Chain (where it reaches an elevation of ~ 1300 m) toward the center of the Madrid Basin, where its attitude becomes horizontal. The top surface envelope of the Páramo Formation shows an upwarping of the carbonate unit reaching a maximum uplift of 500–600 m with respect to the Madrid Basin and forming a ~ 300 km asymmetric dome with the northern flank steeper than the southern one (Scotti et al., 2014). The setting of the Páramo Formation simply indicates that uplift occurred after the deposition of the formation itself, and therefore post-dates the main compressive phases (Scotti et al., 2014), and provides a minimum uplift rate of ~0.2 mm y−1 that possibly increases up to ~0.6 mm y−1 (Scotti et al., 2014). Roughly at the same time the high elevations of the Iberia Chain experience a change from endorheism to exorheism. Collectively, this evidence suggests that the rivers in this region should respond to a tectonic input, the uplift of the Iberian Chain, and variation in base level related to a change from lake to sea. To constrain the timing of the recent uplift of the Iberian Chain, Scotti et al. (2014) compare the morphometry of rivers draining the Chain with that of some rivers draining the eastern flank of the Central System. The latter are close to the equilibrium profiles suggesting that the uplift of the Iberian Chain could be younger than 5 Ma, when a rapid cooling phase started in the Central System (De Bruijne and Andriessen, 2002; Ter Voorde et al., 2004). Scotti et al. (2014) estimated incision rates by dating of calcareous tufas lying on top of strath fluvial terraces along the high courses of the Martín, Ruguilla, and High Tagus rivers (Fig. 1b). By the obtained ages corresponding to MIS 1, 3, and 5, they calculated a similar incision rate (~0.6 mm y−1) from the Late Pleistocene to the present throughout the chain interior (Table 1). This is a further evidence that the rivers have been incising in response to a single regional tectonic event (Scotti et al., 2014) and secondary to the drastic base level fall and increase in river discharge due to the diachronous capture of the endorheic basins that occurred during the Late Pliocene–Late Pleistocene (Gutiérrez et al., 2008) and the lithology of river bedrock, both poorly consolidated sediments (the Martín and Ruguilla rivers) and stronger limestone (the Tajo River).

3. SIGNUM landscape evolution model Numerical modeling of landscape evolution has been widely used to investigate the interactions of tectonic uplift and erosion at regional scale. Modeling of longitudinal river profiles has been used to reconstruct the histories and patterns of regional and local tectonic uplift. In this study we designed a two-step modeling experiment: in the first step we reconstruct the paleotopography of the Iberian Chain during the Pliocene and the Quaternary (Section 3) using a landscape evolution model SIGNUM (Refice et al., 2012; Giachetta et al., 2014) and in the second step, we model the stream profile evolution of rivers draining the Iberian Chain (Section 4). The numerical model, SIGNUM, represents the land surface by a set of points with three-dimensional coordinates, connected by edges to form a triangulated irregular network (TIN). 3.1. Geomorphological processes Landscape evolution models reproduce the changes of the Earth surface through the interaction of three broad classes of processes which we indicate here with the terms of surface uplift, diffusion, and channeling respectively. These surface processes have been implemented in a large variety of ways, and the study of their interactions, rates, and forms, constitutes most of the current research effort in the field. Changes in the elevation of the topographic surface are described by the equation of mass-conservation (e.g. Tucker et al., 2001): ∂z=∂t ¼ U−K Q m Sn −kd ∇2 z

ð1Þ

where z is the elevation of a point on a TIN; t is time; U is rock uplift rate; K is a dimensional coefficient representing bedrock erodibility (e.g., Howard et al., 1994; Sklar and Dietrich, 1998; Whipple and Tucker, 1999; Tucker and Bras, 2000); m and n are dimensionless exponents, generally fixed to 0.5 and 1, respectively (e.g., Tucker and Whipple, 2002); Q is water discharge often approximated by the expression Q = P · A, where P is effective precipitation and A is drainage area; S is local channel slope; kd is a constant of diffusion and ∇2z is

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

53

Fig. 4. Longitudinal profiles of rivers draining the Iberian Chain (see Fig. 3 for location), including inset diagrams of log slope vs. log drainage area (modified from Scotti et al., 2014). The profiles and diagrams were extracted from the SRTM DEM using the Stream Profiler Tool (available at www.geomorphotools.org). Gray arrows indicate the location of major knickpoints, while black arrows show the location of main knickpoints related to rock-type changes. θ and ksn are also indicated for each river. D: dams.

surface curvature. The second term on the right-hand side of Eq. (1) represents the “stream power law” for the fluvial incision model (e.g., Howard, 1994, 1997; Stock and Montgomery, 1999; Whipple and Tucker, 1999). The last term in the right-hand side of Eq. (1) describes hillslope processes (Dietrich et al., 1995; Roering et al., 1999). SIGNUM calculates the height changes at each time step, using a forward-in-time explicit resolution scheme to solve all terms in Table 1 Incision rates calculated from the geochronological data. Samplea

U/Th age

Incision rate mm y−1

TA1 TA2 MA1 MA2

89 ka 105 ka 91 ka 41 ka

0.69 ± 0.04 0.56 ± 0.04 0.66 ± 0.05 0.62 ± 0.03

a

Sampling locations are indicated in Fig. 1 (Scotti et al., 2014).

Eq. (1). Time is simulated as a sequence of discrete time steps and the main program updates all the geomorphological variables (e.g. elevation, slope, and drainage area) at each time step. To preserve the numerical stability in the solutions of Eq. (1), each time step is calculated by each process module taking into account the Courant–Friedrichs–Lewy condition (hereafter CFL condition) (Courant et al., 1967). The minimum value among all the CFL time steps calculated by the process modules is used as a time step for the next iteration. We refer the reader to Refice et al. (2012) for more details on the calculation of the CFL limit for each process module and on the algorithms and structure of SIGNUM. 3.2. Sedimentation in closed depressions The large-size internally drained basins in eastern Iberia and their change to exorheism in the Pliocene–Quaternary represent a firstorder feature of the landscape evolution of the Iberian Chain. Therefore,

54

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

a numerical treatment of closed topographic depressions is required for this study; the role of endorheic basins must be explicitly considered in the drainage network development and its interaction with tectonic uplift in the model simulation. We developed a simplified method to simulate sediment deposition in a closed basin: since the incision takes place if the steepest descent slope is negative, all points on the TIN having a steepest-descent slope Sst ≥ 0 (meaning that all neighbors points have elevation z ≥ zi) are flagged as lake-points. Thus, no incision is allowed by the model in these points and the sediment flux, Qs, from an upstream cell is deposited in the lake-point cell, i. The model checks the elevation of the point, i, and if it is higher than the minimum elevation calculated for all its neighbors, an amount of sediment is redistributed to the neighbors. This method produces a horizontally leveled surface interpolating the filled lake-points. If an internally drained basin is overfilled by sediments, it can be captured by an external river breaching the divide of the closed depression and an outlet of the lake is produced. 3.3. Model setup A 300 km2 square synthetic domain was used to fit the present-day area of the Iberian Chain. The landscape is numerically represented by a randomly sampled and triangulated set of points with x, y, and z coordinates. The spacing of points is 2 ± 0.4 km for all runs. We used the geological and geomorphological data by Scotti et al. (2014) described in Section 2 to implement robust initial conditions for the model and test the adjustment of the simulated topography and hydrography to a regional tectonic uplift. We tested three initial topographies (see Fig. 5): - topography 1 (Fig. 5a): the initial topography is a flat surface that simulates a wide planation surface at sea level; - topography 2 (Fig. 5b): the initial topography is more complex, including lacustrine basins surrounded by low elevated ranges and a planation surface. We consider that the planation surface is connected to the depositional surface of the lacustrine basins (Gutiérrez Elorza and Gracia, 1997), hence they are indistinguishable in Fig. 5. A mountain belt, located at the western boundary, oriented NE–SW and reaches a maximum elevation of ~ 1000 m and a half-width of 50 km, corresponds to the present Central System. Along the northeastern boundary and in the central sector of the synthetic landscape, two parallel ranges, oriented NW–SE, correspond to the modern Aragonian Branch and Castilian Branch of the Iberian Chain. These two ranges, 50 and 70 km wide respectively, and with initial maximum elevation of ~ 300 m, isolate an internally drained lacustrine area corresponding to the present Calatayud and Teruel basins (Armenteros et al., 1989; see Fig. 3). To the E a low topographic barrier, oriented NNE–SSW, with a maximum elevation of ~300 m and a width of 50 km, corresponds to the coastal ranges (Cortés-Gracia and Casas-Sainz, 2000); - topography 3 (Fig. 5c): the initial topography is equal to topography 2, except for a mountain range, ~700 m high, 100 km wide, located in the NE sector of the domain, and corresponding to the present Maestrazgo Mts. Here, planation surfaces are located on the mountain top (Peña et al., 1984; Gutiérrez Elorza and Gracia, 1997). Initial topographies 2 and 3 including the Madrid (High Tagus)/ Guadiana basin were set according to the history of the main tectonic and erosive episodes in the eastern Iberia derived from previous studies (see Section 2). The Duero and Ebro basins, as well as the Mediterranean Sea, were set as open boundaries of the synthetic domain (Fig. 5). The Madrid (High Tagus)/Guadiana basin is represented as a wide flat surface, with average elevation a few meters above sea level. Water and sediment can flow out across the four boundaries of the synthetic domain. Two exceptions are the sectors corresponding to the Central System and the Prebetics.

The regional uplift was simulated by a vertical displacement applied to all the points. We tested two uplift patterns characterized by different spatial distributions of uplift rates to obtain the final topography that best fits the present landscape of the Iberian Chain. First, we use a simple dome-like uplift pattern characterized by a maximum uplift rate in the NE sector of the synthetic domain (Fig. 6a). Second, we use a more complex uplift pattern subdividing the synthetic domain into six sectors (Fig. 7a), in which uplift rates are assumed to be uniform, except for sector 1 characterized by an uplift gradient (Fig. 7b). The six sectors have been identified according to topographic, geological, geomorphological and chronological data. Sector 1 has increasing uplift rates toward the Iberian Chain interior, as described by the attitude of the Páramo Formation: it indicates a maximum uplift of around 500 m with respect to the present Madrid Basin elevation (Scotti et al., 2014 and references therein). This maximum uplift rate is equal to that of sector 2. Sector 3 includes the Central System that has been affected by an uplift of 3.5 ± 0.7 km in the last 7 My (De Bruijne and Andriessen, 2002). Sector 4 corresponds to intermontane basins (the Teruel and CalatayudMontalbán basins and the upstream Turia R. valley; López-Martínez et al., 1987; Anadón et al., 1990; Anadón and Moissenet, 1996; Alcalá et al., 2000; van Dam and Sanz Rubio, 2003) (Fig. 1a). Here the uplift rate is slightly lower than the surrounding areas, generating tectonic subsidence (Gutiérrez et al., 2008). Sector 5 includes the Maestrazgo Mt. area, where topographic and geomorphological data indicate that it is the most uplifted portion of the Iberian Chain. Sector 6 includes the Guadiana Basin on the far end of the dome-like uplift deformation of the Iberian Chain. A power–law relationship between channel slope and drainage area can be derived from Eq. (1) at steady-state (∂z/∂t = 0) and with uniform U and K. Under these conditions, channel concavity is given by the ratio of stream power model parameters m/n (Howard, 1980). Considering that the average river concavity index θ = ~ 0.27 (Scotti et al., 2014), we implemented a channeling model using m = 0.25 − 0.45 and n = 1 as parameters of the stream power model of Eq. (1). We assumed an effective precipitation P = 1000 mm y−1, and P, K, m, n and kd were set uniform in space and time. The value of P results from the mean between the minimum (300 mm y−1) and maximum (1700 mm y−1) precipitation rates estimated for the period from the Miocene to the Pleistocene in the Iberian Peninsula (van Dam and Weltje, 1999; van Dam, 2006; Fernández et al., 2007; Jiménez-Moreno et al., 2010; AEMET-IM, 2011). We performed several calibration runs by iteratively varying the values of the erosion process parameters and uplift rates. Tables 2 and 3 describe the ranges of parameters used to calibrate the erosion processes on real data. We used different simulation time intervals, ranging between a maximum duration of 4 My and a minimum duration of 2 My. The time step also varies depending on the set of erosion parameters and according to the CFL conditions. In order to identify the best values of the geomorphological and tectonic parameters predicting the best-fit topography, we compared first order features, such as the elevations and drainage patterns of the real and simulated topographies. 3.4. Landscape evolution modeling results The best fitting set of erosion parameters and pattern of uplift should satisfy all the following conditions: 1) the values of the elevations of the simulated topography fit to those from the SRTM DEM; 2) the pattern of the drainage network is similar to the present hydrography; 3) the timing of the transition from endorheism to exorheism is consistent with data from literature; and 4) the incision rates reasonably agree with incision rates estimated by 230U–Th datings (Scotti et al., 2014). The model results using the dome-like uplift produced unsatisfactory predictions (e.g., Fig. 6b). Although the synthetic landscape shows a more uplifted area in correspondence to the Maestrazgo Mts. and explains the radial drainage pattern, it presents high mismatch concerning

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

55

m

a Duero Basin

Ebro Basin

m0

Tagus River

km

km

Mediterranean Sea

b

m

Duero Basin

Ebro Basin AB CB

CS

m

MB

Tagus River

GB

km

km Mediterranean Sea

m

c

Duero Basin

Ebro Basin AB CB

CS

m

MB

Tagus River

MM

GB km

km Mediterranean Sea

56

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

a

b

Fig. 6. Setup of uplift pattern 1 and modeling results. a) Asymmetric dome-like uplift pattern with a maximum uplift rate in the north-eastern sector of the synthetic domain. b) Topography obtained after 4 My of landscape evolution using the uplift pattern in Fig. 6a and the initial topography of Fig. 5a.

some topographic features of the Iberian Chain (e.g. intermontane basins) and the timing of drainage reversal. Applying the uplift pattern with six sectors (Fig. 7a, Table 3) to topographies 1 and 2 configurations (Fig. 5a, b), we obtained the model results that fit partially the presentday landscape, but produce a higher mismatch relative to the model run using topography 3 (Fig. 5c). Indeed, the results from the topography 3 configuration appears to be the best including higher elevations in correspondence to the Central System and the Maestrazgo Mts., the flat high standing landscape in the Iberian Chain interior, and the depressions in correspondence to the Madrid (High Tagus), Guadiana,

Calatayud, Jiloca, and Teruel basins (Figs. 1 and 8). This result also shows similarities between the present and modeled drainage networks (Fig. 8). In the Madrid (High Tagus) and Guadiana basins (sector 1 in Fig. 7) the southwestern part of the lake surface changed into an alluvial plain at around 2.8 Ma, since headward incision proceeded from the open boundary in the SW that corresponds to the present location of the Tagus River. The values of the model parameters predicting the best-fit topography are listed in Table 2. The model predicted an average longterm incision rate of ~ 0.22 mm y−1, which can be considered a good

Fig. 5. Three synthetic domains (300 × 300 km in size) showing different initial topographies and representing the Iberian Chain topography during the Late Pliocene (~3 Ma); a): flat surface at sea level; b) and c): lacustrine basins surrounded by low elevated ranges; AB: Aragonian Branch; CB: Castilian Branch; CS: Central System; GB: Guadiana Basin; MB: Madrid– High Tagus Basin; MM: Maestratzgo Mountains.

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

57

a Ebro Basin

Duero Basin

4 3

2

2

Madrid Basin

1

5 4 Mediterranean Sea

6

b

m y-1 Ebro Basin Duero Basin

4

2

3

2 5 4

1

km

6 Tagus Basin

Mediterranean Sea

km

Fig. 7. Setup of uplift pattern 2. a) Topography of the Iberian Chain and surrounding areas subdivided in sectors characterized by different uplift rates. Sector 1: uplift rates increasing from the Madrid (High Tagus) Basin toward the Iberian Chain interior; sector 2: Aragonian Branch; sector 3: Central System; sector 4: intermontane basins (Teruel, Calatayud-Montalbán and the upstream Turia R. valley); sector 5: Maestrazgo Mts.; sector 6: Guadiana Basin. b) Complex uplift pattern implemented in the numerical model.

approximation of the incision rate estimated by 230U–Th datings (Fig. 1b; Scotti et al., 2014). We used the best-fit topography to extract two maps (local relief and incision rate) regarding the spatial distribution of fluvial incision. The synthetic local relief (Fig. 9a) was calculated using the same procedure of the map in Fig. 3 (see Scotti et al., 2014). The values of local relief show that the maximum incision (up to 400 m) occurs along the flanks of the synthetic Iberian Chain, except

for the Madrid (High Tagus)/Guadiana basins sector where the lowest values are located. We generated an incision rate map using the modeled incision rates for the last time step (2–0 ky), in which all the highest values (0.6–1 mm y−1) occur along the Iberian Chain flanks (Fig. 9b). In the Madrid (High Tagus) and Guadiana basins and the intermontane basins the values are very low (0.05–0.2 mm y−1).

58

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

Table 2 Erosion parameters and time intervals used in calibration runs and values implemented for the best-fit topography (Fig. 8). Model parameters

Investigated

kd (m2 y−1) K (m1–2m y−1) m/n Simulation time (My) Long-term erosion rate (mm y−1)

Best-fit topography

Min

Max

1 × 10−3 1 × 10−5 0.25 2.0 0.13

1 × 10−2 3 × 10−4 0.45 4.0 1.1

5 × 10−3 1.25 × 10−4 0.3 3.0 0.22

ei ¼

4. Modeling of river profiles 4.1. Numerical approach To test the reliability of the landscape modeling results using a different source of data and methodological approach, we performed numerical experiments of river longitudinal profile evolution. We intended to determine if K, m, n and uplift rates that produce the bestfit topography from the SIGNUM simulations are similar to those predicted from model analyses of river profiles extracted from the 90 m SRTM DEM. The height variations along a river profile were modeled using the following equation (Whipple and Tucker, 1999, 2002; Tucker and Whipple, 2002; Roberts and White, 2010; Roberts et al., 2012a,b): ∂zðx;t Þ =∂t ¼ U ðx;t Þ −Eðx;t Þ

ð2Þ

describing the changes in elevation at location x along the river profile and at time t, where U represents the bedrock uplift at x, and E describes river incision. E was calculated according to the stream power incision model (e.g. Howard and Kerby, 1983; Howard, 1994; Whipple and Tucker, 1999): E ¼ KAm Sn

ð3Þ

where A is local drainage area and S is channel slope. Incision rates are proportional to the parameter K in Eq. (1). Several studies (e.g. Howard, 1980; Howard and Kerby, 1983; Whipple and Tucker, 2002; Gasparini et al., 2007) have demonstrated that if all independent parameters in Eq. (1) are constant through time, the incision rates tend to progressively balance uplift rates along the river profile up to steady-state condition. We use the detachment-limited model of Eq. (3), which assumes steady, uniform flow and in which sediment transport exceeds sediment supply (Howard, 1994). In the detachment-limited model, incision is a power–law function of the shear stress applied to the bed. This model allows the long-term preservation of knickpoints corresponding to local variations in bedrock erodibility, changes in uplift rate along the stream or sudden base level falls induced by stream captures (Hack, 1973; Howard, 1994; Bishop, 1995; Whipple and Tucker, 1999; Mather, 2000; Stokes et al., 2002).

Table 3 Range of uplift rates used in calibration runs and values implemented for the best-fit topography (Fig. 8). Sectors of the Iberian Chain (see Fig. 7)

1 2 3 4 5 6

We used the spatial pattern of uplift rates that predict the best-fit topography in the SIGNUM experiments to model the values of K to produce knickpoints at the locations of the real ones in our river profile modeling. The values of K estimated from the distribution of the knickpoints were compared with those predicting the best fitting topography from SIGNUM. As an additional check we compared modeled K-values with those determined from measured incision rates from geochronological data (Scotti et al., 2014). The average mismatch between the real and modeled profiles was calculated by:

Investigated uplift rates (mm y−1) Min

Max

0.125–0.4 0.3 0.3 0.15 0.3 0.125

0.3–0.7 1 1 0.4 1 0.5

Best-fit topography (mm y−1)

0.25–0.5 0.5 0.55 0.35 0.55 0.3

XN i¼1

ðzmi −zri Þ=N

ð4Þ

where ei refers to mismatch at the point, i, along the profile, zmi is the modeled height of the point, zri is the height of the point extracted from the DEM, and N is the total number of points. We assumed that the best-fit profile was obtained when the average mismatch, e, between the actual and modeled river profile was lower than the threshold value of 1 m. 4.2. Setup of river profile model We chose four rivers draining different sectors of the Iberian Chain. The Tagus and Tajuña rivers originate at the top of the Iberian Chain and flow along its south-western flank down to the Madrid (High Tagus) basin, an area characterized by the tilting related to the domelike uplift of the Chain (Scotti et al., 2014; sector 1 in the 3D modeling, Fig. 7). The Mijares River originates in the Maestrazgo Mts. (sector 5 in Fig. 7), crosses the eastern flank of the Iberian Chain, and flows into the Mediterranean Sea. The Martín River drains the northern flank of the chain, originating in a small intermontane basin and adjoins the Ebro River at the northern boundary of the Iberian Chain. The Martín River flows into an area characterized by uniform uplift (sector 2 in Fig. 7; Scotti et al., 2014). To model the profiles of the Tagus, Tajuña, Mijares and Martín rivers according to Eq. (2), we extracted the drainage area and elevation data from the 90 m SRTM DEM of the Iberian Chain. The values of m and n were derived from the average θ value calculated from the morphometric analysis of the rivers draining the study area (Scotti et al., 2014). In the modeling, we adopted the same initial conditions of the 3D modeling and then applied a tectonic input to each profile. In more detail, U was set for each river according to the uplift pattern shown in Fig. 7a. ∂z/∂t = U (Eq. (2)) at the lower boundary when z N 0 at x = 0; otherwise we set ∂z/∂t = 0 at x = 0. In Eq. (2) S is variable through time, depending on elevation points, and it was calculated for each point using the nearest downstream points. We tested K, initially set as uniform along each profile, in the range 10−7 to 10−2 (Stock and Montgomery, 1999; Whipple et al., 2000; Kirby and Whipple, 2001). Then we chose the value of K predicting the minimum average mismatch e between the actual and modeled river profiles (Eq. (4)). Successively, K was varied along the profile by trial and error until the average mismatch e became lower than the threshold value of 1 m (Pelletier, 2007). In all simulations, the elevation along the synthetic river profile was given by Eq. (2), solved in a forward explicit, finite-difference scheme. The time step varied dynamically to satisfy the CFL condition (see Section 3.1). The total duration of each iteration is 3 My. The slope exponent n in the stream power law has been set to 1, considering that the incision rate is linearly proportional to the unit stream power (Howard and Kerby, 1983; Whipple and Tucker, 1999). To help the comparison of the resulting K values with those from previous studies, we performed numerical experiments setting the stream power parameter m = 0.3 and 0.4. Indeed, m b 0.3 in the stream power model is inconsistent with basic hydraulic relationships of channel width and bank full discharge, while common values of m fall in the range 0.3 − 0.6 (Stock and Montgomery, 1999; Whipple, 2004).

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

2.5 Ma Duero Basin

59

1.0 Ma Ebro Basin

Tagus River Mediterranean Sea

2.0 Ma

0.5 Ma

1.5 Ma

0 Ma

Fig. 8. Time evolution of the best-fit topography: time steps (500 ky each) of the numerical model showing the progressive development of the best-fit synthetic topography and drainage network. Rivers are represented with line width proportional to their modeled drainage area.

4.3. Results of river profile modeling The best-fit profiles resulting from river profile modeling are shown in Figs. 10 and 11. Fig. 12 shows the results regarding the average mismatch between the real and modeled profiles and K modeled along the profile. The best-fit K values obtained using m = 0.3 are acceptable for the outcropping rock-types in all the modeled rivers (Table 4). Moreover, the value m = 0.3 is expected for rivers in which incision occurs during high discharges (Stock and Montgomery, 1999), such as the Iberian Chain rivers that have a regime ruled by the Mediterranean climate (AEMET-IM, 2011). The good match of modeled and measured bedrock erodibility is strengthened by the little variation of modeled K values in the Tajuña, Tagus, Mijares and Martín rivers, corresponding to sedimentary rocks. The only exception is represented by a 7 km reach along the Martín River profile where metamorphic rocks outcrop (quartzite and slate; IGME, 1994). The K values obtained using m = 0.4 are out of the typical range of the outcropping lithologies (Stock and Montgomery, 1999; Whipple et al., 2000; Kirby and Whipple, 2001; see Table 5). In the upstream part of the Tagus River, at 0–100 km from the divide, the bedrock consists mostly of carbonate rocks, locally fractured and affected by karst corrosion (Gutiérrez Elorza and Gracia, 1997; Casas-Sainz and Cortés-Gracia, 2002). According to both geological and geomorphological data (Scotti et al., 2014 and references therein) and the result of landscape evolution simulations, the uplift rate is constant here. In this portion of the profile, K was assumed to be between 5.8 × 10−5 and 2.4 × 10−4 m0.4 y−1, and the lower K values represent more resistant lithologies among the carbonate rocks (dolomite),

whereas the higher K values correspond to the fractured, more erodible limestone (Fig. 10a). In the middle part of the Tagus River profile (between 100 and 200 km), the higher K values obtained by river modeling (2.0–2.5 × 10− 4 m0.4 y− 1) are related to highly erodible sedimentary rock types. In this part, the Tagus River exits from the high standing flat landscape and flows down to the Madrid (High Tagus) Basin, crossing the area affected by the uplift gradient. A regular increase of K values is obtained in the downstream part of the Tagus River, ranging between 2.7 × 10− 4 m0.4 y− 1 for clayey lithologies (between 300 and 200 km) and 2.0 × 10− 4 m0.4 y−1 for alternating sedimentary rocks (300–490 km) (Fig. 10a). Modeling the Tajuña River, the predicted values of K in the upstream profiles portion (0–100 km) range between 1.3–3.0 × 10−4 m0.4 y−1 for carbonate rock types and 1.5–2.1 × 10−4 m0.4 y−1 for alternating clay and sandstone (Fig. 10b). The downstream portion of the profile (310–400 km) is similar to that of the Tagus River. In the Mijares River modeling (Fig. 11a), where we assume that the uplift is uniform along the profile, the K values vary between 6.0 × 10−5–5.0 × 10−4 m0.4 y−1: the higher variability of K corresponds to the alternation of dolomite, limestone and conglomerate along the Mijares River more than along the Tagus and Tajuña rivers. Most of the irregularities in the Mijares River profile are partially due to NNE– SSW normal faults, related to the opening of the Valencia Through started in the Oligocene. These faults generated small depressions filled with alluvial deposits (Simón et al., 2012 and references therein; Fig. 11a). The modeling for the Martín River predicted values of K between 9.0 × 10−5 and 1.75 × 10−4 m0.4 y−1. The lowest K value was calculated

60

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

a Duero Basin

Ebro Basin

Local Relief (m) 0 - 50 51 - 100 101 - 150 151 - 200 201 - 250 251 - 300 301 - 350 351 - 400

Tagus River

Mediterranean Sea

b Duero Basin

Tagus River

Ebro Basin

Erosion Rate (m y-1)

Mediterranean Sea

Fig. 9. Topographic analysis of the simulated landscape. a) Local relief map draped on a shaded best-fit topography (Fig. 8); b) Map of incision rates calculated for the best-fit topography during the last time-step of simulation (2–0 ky).

for a reach of the upstream part of the Martín River where metamorphic rocks outcrop. The large knickzone between 24 and 31 km of distance from the divide corresponds to this metamorphic bedrock (Fig. 11b). In the intermediate part of the profile between 31 and 75 km, the bedrock is mainly represented by limestone and modeled K varies between 1.05 and 1.25 × 10−4 m0.4 y− 1. Slightly higher values (1.2–1.35 × 10−4 m0.4 y−1) were obtained along the downstream part of the profile, where the outcropping rocks are alternations of sandstone and clay in the Ebro River valley. The long-term average modeled incision rates calculated from the best fitting river profiles (Table 4) are consistent with the incision rates estimated by 230U–Th dating of strath terraces (Table 1; Scotti et al., 2014). 5. Discussion The results of SIGNUM modeling for the Iberian Chain are consistent with the present-day topography and drainage pattern and support

the hypothesis of Plio-Quaternary uplift (Scotti et al., 2014). For example, the uplift pattern used in the landscape evolution model spanning over 3 My allowed us to obtain a synthetic topography presenting many features similar to the modern one (see Figs. 1 and 8 for comparison). According to the modeling of initial topography 3 (Fig. 8), the landscape evolution is mostly driven by a regional trend of the uplift that increases toward the northeastern sector of the Iberian Chain, reaching its local maximum in the Maestrazgo Mts. (sector 5 in Fig. 7a). This uplift pattern led to the capture of the synthetic Madrid (High Tagus) and Guadiana basins and the headward incision of the paleo-Tagus River. Moreover, on the southern flank of the synthetic Iberian Chain (sector 1; Fig. 7), the increasing uplift rate induced a tilting of the deposition surface of the paleo-lake deposits. In the real landscape, this tilting is recorded by the attitude of the Páramo Formation (Scotti et al., 2014 and references therein). The uplift gradient in sector 1 was interrupted by sector 4, which corresponds to

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

61

a Tagus R. modeling Profile extracted from DEM Simulated profile

UNIFORM UPLIFT 0.5 mm y-1 UPLIFT GRADIENT Max 0.5 mm y-1 – Min 0.25 mm y-1

Lithology

b Tajuña R. modeling

Profile extracted from DEM Simulated profile

UPLIFT GRADIENT Max 0.45 mm y-1 Min 0.25 mm y-1

Lithology LEGEND Evaporite, marl, clay, sand

Conglomerate

Conglomerate, sand, clay

Evaporite, sand, clay

Clay

Limestone

Sandstone

Marl, clay, sand, conglomerate

Dam

Fig. 10. Modeling of the longitudinal profiles of the Tagus and Tajuña rivers. a) Longitudinal profile extracted from the SRTM DEM and simulated profile of the Tagus River. b) Longitudinal profile extracted from the DEM and the simulated profile of the Tajuña River. The solid and open arrows indicate the uplift pattern implemented in SIGNUM that predicts the best-fit topography in Fig. 8; the bar below the profile shows the rock types outcropping along the river course, derived from the 1:100.000 geological maps (IGME, 1994).

62

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

a

Mijares R. modeling

Profile extracted from DEM Simulated profile

UNIFORM UPLIFT = 0.55 mm y-1

Lithology

b

Martìn R. modeling

Profile extracted from DEM Simulated profile

UNIFORM UPLIFT = 0.5 mm y-1

Lithology LEGEND Sandstone, clay, evaporite

Conglomerate

Clay and evaporite

Slate, quartzite, calcareous flysh and greywacke

Sandstone

Conglomerate, sandstone and clay

Limestone

Marl

Dam

Fig. 11. Modeling of the longitudinal profiles of the Mijares and Martín rivers. a) Longitudinal profile extracted from the SRTM DEM and the simulated profile of the Mijares River. b) Longitudinal profile extracted from the DEM and the simulated profile of the Martín River. The solid arrows indicate the uplift pattern implemented in SIGNUM that predicts the best-fit topography in Fig. 8; the bar below the profile shows the rock types outcropping along the river course, derived from the 1:100.000 geological maps (IGME, 1994).

the intermontane basins of Calatayud, Jiloca and Teruel (Figs. 1a and 7). To the north, the uplift reached a value of 0.5 mm y−1 in correspondence to the Aragonian branch. In the southern and southeastern

borders of the synthetic domain, the lower uplift rate of the northernmost end of the Prebetics (0.3 mm y−1 in sector 6) and of the southeastern portion of the Iberian Chain (0.35 to 0.45 mm y−1 in sector

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

a

b 0

10

−1

10

−2

−3

10

−4

10

−5

10

−6

10

−7

10

Tagus R.

−8

10

0

1000 800 600 400 200 0 −200 −400 −600 −800 −1000 0

mismatch (m)

10 K (m0.4 y−1)

63

5

10

15

20

25

30

35

Tagus R. 5

10

0

10

10

mismatch (m)

10

K (m

0.4 −1

y )

−3

10

−4

10

−5

10

−6

10

Tajuña R.

10

−8

0

200 0 −200

5

10

15

20

25

30

Tajuña R.

−600 −800 0

35

5

10

0

2000

−1

1500

−2

1000

10

mismatch (m)

10

−3

10

−4

10

−5

10

−6

−7

Mijares R.

10

−8

25

30

35

500 0 −500

5

10

15

20

25

30

35

Mijares R.

−1500

−2000 0

5

10

15

20

25

30

35

Iteration

Iteration 0

10

2000

−1

10

1500

−2

mismatch (m)

10

−3

10

−4

10

−5

10

−6

10

1000 500 0 −500

−1000

−7

Martìn R.

10

5

10

15

20

25

30

Martìn R.

−1500

−8

0

20

−1000

10

0

15

Iteration

10

K (m0.4 y−1)

35

400

Iteration

K (m0.4 y−1)

30

−400

−7

10

25

600

−2

10

20

800

−1

10

15

Iteration

Iteration

35

Iteration

−2000 0

5

10

15

20

25

30

35

Iteration

Fig. 12. Relationship between modeled K and river profile mismatch. a) Plots showing the iterations performed to find the best-fit values of bedrock erodibility K for each river; solid circles indicate the mean of the modeled values of K; bars represent the range of tested K for each iteration. Y-axis is logarithmic to better display the data. b) Plots showing the mismatch between the predicted profile and the profile from the SRTM DEM; average mismatch values corresponding to each simulated range of K in (a) are indicated with the solid circles; bars indicate the range of ±1σ dispersion from the average mismatch.

1) avoided the early capture of the internal basins by rivers flowing to the Mediterranean Sea. The model suggests that the uplifting divides at the boundary between the synthetic Guadiana Basin and the

Mediterranean coast (Prebetics) (Fig. 7) have been preventing the capture of the internally drained basins for ~ 1.5 My since the onset of the uplift. This can be partially related to the small size of the drainage

64

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

Table 4 Best-fit K values, average mismatches, standard deviations and long-term incision rates calculated for the best-fit river profiles. River

Bedrock erodibility K (m0.4 y−1) (range along profile)

Average mismatch e (m) ± 1σ

m, n

Long-term incision rate over 3 My mm y−1 ± 1σ

Tagus Tajuña Mijares Martìn

5.8 × 10−5–2.8 × 10−4 1.3 × 10−4–3.0 × 10−4 6.0 × 10−5–5 × 10−4 9.0 × 10−5–1.75 × 10−4

0.76 ± 6.22 0.93 ± 2.81 0.05 ± 2.98 0.69 ± 3.6

0.3, 1 0.3, 1 0.3, 1 0.3, 1

0.23 ± 0.05 0.2 ± 0.02 0.4 ± 0.15 0.34 ± 0.1

area characterizing the rivers draining the Mediterranean margin of the synthetic domain. After 1.5 Ma, incision overcame the uplifting topographic barriers (Fig. 8), allowing the synthetic Tagus, Turia and Jiloca rivers to start capturing the interior of the Iberian Chain and the Guadiana Basin. The sediment overfilling simulated in the internally drained depressions probably partially contributed to their capture. Around 1–0.5 Ma, the rivers reached the intermontane basins. After the Teruel and Calatayud depressions were captured by the paleo-Turia and paleo-Jalon rivers, an axial stream developed in the central part of the basins, connecting the original centripetal hydrography of the previously closed depression. This late capture of the intermontane basins resulted in the topography of the Iberian Chain interior characterized by low relief areas, and deeply incised valleys on the Mediterranean and Ebro sides (Fig. 8). The role of the uplift pattern is confirmed by the river profile modeling. The concavity of the modeled Tagus, Tajuña, Mijares and Martín rivers is θ = m/n = 0.3, whereas the real ones are 0.42, 0.28, 0.28 and 0.17 respectively. Despite these differences, the mismatch between the modeled and real profiles is very small (e b 1 m; Table 4). We attribute the differences between the Tagus and Tajuña river concavities to the main trunk orientation with respect the uplift gradient. The Tajuña River channel has been always perpendicular to it; indeed, its concavity is closest to 0.3. In contrast, the most upstream portion of the Tagus River is parallel to the uplift gradient, although the downstream portion becomes more perpendicular to it (Figs. 4 and 7). In addition, the bedrock in the upper reach of the Tagus River is more resistant than that in the lower reach, as shown by the range of K calculated from the profile modeling. The high values of ksn of rivers draining the Mediterranean and Ebro sides of the Iberian Chain (Fig. 3) mark the boundary between the deeply incised flanks of the Iberian Chain and the relict high standing surface preserved in the Iberian Chain interior. The location of this boundary is indicated in the swath profiles (Fig. 2), where the local relief shows the difference in river incision between the flanks and interior of the chain. The map of incision rates (Fig. 9b) based on the best-fit landscape shows a spatial distribution similar to ksn (Fig. 3). The higher incision rates that border the relict landscape, indicate that the chain interior will likely be captured by some rivers. This interpretation potentially rules out Quaternary climate change as one of the main drivers of landscape evolution in the study area. Indeed, the climate fluctuations have a time scale shorter than the tectonic forcing (106 y) and, especially in the Quaternary, they have been so rapid (40–100 ky cycle) that significant geomorphological adjustment may have not occurred (Whipple, 2001).

The values of local relief calculated on the synthetic landscape (Fig. 9a) have a range and distribution similar to that from the SRTM DEM (Fig. 2). In both data the highest values are ~400 m and occur in the northern, eastern and north-western flanks of the Iberian Chain and the Central System. The lowest values of local relief occur in the chain interior and the Madrid (High Tagus) and Guadiana basins. The results of SIGNUM modeling (Table 2) indicate that the simulated uplift rates (0.25–0.55 mm y−1) induce an average long-term incision rate of 0.22 mm y−1 over 3 My. The large-scale deformation and relatively low uplift rates producing the best fitting topography suggest that the relief generation of the Iberian Chain is probably due to mantle flow. Various studies indicate a slow velocity anomaly of seismic waves beneath central–eastern Iberia (Wortel and Spakman, 2000; Piromallo and Morelli, 2003; Li et al., 2008; Schmid et al., 2008; Boschi et al., 2009, 2010; Faccenna et al., 2010) A dynamic topography model predicted that this anomaly could have produced a positive topography (uplift) of ~500 m (Boschi et al., 2010; Faccenna et al., 2014). The calculation of the average incision rate for the last 100 ky (t = 100–0 ky) provided a value of 0.35 mm y−1 that is significantly higher than the average long-term incision rate. During this time interval, higher incision rates were noted in the central sector of the chain, ranging between 0.6 and 1 mm y−1 (Fig. 9b), similar to measured changes of incision rate obtained by Scotti et al. (2014). A recent study (Antón et al., 2014) analyzing geomorphic indices in the Duero Basin shows that, while tectonic uplift exerts a first-order control on drainage network evolution, drainage captures influence the further development of river networks by accelerated incision. The increase of modeled incision rates confirms that rivers have been incising the Iberian Chain flanks in response to a regional uplift and secondly to diachronous captures of the small intermontane basins. The long-term incision rates from stream profile modeling and from SIGNUM are very similar for the Tagus and Tajuña rivers. On the contrary the long-term incision rates from the Martín and Mijares profile modeling are slightly higher (Table 4). This difference is probably related to a model limitation: in SIGNUM, the uniform value of K (1.25 × 10−4 m0.4 y− 1) did not reproduce the bedrock variability present in the real Iberian Chain. Indeed, a finer reconstruction of K by river profile modeling and the comparison with geological data (see lithology bars in Figs. 10 and 11) allowed us to obtain values of incision rates more similar to those obtained from radiometric techniques. The incision coefficient of the stream power law calculated for all rock types along the modeled streams, falls in a range of values (5.8 × 10−5–5.0 × 10−4 m0.4 y−1) corresponding to the estimates of similar rock types reported in previous works (Table 5).

Table 5 Mean values of K obtained from this study or the literature. K

m

n

Lithology

References

7 × 10−3 (m0.2 y−1) 2.3 × 10−6 (m0.2 y−1) 2.3 × 10−6 (m0.2 y−1) 4.32 × 10−4 (m0.2 y−1) 5.7 × 10−4 (m0.2 y−1) 1.2 × 10−4 (m0.4 y−1) 2.42 × 10−5 (m0.2 y−1)

0.4 0.4 0.4 0.4 0.4 0.3 0.4

1 1 1 1 1 1 1

Mudstones Metasediments Granitoids Conglomerates, sandstones Sandstones, siltstones Limestones, dolostones, conglomerates, sandstone, mudstones Limestones, dolostones, conglomerates, sandstone, mudstones

Stock and Montgomery (1999) Stock and Montgomery (1999) Stock and Montgomery (1999) Kirby and Whipple (2001) Whipple et al. (2000) This study This study

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

From numerical experiments we defined a pattern of uplift and incision and quantified the rates of these processes. The results allowed us to reconstruct the Plio-Quaternary evolution of the Iberian Chain landscape. Starting from a topography characterized by large-size internally drained basins separated by low elevated ranges, an asymmetric doming uplifted the study area at ~ 3 Ma. This base level lowering induced the headward incision of the rivers, but most of the area remained internally drained. At first only the Tagus River started to capture the High Tagus (Madrid) Basin. Later around 1 Ma the Turia and Jalon rivers started to reach the intermontane basins. During this later piracy, the higher uplift rate of the Maestrazgo Mts. played an important role in preventing rivers from capturing the chain interior. 6. Conclusions This study has quantified the magnitude and trend of tectonic and erosion processes at the scale of the Iberian Chain, central–eastern Spain. The results of the landscape evolution modeling for the Iberian Chain and the best-fit parameters obtained from river profile modeling led to the following conclusions: 1. The Iberian Chain is an intraplate orogen that has been experiencing a recent regional uplift since ~3 Ma. This uplift was probably not uniform in the study area with a maximum in the north-eastern corner of the study area (Maestrazgo Mts.) and decreasing toward the southwestern sector. 2. The uplift rates that yielded the best-fit topography range from 0.25 to 0.55 mm y−1. 3. The rivers draining the flanks of the Iberian Chain have been responding to the regional tectonic event with an average longterm incision rate of 0.22 mm y−1, which is in accordance with the Late Pleistocene incision rates of ~ 0.6 mm y− 1 inferred from geochronological data by Scotti et al. (2014). 4. The resulting best-fit stream power parameters (m = 0.3, n = 1) are likely, and the bedrock erodibility K (5.8 × 10−5 − 5.0 × 10−4 m0.4 y−1) are similar to the values for similar lithologies reported in previous studies. 5. The capture of internally drained areas of the Iberian Chain was diachronous because of the different uplift rates. Acknowledgments We thank the Editor Takashi Oguchi, Martin Stokes and an anonymous reviewer for their constructive review, and Dr. Sean F. Gallen for improving the English text. References Agencia Estatal de Meteorología (AEMET-IM), 2011. Atlas climático de España y Portugal (Iberian Climate Atlas). Agencia Estatal de Meteorología, Ministerio de Medio Ambiente y Medio Rural y Marino, Instituto de Meteorologia de Portugal http:// www.aemet.es/es/serviciosclimaticos/datosclimatologicos/atlas_climatico. Alcalá, L., Alonso-Zarza, A.M., Álvarez Sierra, M.A., Azanza, B., Calvo, J.P., Cañaveras, J.C., van Dam, J.A., Garcés, M., Krijgsman, W., van der Meulen, A.J., Morale, J., PeláezCampomanes, P., Pérez Gonzalez, A., Sánchez Moral, S., Sancho, R., Sanz Rubio, E., 2000. El registro sedimentario y faunístico de las cuencas de Calatayud-Daroca y Teruel. Evolución paleoambiental y pal eoclimática durante el Neógeno. Rev. Soc. Geol. Esp. 13, 323–343. Alonso-Zarza, A., Calvo, J., 2000. Palustrine sedimentation in an episodically subsiding basin: the Miocene of the northern Teruel Graben (Spain). Palaeogeogr. Palaeoclimatol. Palaeoecol. 160, 1–21. Álvaro, M., Capote, R., Vegas, R., 1979. Un modelo de evolución geotectónica para la Cadena Celtibérica. Acta Geol. Hisp. 14, 172–177. Anadón, P., Moissenet, E., 1996. Neogene basins in the Eastern Iberian Range. In: Friend, P., Dabrio, C. (Eds.), Tertiary Basins of Spain, the Stratigraphic Record of Crustal Kinematics. Cambridge University Press, Cambridge, pp. 68–76. Anadón, P., Moissenet, P., Simón, J., 1990. The Neogene grabens of the Eastern Iberian Chain (Eastern Spain). Paleontol. Evol. Mem. Espec. 2, 97–130. Antón, L., De Vicente, G., Muñoz-Martín, A., Stokes, M., 2014. Using river long profiles and geomorphic indices to evaluate the geomorphological signature of continental

65

scale drainage capture, Duero basin (NW Iberia). Geomorphology 206, 250–261. http://dx.doi.org/10.1016/j.geomorph.2013.09.028. Armenteros, I., Dabrio, C.J., Guisado, R., Sànchez de Vega, A., 1989. Megasecuencias sedimentarias del terciario del borde oriental de la Cuenca de Almazán (Soria– Zaragoza). Stud. Geol. Salmanticensia 5, 107–127 (Special Volume). Babault, J., Loget, N., Van Den Driessche, J., Castelltort, S., Bonnet, S., Davy, P., 2006. Did the Ebro basin connect to the Mediterranean before the Messinian salinity crisis? Geomorphology 81, 155–165. Birot, P., 1959. Esquisse morphologique des Monts Celtibériques orientaux. Bull. Com. Trav. Hist. Sci. Sect. Géog. 1 (72), 101–130. Bishop, P., 1995. Drainage rearrangement by river capture, beheading and diversion. Prog. Phys. Geogr. 19, 449–473. Boschi, L., Fry, B., Ekström, G., Giardini, D., 2009. The European upper mantle as seen by surface waves. Surv. Geophys. 30, 463–501. http://dx.doi.org/10.1007/s10712-0099066-2. Boschi, L., Faccenna, C., Becker, T.W., 2010. Mantle structure and dynamic topography in the Mediterranean Basin. Geophys. Res. Lett. 37, L20303. http://dx.doi.org/10.1029/ 2010GL045001. Braun, J., Sambridge, M., 1997. Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretization. Basin Res. 9, 27–52. Burov, E., Toussaint, G., 2007. Surface processes and tectonics: forcing of continental subduction and deep processes. Glob. Planet. Chang. Spec. Vol. Topogr. Eur. 58, 141–164. Calvo, J.P., 2004. Rasgos comunes de las Cuencas cenozoicas. In: Vera, J.A. (Ed.), Geología de España. SGE, IGME, Madrid, pp. 584–586. Calvo Hernández, J.M., 1993. Cinemática de las fallas discontinuas en el sector central de la Cordillera Ibérica (PhD Thesis). University of Zaragoza, Spain. Calvo, J.P., Daams, R., Morales, J., López-Martínez, N., Agustí, J., Anadón, P., Armenteros, I., Cabrera, L., Civis, J., Corrochano, A., Díaz-Molina, M., Elizaga, E., Hoyos, M., MartínSuárez, E., Martínez, J., Moissenet, E., Muñoz, A., Pérez-García, A., Pérez-González, A., Portero, J.M., Robles, F., Santisteban, C., Torres, T., Van der Meulen, A.J., Vera, J.A., Mein, P., 1993. Up-to-date Spanish continental Neogene synthesis and paleoclimatic interpretation. Rev. Soc. Geol. Esp. 6 (3–4), 29–40. Capolongo, D., Giachetta, E., Refice, A., 2011. Numerical framework for geomorphological experiments. Geogr. Fis. Din. Quat. 34, 75–80. http://dx.doi.org/10.4461/GFDQ.2011. 34.8. Carretier, S., Lucazeau, F., 2005. How does alluvial sedimentation at range fronts modify the erosional dynamics of mountain catchments? Basin Res. 17 (3), 361–381. Casas-Sainz, A.M., Cortés-Gracia, A.L., 2002. Cenozoic landscape development within the Central Iberian Chain, Spain. Geomorphology 44, 19–46. Casas-Sainz, A., De Vicente, G., 2009. On the tectonic origin of Iberian topography. Tectonophysics 474, 214–235. Casas-Sainz, A.M., Faccenna, C., 2001. Tertiary compressional deformation of the Iberian plate. Terra Nova 13, 281–288. Casas-Sainz, A.M., Cortés, A., Gapais, D., Nalpas, T., Román, T., 1998. Modelización analógica de estructuras asociadas a compresión oblicua y transpresión. Ejemplos del NE peninsular. Rev. Soc. Geol. Esp. 11 (3–4), 137–150. Clauzon, G., Suc, J.P., Gautier, F., Berger, A., Loutre, M.F., 1996. Alternate interpretation of the Messinian salinity crisis: controversy resolved? Geology 24, 363–366. Cloetingh, S., Ziegler, P., Bogaard, P., Andriessen, P., Artemieva, I., Bada, G., Balen van, R., Beekman, F., Ben-Avraham, Z., Brun, J.P., Bunge, H.P., Burov, E., Crabonell, R., Facenna, C., Friedrich, A., Gallart, C., Green, A., Heidbach, O., Jones, A., Matenco, L., Mosar, J., Oncken, O., Pascal, C., Peters, G., Sliaupa, S., Soesoo, A., Spakman, W., Stephenson, R., Thybo, H., Torsvik, T., Vicente de, G., Wenzel, F., Wortel, M., 2007. TOPO-EUROPE: the geoscience of coupled deep Earth-surface processes. Glob. Planet. Chang. 58, 1–118. Cortés-Gracia, A.L., Casas-Sainz, A.M., 1996. Deformación alpina de zócalo y cobertera en el borde norte de la Cordillera Ibérica. Rev. Soc. Geol. Esp. 9 (1–2), 51–66. Cortés-Gracia, A.L., Casas-Sainz, A.M., 2000. ¿Tiene el sistema de fosas de Teruel origen extensional? Rev. Soc. Geol. Esp. 13 (3–4), 445–470. Courant, R., Friedrichs, K., Lewy, H., 1967. On the partial difference equations of mathematical physics. IBM J. 215–234. De Bruijne, C., Andriessen, P., 2002. Far field effects of alpine plate tectonismin the Iberian microplate recorded by fault-related denudation in the Spanish Central System. Tectonophysics 349, 161–184. De Vicente, G., Vegas, R., Muñoz Martín, A., Silva, P.G., Andriessen, P., Cloetingh, S., González Casado, J.M., VanWees, J.D., Álvarez, J., Carbó, A., Olaiz, A., 2007. Cenozoic thick-skinned deformation and topography evolution of the Spanish Central System. Glob. Planet. Chang. 58, 335–381. Densmore, A., Ellis, M., Anderson, R., 1998. Landsliding and the evolution of normal-faultbounded mountains. J. Geophys. Res. 103 (B7), 15203–15219. Dietrich, W.E., Reiss, R., Hsu, M., Montgomery, D.R., 1995. A process-based model for colluvial soil depth and shallow landsliding using digital elevation data. Hydrol. Process. 9, 383–400. Faccenna, C., Becker, T.W., Lallemand, S., Lagabrielle, Y., Funiciello, F., Piromallo, C., 2010. Subduction-triggered magmatic pulses: a new class of plumes? Earth Planet. Sci. Lett. 299, 54–68. http://dx.doi.org/10.1016/j.epsl.2010.08.012. Faccenna, C., Becker, T.W., Auer, L., Billi, A., Boschi, L., Brun, J.P., Capitanio, F.A., Funiciello, F., Horvàth, F., Jolivet, L., Piromallo, C., Royden, L., Rossetti, F., Serpelloni, E., 2014. Mantle dynamics in the Mediterranean. Rev. Geophys. 52, 283–332. http://dx.doi. org/10.1002/2013RG000444. Fernández, M.H., Álvarez Sierra, M.A., Peláez-Campomanes, P., 2007. Bioclimatic analysis of rodent palaeofaunas reveals severe climatic changes in Southwestern Europe during the Plio-Pleistocene. Palaeogeogr. Palaeoclimatol. Palaeoecol. 251 (3–4), 500–526. http://dx.doi.org/10.1016/j.palaeo.2007.04.015. Ferreiro, E., Ruiz, V., Lendínez, A., Lago, M., Meléndez, A., Pardo, G., Ardevol, L., Villena, J., Hernández, A., Alvaro, M., Gómez, J.J., Carls, P., 1991. Mapa y memoria explicativa

66

E. Giachetta et al. / Geomorphology 246 (2015) 48–67

de la Hoja 40 (Daroca) del Mapa Geológico Nacional a escala 1:200.000, Instituto Tecnológico y Geominero de España, 239 p. Flint, J.J., 1974. Stream gradient as a function of order, magnitude, and discharge. Water Resour. Res. 10, 969–973. Gasparini, N.M., Whipple, K.X., Bras, R.L., 2007. Predictions of steady-state and transient landscape morphology using sediment-flux dependent river incision models. J. Geophys. Res. 112, F03S09. http://dx.doi.org/10.1029/2006JF000567. Giachetta, E., Refice, A., Capolongo, D., Gasparini, N.M., Pazzaglia, F.J., 2014. Orogen-scale drainage network evolution and response to erodibility changes: insights from numerical experiments. Earth Surf. Process. Landf. 39, 1259–1268. http://dx.doi.org/ 10.1002/esp.3579. Gracia, F.J., Gutiérrez, F., Gutiérrez, M., 2003. The Jiloca karst polje-tectonic graben (Iberian Range, NE Spain). Geomorphology 52, 215–231. Gracia-Prieto, F.J., Gutiérrez Elorza, M., Leránoz Istúriz, B., 1988. Las superficies de erosión neógenas en el sector central de la Cordillera Ibérica. Rev. Soc. Geol. Esp. 1 (1–2), 135–142. Guimerà, J., Alvaro, M., 1990. Structure et évolution de la compression alpine dans la Chaîne Ibérique et al. Chaîne Côtière Catalane. Bull. Soc. Geol. Fr. 8 (6), 339–348. Guimerà, J., Más, R., Alonso, A., 2004. Intraplate deformation in the NW Iberian Chain: Mesozoic extension and contractional inversion. J. Geol. Soc. Lond. 16, 291–303. Gutiérrez Elorza, M., Gracia, F.J., 1997. Environmental interpretation and evolution of the Tertiary erosion surfaces in the Iberian Range (Spain). In: Widdowson, M. (Ed.), Palaeosurfaces: Recognition, Reconstruction and Palaeoenvironmental Interpretation. Geological Society, London, Special Publication 120, pp. 147–158. Gutiérrez, F., Gracia, J., Gutiérrez, M., 1996. Consideraciones sobre el final del relleno endorreico de las fosas de Calatayud y Teruel y su paso al exorreísmo. Implicaciones morfo-estratigráficas y estructurales. IV Reunión Nacional de Geomorfología, S.E.G., O Castro (A Coruña), Spain, pp. 23–43. Gutiérrez, F., Gutiérrez, M., Gracia, F.J., McCalpin, J.P., Lucha, P., Guerrero, J., 2008. PlioQuaternary extensional seismotectonics and drainage network development in the central sector of the Iberian Range (NE Spain). Geomorphology 102 (1), 21–42. Gutiérrez, F., Masana, E., González, A., Guerrero, J., Lucha, P., McCalpin, J.P., 2009. Late Quaternary paleoseismic evidence on the Munébrega Half-graben fault (Iberian Range, Spain). Int. J. Earth Sci. 98, 1691–1703. http://dx.doi.org/10.1007/s00531008-0319-y. Gutiérrez, F., Gracia, F.J., Gutiérrez, M., Lucha, P., Guerrero, J., Carbonel, D., Galve, J.P., 2012. A review on Quaternary tectonic and nontectonic faults in the central sector of the Iberian Chain, NE Spain. Fallas cuaternarias tectónicas y gravitacionales en el sector central de la Cordillera Ibérica, NE de España. J. Iber. Geol. 38 (1), 145–160. Hack, J.T., 1957. Studies of longitudinal profiles in Virginia and Maryland. U. S. Geol. Surv. Prof. Pap. 294 (B), 45–97. Hack, J.T., 1973. Stream-profile analysis and stream-gradient index. J. Res. U. S. Geol. Surv. 1, 421–429. Hancock, G.R., Coulthard, T.J., Martinez, C., Kalma, J.D., 2011. An evaluation of landscape evolution models to simulate decadal and centennial scale soil erosion in grassland catchments. J. Hydrol. 398 (3), 171–183. Hovius, N., 2000. Macro-scale process systems of mountain belt erosion. In: Summerfield, M.A. (Ed.), Geomorphology and Global Tectonics. John Wiley and Sons, London, pp. 77–105. Howard, A.D., 1980. Thresholds in river regime. In: Coates, D.R., Vitek, J.D. (Eds.), The Concept of Geomorphic Thresholds. Allen and Unwin, Boston, pp. 227–258 (Ch. 11). Howard, A.D., 1994. A detachment-limited model of drainage basin evolution. Water Resour. Res. 30 (7), 2261–2285. Howard, A.D., 1997. Badland morphology and evolution: interpretation using a simulation model. Earth Surf. Process. Landf. 22, 211–227. Howard, A.D., Kerby, G., 1983. Channel changes in badlands. Geol. Soc. Am. Bull. 94, 739–752. Howard, A.D., Dietrich, W.E., Seidl, M.A., 1994. Modeling fluvial erosion on regional to continental scales. J. Geophys. Res. 99, 13971–13987. http://dx.doi.org/10.1029/ 94JB00744. IGME (Ed.), 1994. Mapa Geológico de la Península Ibérica, Baleares y Canarias, 1: 1.000.000. Madrid. Instituto Tecnológico y Geominero de España. Jiménez-Moreno, G., Fauquette, S., Suc, J.P., 2010. Miocene to Pliocene vegetation reconstruction and climate estimates in the Iberian Peninsula from pollen data. Rev. Palaeobot. Palynol. 162 (3), 403–415. http://dx.doi.org/10.1016/j.revpalbo.2009.08. 001. Kirby, E., Whipple, K., 2001. Quantifying differential rock-uplift rates via stream profile analysis. Geology 29, 415–418. Lafuente, P., Arlegui, L.E., Casado, I., Ezquerro, L., Liesa, C.L., Pueyo, Ó., Simón, J.L., 2011a. Geometría y cinemática de la zona de relevo entre las fallas Neogéno-Cuaternarias del Concud y Teruel (Cordillera Ibérica). Rev. Soc. Geol. Esp. 24 (1–2), 117–133. Lafuente, P., Arlegui, L.E., Liesa, C.L., Simón, J.L., 2011b. Paleo-seismological analysis of an intraplate extensional structure: the Concud fault (Iberian Chain, Eastern Spain). Int. J. Earth Sci. 100, 1713–1732. http://dx.doi.org/10.1007/s00531-010-0542-1. Li, C., van der Hilst, R.D., Engdahl, E.R., Burdick, S., 2008. A new global model for P-wave speed variations in Earth's mantle. Geochem. Geophys. Geosyst. 9, Q05018. http://dx.doi.org/10.1029/2007GC001806 Q05018. Loget, N., Van Den Driessche, J., 2009. Wave train model for knickpoint migration. Geomorphology 106, 376–382. López-Martínez, N., Agustí, J., Cabrera, L., Calvo, J., Civis, J., Corrochano, A., Daams, R., Díaz, M., Elizaga, E., Hoyos, M., Martínez, J., Morales, J., Portero, J.M., Robles, F., Santisteban, C., Torres, T., 1987. Approach to the Spanish continental Neogene synthesis and paleoclimatic interpretation. Ann. Inst. Geol. Public. Hung. 70, 383–391. Martín-Serrano, A., 1991. La definición y el encajamiento de la red fluvial actual sobre el Macizo Hespérico en elmarco de su geodinámica alpina. Rev. Soc. Geol. Esp. 4, 337–351.

Mather, A.E., 2000. Adjustment of a drainage network to capture induced base-level change: an example from the Sorbas basin, SE Spain. Geomorphology 34, 271–289. Muñoz, J.A., 1992. Evolution of a continental collision belt: ECORS-Pyrenees crustal balanced cross-section. In: McClay, K.R. (Ed.), Thrust tectonics. Springer, Netherlands, pp. 235–246. Muñoz-Martín, A., De Vicente, G., 1998. Origen y relación entre las deformaciones y esfuerzos alpinos de la zona centro-oriental de la Península Ibérica. Rev. Soc. Geol. Esp. 11, 57–70. Pelletier, J.D., 2007. Numerical modeling of the Cenozoic fluvial evolution of the southern Sierra Nevada, California. Earth Planet. Sci. Lett. 259, 85–96. http://dx.doi.org/10. 1016/j.epsl.2007.04.030. Peña, J.L., Gutiérrez, M., Ibáñez, M.J., Lozano, M.V., Rodríguez, J., Sánchez, M., Simón, J.L., Soriano, A., Yetano, M., 1984. Geomorfología de la Provincia de Teruel. Instituto de Estudios Turolenses, Excma. Dip., Provincial de Teruel (Spain), p. 149. Perea, H., Masana, E., Santanach, P., 2012. An active zone characterized by slow normal faults, the northwestern margin of the València trough (NE Iberia): a review. J. Iber. Geol. 38 (1), 31–52. Piromallo, C., Morelli, A., 2003. P-wave tomography of the mantle under the Alpine‐Mediterranean area. J. Geophys. Res. 108 (B2). http://dx.doi.org/10.1029/2002JB001757. Refice, A., Giachetta, E., Capolongo, D., 2012. SIGNUM: a Matlab, TIN-based landscape evolution model. Comput. Geosci. 45, 293–303. http://dx.doi.org/10.1016/j.cageo.2011. 11.013. Roberts, G.G., White, N., 2010. Estimating uplift rate histories from river profiles using African examples. J. Geophys. Res. 115, B02406. http://dx.doi.org/10.1029/ 2009JB006692. Roberts, G.G., Paul, J.D., White, N., Winterbourne, J., 2012a. Temporal and spatial evolution of dynamic support from river profiles: a framework for Madagascar. Geochem. Geophys. Geosyst. 13, Q04004. http://dx.doi.org/10.1029/2012GC004040. Roberts, G.G., White, N.J., Martin-Brandis, G.L., Crosby, A.G., 2012b. An uplift history of the Colorado Plateau and its surroundings from inverse modeling of longitudinal river profiles. Tectonics 31, TC4022. http://dx.doi.org/10.1029/2012TC003107. Roca, E., Guimerà, J., 1992. The Neogene structure of the eastern Iberian margin: structural constraints on the crustal evolution of the Valencia trough (western Mediterranenan). Tectonophysics 203, 203–218. Rodríguez-Pascua, M., De Vicente, G., 1998. Análisis de paleoesfuerzos en cantos de depósitos conglomeráticos terciarios de la cuenca de Zaorejas (rama castellana de la Cordillera Ibérica). Rev. Soc. Geol. Esp. 11, 169–180. Roe, G.H., Stolar, D.R., Willett, S.D., 2006. Response of a steady-state critical wedge orogen to changes in climate and tectonic forcing. In: Willett, S.D., Hovius, N., Brandon, M.T., Fisher, D.M. (Eds.), Tectonics, Climate, and Landscape Evolution. Geological Society of America Special Paper 398, Penrose Conference Series. Geological Society of America, Boulder, CO, pp. 227–239. Roering, J.J., Kirchner, J.W., Dietrich, W.E., 1999. Evidence for non-linear, diffusive sediment transport on hillslopes and implications for landscape morphology. Water Resour. Res. 35 (3), 853–870. Rubio, J.C., Simon, J.L., 2007. Tectonic subsidence v. erosional lowering in a controversial intramontane depression: the Jiloca basin (Iberian Chain, Spain). Geol. Mag. 144, 127–141. Schmid, C., van der Lee, S., VanDecar, J.C., Engdahl, E.R., Giardini, D., 2008. Three‐ dimensional S velocity of the mantle in the Africa–Eurasia plate boundary region from phase arrival times and regional waveforms. J. Geophys. Res. 113, B03306. http://dx.doi.org/10.1029/2005JB004193. Scotti, V.N., Molin, P., Faccenna, C., Soligo, M., Casas-Sainz, A., 2014. The influence of surface and tectonic processes on landscape evolution of the Iberian Chain (Spain): Quantitative geomorphological analysis and geochronology. Geomorphology 206, 37–57. http://dx.doi.org/10.1016/j.geomorph.2013.09.017. Simón, J.L., 1984. Compresión y distensión alpinas en la Cadena Ibérica oriental. Instituto de Estudios Turolenses, Teruel (Spain). Simón, J.L., 1989. Late Cenozoic stress field and fracturing in the Iberian Chain and Ebro Basin (Spain). J. Struct. Geol. 11, 285–294. Simón, J.L., Arlegui, L.E., Lafuente, P., Liesa, C.L., 2012. Active extensional faults in the central–eastern Iberian Chain, Spain. J. Iber. Geol. 38, 127–144. Simón, J.L., Pérez-Cueva, A.J., Calvo-Cases, A., 2013. Tectonic beheading of fluvial valleys in the Maestrat grabens (eastern Spain): insights into slip rates of Pleistocene extensional faults. Tectonophysics 593, 73–84. http://dx.doi.org/10.1016/j.tecto.2013.02.026. Sklar, L.S., Dietrich, W.E., 1998. River longitudinal profiles and bedrock incision models: stream power and the influence of sediment supply. In: Tinkler, K.J., Wohl, E.E. (Eds.), Rivers over rock: Fluvial processes in bedrock channels. American Geophysical Union, Washington, D.C., pp. 237–260. Solé Sabarís, L., 1979. La Meseta. In: Teràn, D. (Ed.), Geografía de España. Ariel, Madrid, pp. 42–62. Stock, J.D., Montgomery, D.R., 1999. Geologic constraints on bedrock river incision using the stream power law. J. Geophys. Res. 104, 4983–4993. Stokes, M., Mather, A.E., Harvey, A.M., 2002. Quantification of river-capture-induced baselevel changes and landscape development, Sorbas basin, SE Spain. In: Jones, S.J., Frostick, L.E. (Eds.), Sediment Flux to Basins: Causes, Controls and Consequences. Special Publications 191. Geological Society, London, pp. 23–35. Stolar, D.B., Willett, S.D., Roe, G.H., 2006. Climatic and tectonic forcing of a critical orogen. In: Willett, S.D., Hovius, N., Brandon, M.T., Fisher, D.M. (Eds.), Tectonics, Climate, and Landscape Evolution: Geological Society of America Special Paper 398, Penrose Conference Series. Geological Society of America, Boulder, CO, pp. 241–250. Teixell, A., 1998. Crustal structure and orogenic material budget in the west central Pyrenees. Tectonics 17 (3), 395–406. Ter Voorde, M., De Bruijne, C., Cloetingh, S., Andriessen, P., 2004. Thermal consequences of thrust faulting: simultaneous versus successive fault activation and exhumation Earth Planet. Sci. Lett. 223, 397–415.

E. Giachetta et al. / Geomorphology 246 (2015) 48–67 Tucker, G.E., Bras, R.L., 2000. A stochastic approach to modeling the role of rainfall variability in drainage basin evolution. Water Resour. Res. 36 (7), 1953–1964. Tucker, G.E., Slingerland, R.L., 1994. Erosional dynamics, flexural isostasy, and long-lived escarpments: a numerical modeling study. J. Geophys. Res. 99 (B6), 12229–12243. Tucker, G.E., Whipple, K.X., 2002. Topographic outcomes predicted by stream erosion models: sensitivity analysis and intermodel comparison. J. Geophys. Res. 107, 2179. http://dx.doi.org/10.1029/2001JB000162. Tucker, G.E., Lancaster, S.T., Gasparini, N.M., Bras, R.L., 2001. The Channel-Hillslope Integrated Landscape Development (CHILD) Model. In: Harmon, R.S., Doe III, W.W. (Eds.), Landscape Erosion and Evolution Modeling. Kluwer Academic/Plenum Publishers, New York, NY, pp. 349–388. van Dam, J.A., 2006. Geographic and temporal patterns in the late Neogene (12–3 Ma) aridification of Europe. The use of small mammals as paleoprecipitation proxies. Palaeogeogr. Palaeoclimatol. Palaeoecol. 238, 190–218. van Dam, J., Sanz Rubio, E., 2003. Late Miocene and Pliocene small mammals from the Calatayud Basin (Central Spain). Coloquios Paleontol. 1, 115–126. van Dam, J.A., Weltje, G.J., 1999. Reconstruction of the late Miocene climate of Spain using rodent paleocommunity successions: an application of end-member modelling. Palaeogeogr. Palaeoclimatol. Palaeoecol. 151, 267–305. Vegas, R., 1992. The Valencia trough and the origin of the western Mediterranean basins. Tectonophysics 203, 249–261. Vergés, J., Fernández, M., Martínez, A., 2002. The Pyrenean orogen: pre-, syn-, and postcollisional evolution. In: Rosenbaum, G., Gordon, L. (Eds.), Reconstruction of the evolution of the Alpine-Himalayan orogen. Journal of the Virtual Explorer, Electronic Edition vol. 8, pp. 55–74. http://dx.doi.org/10.3809/jvirtex.2002.00058. Waelbroeck, C., Labeyrie, L., Michel, E., Duplessy, J.C., McManus, J.F., Lambeck, K., Balbon, E., Labracherie, M., 2002. Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records. Quat. Sci. Rev. 21, 295–305. Wegmann, K.W., Zurek, B.D., Regalla, C.A., Bilardello, D., Wollenberg, J.L., Kopczynski, S.E., Ziemann, J.M., Haight, S.L., Apgar, J.D., Zhao, C., Pazzaglia, F.J., 2007. Position of the

67

Snake River watershed divide as an indicator of geodynamic processes in the greater Yellowstone region, western North America. Geosphere 3, 272–281. http://dx.doi. org/10.1130/GES00083.1. Whipple, K.X., 2001. Fluvial landscape response time: how plausible is steady state denudation? Am. J. Sci. 301, 313–325. http://dx.doi.org/10.2475/ajs.301.4-5.313. Whipple, K.X., 2004. Bedrock rivers and the geomorphology of active orogens. Annu. Rev. Earth Planet. Sci. 32, 151–185. Whipple, K.X., Tucker, G.E., 1999. Dynamics of the stream-power river incision model: implications for height limits of mountain ranges, landscape response time-scales, and research needs. J. Geophys. Res. 104 (B8), 17661–17674. Whipple, K.X., Tucker, G.E., 2002. Implications of sediment-flux dependent river incision models for landscape evolution. J. Geophys. Res. 107, B2. http://dx.doi.org/10.1029/ 2000JB000044. Whipple, K.X., Snyder, N., Dollenmayer, K., 2000. Rates and processes of bedrock incision by the Upper Ukak River since the 1912 Novarupta ash flow in the Valley of Ten Thousand Smokes, Alaska. Geology 28, 835–838. Wilcox, R.E., Harding, T.P., Seely, D.R., 1973. Basic wrench tectonics. Assoc. Pet. Geol. Bull. 57 (1), 74–96. Willett, S.D., 1999. Orogeny and orography: the effects of erosion on the structure of mountain belts. J. Geophys. Res. 104, 28957–28981. http://dx.doi.org/10.1029/ 1999JB900248. Willett, S.D., Slingerland, R., Hovius, N., 2001. Uplift, shortening, and steady state topography in active mountain belts. Am. J. Sci. 301, 455–485. Wortel, M.J.R., Spakman, W., 2000. Subduction and slab detachment in the Mediterranean– Carpathian region. Science 290, 1910–1917. http://dx.doi.org/10.1126/science.290. 5498.1910. Wright, V.P., Alonso-Zarza, A.M., Sanz, M.E., Calvo, J.P., 1997. Diagenesis of Late Miocene micritic lacustrine carbonates, Madrid Basin, Spain. Sediment. Geol. 114, 81295.