The main lithologies which outcrop in the Lozoya watershed are the Paleozoic rocks of the Iberian Massif (Figure 1b). The oldest rocks consist of metasediments, paragneiss, quartzites and marbles of Precambrian age, situated to the NE (Arenas et al., 1991; Bellido et al., 1991a, 1991b). In the eastern segment, schists, quartzites and shales of Ordovician and Silurian age outcrop (Capote & Fernández, 1975), whereas the south is characterized by pre-Permian granites and adamellites. All these lithologies were affected by penetrative schistosity and folding during the Variscan orogeny (Arenas et al., 1991; Bellido et al., 1991a, 1991b). During the Permian, dikes of porphyry, lamprophyre, quartz, microdiorite and aplite were emplaced in the SE segment. Mesozoic deposits are located in contact with the Madrid Cenozoic Basin and in the Upper Lozoya Valley (pop-down depression) (NE segment). The base of these Cretaceous deposits is composed of fluvial sediments, the so-called Facies Utrillas (clays, gravels, lutites, carbonates, sandstones and dolomites), overlain by sands, lutites, carbonates, sandstones and dolomites deposited during the Turonian-Santonian marine transgression (Aznar et al., 1995). Above the Cretaceous sequence, there are lake deposits and gypsum. At the beginning of the Alpine movements (Paleocene), conglomerates with Mesozoic calcareous boulders were deposited. This episode corresponds to a period of thrusting before the arkosic filling of the Cenozoic basins (Capote et al., 1987). In the Pinilla del Valle area and to the west, Miocene detritic alluvial sediments are present. The Miocene is characterized by thrusting generating NE-SW structures and folding of the Mesozoic and Cenozoic deposits in the Iberian Mountain range (Vicente et al., 2007). Finally, several impulses of vertical uplifting as a result of regional extensional stresses have been defined in different parts of the SCS during the Pliocene and up to the present day (De Bruijne & Andrienssen, 2002; Vicente et al., 2007). However, in this sector, there is no evidence for neotectonic activity or historic seismicity (Baena-Perez et al., 1998; http://www.ign.es/ign/layoutIn/sismoTerremotosEspana.do). Quaternary landforms in the Lozoya watershed are mainly erosive: deposits located and associated with fluvial and glacial processes, colluviums, peat formations, solifluction formations and karstic environments have been described in several studies (Pedraza, 1994; Torres et al., 1995; Torres et al., 2005; Arsuaga et al., 2006; Pedraza & Carrasco 2005; Pérez-González et al., 2010; Palacios et al., 2011).

The main landforms identified in the Lozoya watershed can be grouped as erosive surfaces, polygenetic residual features (inselbergs, glacis), fluvial morphologies (e.g. strath terraces, rock terraces or alluvial fans), glacial landforms (e.g. cirques, moraines), karstic features (e.g. caves, dolines, karren, karstic infillings), and gravity forms (Karampaglidis et al., 2014a).

Our work is focused on the Quaternary fluvial landforms (Figures 2, 3 & 4), which form an erosive terrace staircase arrangement (Karampaglidis et al., 2014a). According to the characteristics and distribution of these glacis and terraces, together with the morphostructural and bedrock characteristics of the valley, we have divided the Lozoya watershed into four sectors: the Pinilla del Valle sector, the Buitrago de Lozoya sector, the Atazar sector and the Cerro de la Oliva sector (Figure 3, 4 & 5; Karampaglidis et al., 2011).

The first sector of the watershed (Pinilla del Valle) is a NE-SW pop-down tectonic depression (Warburton & Alvarez, 1989; Vicente et al., 2007), which determines the slope distribution and the main direction of the Lozoya River in this area (Figure 5). The geology here is characterized by soft Mesozoic and Cenozoic sedimentary rocks. In this sector, the fluvial style of the Lozoya River changes from a shallow gravel braided facies model to a meandering gravel-bed facies model. The river load is characterized by clastic rocks and polygenic composition (metamorphic and igneous rocks: Figure 1). In this area, the earliest terraces and glacis are erosive, degraded and with reduced dimensions, developing on the Paleozoic metamorphic rocks between T2 (+190-195 m) and T15 (+62-64 m). Below them, rock terraces and glacis between T16 (+50-55 m) and T20 (+17-20 m) are preserved on Mesozoic and Cenozoic rocks. This sector is one of the few areas where fluvial deposits are located, composed of cobbles, coarse-medium gravel, sands and muds corresponding to strath terraces, alluvial fans and fluvial deposits preserved in a karstic system (Pinilla del Valle caves) (Figure 3). In this karstic system, we have observed three episodes of fluvial aggradation situated at +23-25 m (T19), +11-14 m (T21), and +7 (T22) (Karampaglidis et al., 2014a). The second episode (T21, +11-14 m) was dated by thermoluminescence (TL) on quartz grains at 140.4 ± 11.3 ka B.P., while the deposits at +7 m (T22) have an age of 59.4 ± 4.7 ka B.P. (Thermoluminescence Laboratory, Universidad Autónoma de Madrid; Pérez-González et al., 2010). Furthermore, two other levels of strath terraces (T23, +3-5 m & T24, +1-2 m) were identified in this sector. The deposits of T23 were dated at 17,950 ± 1.1 ka B.P. using accelerator mass spectrometry (AMS) in organic sediments (Karampaglidis et al., 2014a). The stream load energy may increase as a result of climate change and/or tectonic activity (Merrits et al., 1994; Pazzaglia et al., 1998). In this sector, local fluvial deposit formation of the last two terrace levels (T22, +7 m & T23 +3-5 m) and valley widening coincide with some of the well-documented periods of local glaciation (19 ka) or cold periods (~60 ka) (Carrasco et al., 2016).

In the second sector (Buitrago de Lozoya), the course of the Lozoya River can be classified as an incised meander bedrock river (sinuosity of 2.71; Karampaglidis et al., 2011), entrenched from T18 (+30-35 m), showing a NE-SW general direction following the NE-SW trends of the SCS pop-up structure (Vicente et al., 2007). In this sector, most of the preserved terraces are erosive and the oldest terraces and glacis are poorly preserved. The surface of the terraces and glacis is characterized by rock-weathering material, such as angular blocks and clasts, with a fine soil covering of 3-4 cm thickness. Occasionally, low terraces (T23, +3-4 m, & T24, +1-2 m) contain 40-50 cm of rounded oligomictic boulders and cobbles, sands and silts (Karampaglidis et al., 2014a). Towards the north, river captures were identified associated with the evolution of the Madarquillo River (Figure 1). Part of this river is the Hontanar stream which presents an unexpected mature-stage meander higher than the local base level of the Valle stream, a short local linear stream situated on the structural contact between gneiss and metasediments (Paleozoic) (Figure 6). The second capture is related to the parallel valley of the Retamar stream (Figure 6). This is a hanging fluvial valley captured from the Valle stream. These captures are related to greater incision, produced when the Valle stream incised the structural contact and/or by fault reactivation, and probably occurred after T13, +80-85 m (Early Pleistocene), since this is the last common rock terrace (Karampaglidis et al., 2011; Karampaglidis, 2015).

In the Atazar sector, the morphology of the Lozoya River is controlled by E-W and N-S fractures (Figure 1 & 2a), and this explains its N-S direction of flow in the southern part of this sector. The valley describes incised meanders entrenched from T13 (+80-85 m), displaying a sinuosity of 2.68 (Figure 5; Karampaglidis et al., 2011). In this sector, rock terraces and glacis sequences are located on the Paleozoic rocks and on Neogene deposits (Figure 1 & 3). The youngest Neogene deposits appear in this sector, forming a surface perched about +220 m above the present river channel (Figure 1, 3 & 4). Below this level, fluvial rock terraces are preserved on the Pliocene deposits (from T1, +200-205 m to T7, +142-144 m), and on Paleozoic and Permian rocks, where the granites are often affected by weathering processes (Karampaglidis et al., 2014a). In the lower part of this sector, a fluvial capture was recognized (Figure 6), changing the direction of the Lozoya River. During the formation of the first terraces between T1 (+200-205 m) and T4 (+170 m), the Lozoya River drained towards the south through the El Berrueco Pass (Figure 2 & 6). In the subsequent phase, the river and its tributaries were captured, changing the flow direction to the east (Figure 6), where the lower sector develops (Karampaglidis et al., 2011). This capture was recognized during the morphostatistical reconstruction of the Lozoya River base levels and from geomorphological evidence. Large planation surfaces are suspended above the current base level and the present local tributaries cannot explain their formation (Karampaglidis et al., 2011). This capture probably occurred by local fault reactivation from southern border activity, Early Pleistocene climatic fluctuations and/or headward erosion, which is related with higher incision rates and was produced during the Lozoya River’s incision into the regional structural contact between gneisses, granites and black slates (Figure 1).

The lowest sector, or Cerro de la Oliva sector, is characterized by higher incision (Figure 5), forming canyons developed from T5 (+160-165 m), with the highest sinuosity (3.36; Karampaglidis et al., 2011), where fluvial meandering and meander cutoffs have been mapped. In this area, the rock terraces are located on Paleozoic black shales and Neogene deposits. In addition, strath terraces were identified at T16 (+50-55 m), T18 (+30-35 m), T19 (+25 m), T22 (+6-8 m) and T23 (+3-5 m) levels (Figure 4; Karampaglidis et al., 2014a). Fluvial deposits are also preserved in karstic systems (Reguerillo Cave), hanging +140 m above the actual base level of the modern Lozoya River (Figure 2). This deposit contains polymict (igneous and metamorphic rocks) sands, silts and clays, deposited by the Lozoya River (Torres et al., 2005). Those authors have undertaken various measurements of residual magnetism in these deposits, characterized by a sequence of five samples of normal polarity. This magnetostratigraphic sequence and the position of this fluvial level in the regional sequence (Figure 7) suggest an Olduvai subchron age (1.77 ka to 1.95 ka) for these fluvial deposits. This karstic system also contains fluvial deposits in a level at +45 m, below a speleothem dated by ESR at 981 ± 76 ka B.P. (Torres et al., 2005). This age could suggest that this terrace is older, although the date is inconsistent with the rest of the chronological trend observed for the sequence in the MCB (Santonja & Pérez-González, 2001; Ordoñez et al., 1990; Ortiz et al., 2009, Silva et al., 2017), which placed terraces at +42-46 m in the Middle Pleistocene. In addition, the paleontological remains found there do not support this age (Figure 7; Sesé & Ruiz, 1992; Morales et al., 1993).

The chronological datasets available for the terraces of the Lozoya valley, and the correlation of this sequence with the main watercourses of the basin (Henares, Jarama and Tajo rivers, Figure 7) (Pérez-González, 1994; Santisteban & Schulte, 2007; Benito-Calvo & Pérez-Gonzalez, 2010; Silva et al., 2017; Karampaglidis et al., 2020) allow us to place the terraces perched up to +3-5 m in the Holocene, and between this and +13-15 m in the Upper Pleistocene. In the Middle Pleistocene, the terraces would be those located from T21 (+11-14 m) up to T16 (+50-55 m), since the level T15 (+60-65 m) has been placed in the Matuyama chron (Pérez-González et al., 2012; Silva et al., 2017), showing reversal of magnetostratigraphic polarities in the Jarama valley. These chronologies are consistent with ages obtained in other valleys in Central Spain (Arlanzón River, Duero Basin), where terraces at +70-78 m have been dated at about 1.14 Ma, levels at +60-67 m at 0.85 Ma, and terraces at +50-54 m at 0.67 Ma (Moreno et al., 2012). The sequence of negative polarities in the Jarama valley could place the terraces T9 (+120-124 m), T10 (+109-114 m) and T11 (+100-104 m) between 1.24-1.77 Ma (Figure 7), whereas the normal polarities obtained in the Lozoya fluvial terrace located in the Reguerillo Cave (Torres et al., 2005), could place terrace T7 (+142-144 m) in the Olduvai subchron, providing older Early Pleistocene ages for the upper terraces. Finally, T1 (+200-205 m), T2 (+190-195 m) and T3 (+180-184 m) could be assigned to the beginning of the Early Pleistocene where terraces close to the junction with the Jarama River are dated between 2.58 Ma and 2.3 Ma (Karampaglidis et al., 2020).

Asymmetry of drainage basins can be used in homogeneous areas to identify rivers’ lateral migrations driven by lithology-structural changes, tectonic activity and climatic changes (Cox, 1994; Garrote et al., 2008; Mink et al., 2014).

where Da is the distance from the drainage basin midline to the active channel and Dd is the distance from the basin midline to the basin divide.

T-index values vary from 0 to 1, where values close to 0 indicate perfect symmetry whereas those nearer 1 mean maximum asymmetry. The T-index can detect basin asymmetry by quantity (0 to 1) and the direction of lateral stream migrations (Cox., 1994; Garrote et al., 2008; Vijith et al., 2017). In this study, the T-index was applied only to the total length of the Lozoya River (Figure 5).

In order to assess lithological control over all these morphometric parameters, the discriminant analysis statistical technique was also applied.

This is a shape index for the valley cross-section and describes the degree of maturity of a valley (Bull & McFadden, 1977).

where Vfw: width of valley floor, Eld: elevation of the left part of the valley, Erd: elevation of the right part of the valley, and Esc: elevation of the valley floor. For Vfw, we have applied the maximum distance between the preserved fluvial terraces on each side in transverse profiles. For Esc, we have utilized the height of the fluvial base level, and for Eld and Erd, the difference in elevation between the terraces considered in each period (Figure 8 & 9).

Typically, this index has been used to distinguish V-shaped valleys (Vf values close to 0), usually associated with linear and rapid fluvial incision related to rapid uplift, from U-shaped and flat-floored valleys (high Vf values), related to less incision, and favoring infilling and sedimentation processes (Seong et al., 2008; Pedrera et al., 2009). In our study, Vf values were calculated in nine sections along the Lozoya River (Figure 8). Our goal was to apply this index along the entire longitudinal profile in all the areas. Based on reconstruction of the Lozoya fluvial base levels defined by the terrace sequence, we were able to calculate the Vf index for five different periods (P1, P2, P3, P4 & P5) of valley formation, developed between five fluvial terraces (P1: T1, T2, T3, T4 & T5; P2: T6, T7, T8, T9 & T10; P3: T11, T12, T13, T14 & T15; P4: T16, T17, T18, T19, & T20; P5: T21, T22, T23, T24 & T25; Figure 9). Next, we related this data to the river incision rates extracted from the Lozoya River terrace sequence (P1≈0.055 mm/a, P2≈0.07 mm/a, P3≈0.065 mm/a, P4≈0.006 mm/a & P5≈0.1 mm/a). In cases where the terraces were not preserved along the cross-section, we interpolated their values from the closest preserved terraces (Benito-Calvo & Pérez-González; 2008). This method allows us to analyze the spatio-temporal incision controlled mainly by bedrock characteristics, uplift processes and climatic conditions (Figure 9; Table 2).

The SL index has been widely used as an xy-profile to identify anomalies in stream profiles (Hack, 1973; Silva et al., 1988; Fernández-García & Garzón Heydt, 1994; Garzón Heydt et al., 2002; Garzon et al., 2012). In this study, we created an SL map using 5-m resolution DEM and GIS, following the method proposed by Font et al., (2010). From the DEM, we extracted a slope percentage map, which was later used to calculate the SL index along the longitudinal profiles of the drainage network. These data were used to interpolate a continuous stream length–gradient surface by applying ordinary kriging (Figure 10).

The longitudinal profiles of streams have always been considered very important data for hydrological and geomorphological studies. The profile convexities are identified as slope changes or escarpments in the streams, and are associated with disequilibrium and nonsteady states (Hack, 1973). They are usually correlated with base level lowering, lithological variations, structural trending, sediment inputs, nonfluvial processes and human modifications (Richards, 1982; Knighton, 1998; Burbank & Anderson, 2001; Schumm, 2005; Gasparini et al., 2007; Phillips & Lutz, 2008; Benito-Calvo et al., 2015). On the other hand, concave-up shapes in longitudinal stream profiles are usually associated with flat and steady-state areas (Leopold, 1994; Sinha & Parker, 1996; Morris & Williams, 1997; Smith et al., 2000; Larue, 2008; Phillips & Lutz, 2008). This assumption is based on stream power theory in that upstream the erosion and slope are increasing while downstream the gradient and discharge energy are decreasing (Smith et al., 2000; Snyder et al., 2000; Roe et al., 2002; Duvall et al., 2004; Goldrick & Bishop, 2007; Phillips & Lutz, 2008). However, we have to be careful with this notion because nonsteady conditions can also produce a smooth concave longitudinal profile (Snow & Slingerland, 1987; Ohmori, 1991; Sinha and Parker, 1996; Whipple, 2004). The analysis of the profile allowed us to identify the major irregularities in the stream longitudinal profile (Figure 11).

The analysis of longitudinal profile irregularities was extended to longitudinal profiles of the Lozoya terraces, to check for the persistence of the anomalies through time (Figures 10 & 11). To do so, the terrace longitudinal profiles were reconstructed using the n ^th degree polynomial functions (Benito-Calvo et al., 2008), which displayed the best fit (R ² > 0.99). Additionally, to examine whether these knickpoints are pre–Quaternary, we reconstructed the bedrock paleotopography of an extensive and well-preserved Turolian coarse-grained alluvial fan system (Profile G3; Figures 1 & 10; Bellido et al., 1991b; Portero et al., 1990; Karampaglidis et al., 2020). These deposits lie in the Atazar and Ponton de la Oliva sectors, parallel to the Lozoya River and always in the NE part of the modern flow (Figure 1). To reconstruct the paleoprofile we used the contact of the deposits with the Paleozoic basement and Early Miocene sediments as a reference point (Figures 1 & 10).

The Concavity Index was used to analyze the general shape of the longitudinal profiles (Table 3). It was computed on the basis of deviations from a straight line (Phillips & Lutz, 2008):

where Hi is the elevation at distance i, Hi* is the elevation along a straight line from the uppermost to lowermost point along the stream line at horizontal distance i, and N is the total number of measurement points. Negative CI values indicate convexity, whereas positive CI values indicate concavity. Convexities in long profiles are typically associated with base level lowering, litho-structural controls, sediment inputs, nonfluvial erosion processes, bedload fluctuations and human modifications, while smooth concave-up long profiles have long been considered as steady-state equilibrium forms (Phillips & Lutz, 2008).

The origin of asymmetric basins could be associated with lithological-structural variation and/or tectonic activity. The T-index is effective for highlighting tilted areas (Garrote et al., 2008). In our case, the transverse topographic asymmetry values display a clear spatial distribution (Figure 5). Small to medium magnitudes are identified in the Pinilla del Valle and Buitrago de Lozoya areas with three prominent directions of river migration: SW, NW and SE. On the other hand, in the Atazar area the T-index values are higher, between 0.3 and 0.5, with only one prominent migration direction (E). Finally, the highest values (>0.5) are identified in the Ponton de la Oliva area, with north and east lateral migrations of the river direction (Figure 5).

In the case of the high asymmetry values in the Atazar area, no geological evidence of neotectonic activity on the preserved Turolian sediments or on the Quaternary landforms has been identified (Portero et al., 1990; Bellido et al., 1991a; Andeweg et al., 1999; Karampaglidis, 2015). Probably, the high values are related to lithological changes delimited by the Lozoya River, located between gentle western slopes on Permian granites and rugged terrain developed on Paleozoic metasediments on the eastern side (Figure 1 & 4). This is a very common feature in the MCB (Garrote et al., 2008), where Miocene detrital sedimentary rocks with extended terraces are located on the western margin, whereas gypsum escarpments characterize the eastern margin (Silva et al., 1988). Pontón de la Oliva is the other sector with high values. In this case the lithology is homogeneous (Figure 1) and high T values cannot be explained by bedrock lithological changes. Just as in the previous sector, the left margins of the Lozoya River are characterized by Turolian detrital sediments with no tectonic evidence of large block tilting (Portero et al., 1990; Bellido et al., 1991a; Andeweg et al., 1999; Karampaglidis, 2015). In the MCB, Garrote et al., 2008 have related the basin asymmetry of the Jarama, Henares, Tajuña and Tajo rivers to the principal morpho-structural patterns of the basin, influenced by the crustal undulations previously defined for the basin. At this point, it is worth mentioning that this area stands out for the Atazar and Tortuero backthrust structures with N-S trend and left-lateral strike-slip movement up to the Somosierra Pass Fault (Figures 1 & 4; Vicente et al., 2007). The Lozoya River general trend in this area is ENE, downcutting parallel to the Atazar backthrust structure, which is a relatively depressed area and almost perpendicular to the NNE direction of the center of the Madrid basin (Figures 1 & 4). Similarly, to Garrote et al., 2008, we can attribute the high T values to previous morpho-structural trends in response to Alpine tectonics and to the old Lozoya River capture (Karampaglidis et al., 2011). The latter would have produced a reorganization of the local drainage in this sector. The stream capture of an asymmetric basin could lead to an acceleration of the river’s incision without external forcings (Prince et al., 2011) and when entrenchment continued, the Lozoya River would have incised and captured the local streams running SW, the same direction as the center of the MCB (Figures 2(b) & 4) (Karampaglidis et al., 2011).

The incision rates and Vf index are very useful tools for detecting areas with base level lowering. In order to analyze the temporal and spatial variations, both vertically and horizontally, of the Lozoya -valley, first of all, we calculated incision rates for the whole basin during the five different periods (P1, P2, P3, P4 & P5) and secondly, we applied the Vf index in the different sectors and to different stages (P1, P2, P3, P4 & P5) of the valley development (Figure 8; Table 2). The incision rates vary from 0.06 to 0.1 mm/a and two periods of high acceleration are the first thing we notice (Figure 9; Table 4). For the Vf index the lowest values (close to 0) are located from the D-D’ to I-I’ cross-sections. In these areas the valley shape is asymmetrical and steep (Figures 2(a) & 7). Overall, the Lozoya valley shows a general tendency to develop a more pronounced V-shaped geometry with time. This is an obvious behavior in a downcutting staircase valley, where uplift processes are accepted as triggering incision and as generating the differences in height between fluvial levels (Westaway et al., 2009).

This estimation of the incision rates dynamic during the Quaternary permitted us to compare the long-term incision rates of the Lozoya River with other studies of staircase fluvial sequences such as: Duero (Central-North Spain; Silva et al., 2017; Rodríguez-Rodríguez et al., 2020), Upper-Lower Tajo (Central-South Spain; Cunha et al., 2008; Pérez-González 1994; Silva et al., 2017; Karampaglidis et al., 2020) and Allier (Central-North France; Pastre 2005). For P1 up to P2, the Lozoya River shows high incision rates (0.15 mm/a) and a large Vf index, between 180 and 80. The highest values of Vf are located in the Buitrago de Lozoya and Atazar areas. Similar incision rates for this period have identified for the Jarama River extended alluvial plains (~0.15 mm/a) and related with climatic changes and the capture of the MCB from the Atlantic drainage system (Karampaglidis et al., 2020). This anomaly is not detected for the Lower Tajo basin (Figure 9). For this period, extended erosive plains have been identified in the Buitrago de Lozoya and Atazar areas (Karampaglidis et al., 2014a). This probably explain the high values of Vf at the start of the P1 period. For the periods P2 and P3, the higher Vf values are identified in the Buitrago de Lozoya area, between 80 and 20, while the lower ones are in the Ponton de la Oliva area, between 20 and 5 (Figure 9). In this period, the Buitrago de Lozoya and Atazar areas are still asymmetric, while the Ponton de la Oliva sector is characterized by steep valleys (Figure 8). Asymmetries and steep valleys, applying the Vf index, are also observed in the Duero basin on the metamorphic rocks and granites and attributed to lithological variability and drainage capture (Antón et al., 2014). During these periods, valleys are still steepening, while the incision rates are falling from 0.09 mm/a to 0.06 mm/a (Figure 9). This indicates no abrupt changes in incision rate. For these periods, this trend has been identified in other large tectonically quiescent Atlantic drainage basins such as the Duero (Central-North Spain), Upper-Lower Tajo (Central-South Spain) and Allier (Central-North France) with values ranging between 0.05 and 0.075 mm/a (Figure 9; Table 4). Finally, for P4 and P5 an acceleration of the incision rates up to 0.175 mm/a has been detected and coincides with the falling trend in Vf (Figure 9). During these periods, the Buitrago de Lozoya, Atazar and Ponton de la Oliva areas are attaining maximum valley steepness (Figure 8). We can easily observe the same trend in the other basins, with the Duero and Allier basins at 0.1 mm/a and 0.115 mm/a respectively. Silva et al., 2017 related this phenomenon with Middle Pleistocene transition driven by climatic fluctuations. For the Lower Tajo basin, the same pattern is repeated but with a much more dynamic rise in values (up to 0.6 mm/a). Cunha et al., 2008 suggested that these high values are driven by an episode of crustal uplift while Silva et al., 2017 proposes eustatically forced controls.

However, not all the sectors show exactly the same dynamic: the Vf curve in the Upper Lozoya Valley (Pinilla del Valle sector, Figure 9) shows a different pattern for the periods P2 and P4. Vf values are stable and stand between 50 and 40. Most likely, the pop-down NE-SW morphology of the SCS crystalline basement, delimited by lithological contrast between Paleozoic gneisses and Cretaceous dolomites to the south and Paleozoic gneisses with Tertiary detrital sediments to the north (Figure 1) maintained stable shape ratios from the Early Pleistocene (P2) to the Middle Pleistocene (P4). The lithological contrast between high rock strength gneisses with highly soluble and fractured limestones and low strength sediments probably favors valley opening, lateral planation and smoothing from weathering and erosion. Schanz & Montgomery, 2016 observed that slaking weathering and crystalline transported material as bedload over weaker sedimentary units acts as an abrasion tool and foments rapid valley widening. Seen in this light, the widening of the valley could generate large values for Vfw, and this keeps Vf stable and yields a U-shape effect. During the Middle Pleistocene an extended karstic system developed on the local Cretaceous dolomites which could foster this (Pérez-González et al., 2010). Subsequently, in the Upper Pleistocene (P5), the valley shape changed sharply, coinciding with the development of the youngest strath terrace and small alluvial fans. In this case, Vf presents the same dynamic as in the other sectors.

The interpretation of the Vf and T-index values was combined with the information provided by the stream longitudinal profile analysis. The SL index represents the slope variation along the stream profiles (Hack, 1973; Font et al., 2010). Whereas high values are normally located in upstream areas and could perfectly well be related to drastic changes in relief due to erosion, further downstream, they are very commonly associated with knickpoint areas (Hack, 1973; Garrote et al., 2002; Pérez-Peña et al., 2009; Font et al., 2010; Garzón Heydt et al., 2012; Antón et al., 2014; Soria-Jauregui et al., 2018). In the MCB, Silva et al., 1988 related high SL values for the Jarama and Manzanares rivers with bedrock variability and neotectonic activity associated to halokinetic and subkarst dissolution processes of the bedrock, together with the reactivation of large basement faults. Other authors such as Pérez-Peña et al., 2009 demonstrate that in the Granada Basin, high SL values could be related to the large drainage areas, while in the Duero Basin, Antón et al. (2014) attribute them to drainage reorganization prompted by the capture of the endorheic basin from the exorheic Atlantic drainage system. In the present study, the stream length–gradient index (SL) was calculated including the main stream profiles of the Lozoya watershed, interpolating to create a map for the entire study area (Figure 10). This map shows certain zones in the lower part of the Lozoya valley with high SL values (Figure 10). The first anomalous zone (B1; Figure 10), defined by intermediate SL values, was where the Lozoya River abandoned the upper valley pop-down depression, cutting a NE-SW structural alignment (profile C-C´; Figure 8). In this case it seems to be influenced by the Alpine structure trends and bedrock lithological contrast between Cretaceous limestones and Paleozoic gneisses (Figure 1& 10). Towards the north, Lozoya tributaries have high SL values (zones B2 and B3) coinciding with canyons incised following a N-S band of schist and slates in structural contact with gneisses (Figure 1). Bedrock lithology contrast and structural contact between metasediments and metaplutonic rocks can generate weakened bedrock areas with differing response to erosion by local streams (Antón et al., 2014). In the Atazar sector, the SL index allowed us to distinguish a large zone with steep gradients along the Lozoya canyon carved into the lithological-structural contact between the orthogneisses and leucogneisses, and fractured metasediments and granites (zone A1, Figures 1 & 9). Below the Atazar Dam (Cerro de la Oliva sector), the Lozoya River again shows high SL values (zone C1). This coincides with the elbow capture where the old Lozoya tributaries and the Lozoya River were captured (Karampaglidis et al., 2011), causing an increase in gradient here and adjustment of the tributaries to the new base level (zone C2; Figure 10). As we can observe, the SL index is a powerful tool for identifying topographical breaks along a stream (knickpoints), but further considerations are needed to interpret it.

Recognition of the most representative anomalies was facilitated by analysis of the slope and curvature variability along the Lozoya longitudinal profile (Figure 11). The current Lozoya longitudinal profile shows the normal tendency to decrease in concavity downstream (Table 3), although with several anomalies. The headwater of the Lozoya River has several minor irregularities, shown by a high variability in the slope and curvature profiles. This variability profile decreases sharply in the pop-down Pinilla del Valle depression, and rises again in the Middle and Lower Lozoya Valley (Buitrago de Lozoya, Atazar and Cerro de la Oliva sectors). Within these sectors, slope and curvature profiles show several major knickpoints, which coincide with the high SL zones detected for the Lozoya River and its entrenchments (zones B1, B2, A1 & C1; Figure 10 & 11). Some of these knickpoints appear related to faults (B1, B2 & A1), and differential uplift associated with Quaternary fault reactivation could be suggested. To assess this point, we analyzed the persistence of these knickpoints over time, by reconstructing the longitudinal profiles of the Lozoya terraces and the paleotopography where the Turolian detrital sediments have been deposited (Figures 10 & 11). Data from terraces have coarser resolution than the modern Lozoya longitudinal profile, but the trend of the terrace profiles and the terrace height data were useful in providing insights into the persistence of the major knickpoints during the downcutting process. The convex morphology of the B1 knickpoint can be recognized in the terrace profiles developed from the Lower-Middle Pleistocene to the present day (T13-T24), whereas earlier terrace profiles do not provide clear evidence to support the existence of this knickpoint. Nevertheless, the straight or concave morphologies shown by the oldest terrace profiles in the B1 zone might be artifacts of the limited data available to interpolate the oldest terraces. This situation is similar for the knickpoint A1, easily identifiable in terrace profiles T13-T24, but with a very weak signal in older terraces. For example, in terrace profiles T2-T5, the general trend of the interpolated curve does not reflect knickpoint A1, yet it can be recognized in the local data shown by the terrace heights (Figure 12). Further, the reconstructed GS3 profile suggests that the presence of the A1 main trend is clearly pre-Quaternary, since at least the Turolian, and the Lozoya River is smoothing it out (Figure 11). In addition, it should be noted that there is no evidence of postfaulting or dipping for these deposits (Portero et al., 1990; Bellido et al., 1991a; Andeweg et al. 1999; Karampaglidis, 2015). This could suggest that the A1 knickpoint origin is related to a response to Alpine tectonics and/or neotectonic reactivation from Late Miocene active thrusting, whereas the southern border of the SCS was covered by Turolian sediments (Andeweg et al. 1999; Vicente et al., 2007; Karampaglidis et al., 2020).

Analysis of the stream longitudinal profiles indicates that the main knickpoints detected have an Alpine origin associated with the trends of the SCS principal structures in this area and controlled by lithological contrast and stream capture. Additionally, reconstructions of terrace profiles show no significant differential movements associated with fault reactivation from at least the Early-Middle Pleistocene boundary to the present day.

Combining the Lozoya terrace records and geomorphometric indices we observe that the Lozoya drainage basin evolution seems to be conditioned by Alpine morphostructural trends, lithological contrast, river captures assessed both at the local and basin level, climatic fluctuations and denudational uplift. Thus, at this point the mix of these patterns could explain the Lozoya basin’s erosive staircase style and the poor preservation of the terrace alluvial cover in some sectors. The shape, size and preservation of bedrock terraces depend on: climatic variations which configure sediment supply and river discharge (Hancock & Anderson, 2002; Turowski et al., 2007, 2008; Yanites & Tucker, 2010); base level lowering caused by sea level changes or tectonics (Merritts et al., 1994) or river capture (Stokes et al., 2017); lithological contrasts (Schanz & Montgomery, 2016; Stokes et al., 2017) and internal forcings such as meander cutoffs (Finnegan & Dietrich, 2011).

In the Lozoya drainage basin, a long staircase sequence was identified (24 levels), with 15 Early Pleistocene, 5 Middle Pleistocene and 4 Late Pleistocene levels (Karampaglidis et al., 2014a). The vertical separation between the terraces is ~10-15 m for the Early Pleistocene ones, while for the Middle-Late ones it is lower: ~3-7 m. The sequence is mainly composed of strath and bedrock terraces. Bedrock, or strath, river terraces are formed when there is a change in ratio from vertical incision to lateral erosion (Hancock & Anderson, 2002). When rivers with an alluvial cover migrate across valley bottoms, the underlying bedrock is eroded into a planar surface called a strath (Personius et al., 1993; Wegmann & Pazzaglia, 2002; Fuller et al., 2009; Finnegan & Balco, 2013; Pazzaglia, 2013; Schanz & Montgomery; 2016). In the Buitrago de Lozoya and Atazar areas, the river cut through Paleozoic gneisses and granites, and the identified terraces are mainly erosive, with poor preservation of the alluvial cover. Meanwhile, in the Upper Lozoya and Ponton de la Oliva sectors strath terraces 2-4 m thick are recognized, with massive structure and oligomictic material. In these sectors, the Lozoya River cuts through Cretaceous dolomites or Miocene detritals sediments and black slates respectively. So, it is obvious that the bedrock geology influenced the preservation of the alluvial cover. Indeed, this is very common: Montgomery (2004) observed that alluvial cover formation and/or preservation is easier on less resistant lithologies like sandstone and siltstone than more resistant ones such as quartzite and basalt. Along the same lines, Stokes et al. (2017) contend that underlying bedrock geology and its stratigraphic and structural configuration play a role just as important as other mechanisms (lowering and climate) in terrace formation and/or preservation. Furthermore, Schanz & Montgomery (2016) showed that bedrock lithology is a very important factor in the rates of lateral erosion influencing floodplain width and the potential for bedrock planation and bedrock cover preservation. In the same study, they notice the importance of physical weathering from wetting and drying in alluvial cover preservation. In our case, in the Atazar and Buitrago de Lozoya sectors, the Paleozoic granites and gneisses bedrock has a high mica content and is prone to physical weathering (Bellido et al., 1991b; Karampaglidis et al., 2014a). In addition, as we have seen, in these areas we identified three pre-Quaternary knickpoints (A1, B1 & B2; Figure 10) which break the Lozoya River’s longitudinal profile abruptly and thus produce acceleration and higher stream power, so that stronger vertical incision could take place locally. In addition, Schanz & Montgomery (2016) observed that strong vertical incision creates erosive point banks perched above the water table and exposed to mechanical and chemical weathering processes. That means that on rocks prone to weathering during high flows with crystalline bedload, the weathered material is rapidly washed away, perhaps resulting in centimeters of lateral erosion in a single event (Montgomery, 2004; Stock et al., 2005; Collins et al., 2016). Thus, fast-weathering lithology with existing morphostructural longitudinal profile anomalies could explain the poor preservation of the alluvial bedrock cover in these areas. At this point it should be mentioned that more research along these lines is needed, and a detailed study of erosion rates and surface dating will probably help us improve our understanding.