Identification and geomorphic characterization of fluvial knickzones in bedrock rivers from Courel Mountains Geopark

Vol.:(0123456789) 1 3 Environmental Earth Sciences (2023) 82:475 https://doi.org/10.1007/s12665-023-11098-5 ORIGINAL ARTICLE Identification andgeomorphic characterization offluvial knickzones inbedrock rivers fromCourel Mountains Geopark HoracioGarcía1,2 Received: 27 March 2023 / Accepted: 9 August 2023 / Published online: 20 September 2023 © The Author(s) 2023 Abstract The gradient characteristics of Courel Mountains Geopark bedrock rivers were examined. Unlike alluvial rivers, bedrock rivers have been the great forgotten of fluvial geomorphology globally. Based on the decreasing rate of gradient with increasing measurement length, a relative steepness was obtained as indicator of knickzone. Supported by GIS techniques and DEMs, the changes in slope along the longitudinal profile of the rivers were detected. The number of the extracted knickzones rises to 325, which means a frequency of knickzones of 0.467 km−1. The total length of the knickzones is 285km, representing about half of the drainage network as knickzone (47%). The mean height, the length, and the gradient of all the knickzones were ~ 110m, ~ 880m, and 0.178m·m−1, respectively. There is no distribution pattern and the knickzones are everywhere, although they are more present in reaches with NW–SE direction and order 1. Several environmental factors were crossed to know more about the occurrence and knickzones characteristics, suggesting that density and direction of fractures regulate the number and the trajectory of the knickzones, while the lithology controls the singularity of the forms. The geomorphological and the topographical characteristics of the bedrock rivers make them high ecological, scenic, landscape, and recreational value. Findings from this study can be also used by managers to develop and/or improve strategies for conservation, valorisation, and how to approach the tourist who visits the Geopark. Scientific tourism can offer a unique and educational travel experience, allowing participants to learn about bedrock rivers and knickzones. Keywords Fluvial geomorphology· Knickpoint· GIS· Geotourism· Galicia Introduction Forms and process in bedrock rivers reflect the interactions between erosive dynamic and the resistance of the substrate through flow. Under a theoretical framework, bedrock rivers are characterized by erosion and transport taxes more important than production (i.e., capacity to generate sediments) and storage (Wohl 2015). Conversely, in alluvial rivers, production and storage are more present than erosion and transport. On these theoretical conceptions, there are a multitude of intermediate situations derived from the energy and granulometric variety of bedrock rivers (Ortega 2010). The geomorphological and the topographical singularities of the bedrock rivers make them high ecological, scenic, landscape, and recreational values (García etal. 2021; Ollero 2017). Unlike alluvial rivers, bedrock rivers have been the great forgotten of fluvial geomorphology globally (GarzónHeydt 2010). Longitudinal profiles have been one of the most used tools for the characterization of bedrock rivers due to their ability to describe their geomorphological evolution through, basically, detecting changes in slope. Following the proposal launched by Hayakawa and Oguchi (2009), for this work, I define these geomorphological attention points as knickzones: “locally steep riverbed segment”. This definition also encompasses the term knickpoint, generally defined as referred to the exact location where there is a sharp change in channel slope (downstream increases) (Bierman and Montgomery 2013; Tinkler 2004). In geometric terms, a knickpoint is a single record represented by a point while a knickzone refers to a linear entity. In the first case, the topographic change location is recorded; in the second a change transition. * Horacio García [email protected] 1 Department ofGeography, University ofSantiago de Compostela, Galicia, Spain 2 Ambiosol Research Group, University ofSantiago de Compostela, Galicia, Spain Environmental Earth Sciences (2023) 82:475 1 3 475 Page 2 of 19 Bedrock rivers from Courel Mountains Geopark (onwards, CMG) (Galicia, NW of Iberian Peninsula) are one of the geo-assets that favored its designation as UNESCO Global Geoprak in 2019. However, for this region from Spain, there are still few works on both fluvial geomorphological dynamics and bedrock rivers although some of the descriptive types stand out (Pérez-Alberti 1982, 1985, 2000), others on granitic forms (Álvarez-Vázquez and De Uña-Álvarez 2021; De Uña-Álvarez etal. 2014, 2009), or others that include bedrock rivers as part of classification processes (García, 2014). The reference work on bedrock rivers of the Iberian Peninsula (Ortega and Durán, 2010) does not include Galicia, possibly due to this lack of studies on bedrock rivers (despite having a rich and varied representation). This gap, not remedied in the last decade, must be a challenge and a need to which this article tries to contribute. The main objective of this paper is to identify knickzones from CMG bedrock rivers following the GIS procedure proposed by Hayakawa and Oguchi (2006) for Japanese mountain rivers. To achieve this aim, it is necessary to set the rivers of the study area a topographical criterion to define the knickzones in terms of stream gradient threshold and reach length. Subsequently, the distribution control factors and the geomorphological characteristics of the knickzones are analyzed. By geomorphological characterization, I mean both topographical, tectonic, and lithological aspects. Study area CMG is in the Northwest of the Iberian Peninsula (Fig.1), inside of the limits of Galicia region (Spain). The total area of the study zone is ~ 600 km2. 73% belongs to the CMG, which one is not included at whole (southern zone—more or less south of Quiroga village—was removed because the predominance of Cenozoic rocks and lower slopes). In order not to divide the basins that partially form part of the CMG, they were fully incorporated as part of the study area. These zones represent the remaining 27% and correspond to the “surrounding area” (north and west flank of the Lor-Lóuzara river basin). From here, I will only use the acronym CMG to refer to the study area, which one is made up of the river watersheds: Lor, Parteme, Quiroga, Castelo das Eiras, Soldón, Casas, and Selmo (only river head) (Fig.1). Fig. 1 Location of the study area and main physical characteristics. Geological map shows slate and quartzite as the lithology dominant, but there is also the presence of another metamorphic (e.g., schist) or sedimentary rocks (e.g., sandstone). Fracture is representing only the main Environmental Earth Sciences (2023) 82:475 1 3 Page 3 of 19 475 Horacio etal. (2018) developed a classification for Galicia based on areas with similar topography and lithology and, in consequence, analogous geomorphological processes, named lithotopo units (see also Montgomery 1996). According to these authors, CMG belongs to the lithotopo unit named A3–A2, characterized by an area with compact materials (metamorphic rocks, mainly slate and quartzite bands), marked orientation and a relief with high slope and roughness. The average elevation of CMG is close to 900m (~ 225min, ~ 1650 max, ~ 1400 range), with an average slope of ~ 25degrees (°). Table1 contains detailed elevation and slope information by basin and on the fluvial network, which is dominated by a dendritic pattern, but with some subdendritic connotation or even areas with rectangular patterns. The average length of the main river by basin is 19.2km, oscillating on a length range from ~ 4 to 56km. The hydrological regime is a mix between the oceanic pluvial regime and the Mediterranean regime, except for the most mountainous areas that have some snowfall. The average precipitation ranges between 1100 and 1400l m−2 per year, and an annual average temperature of ~ 10–12ºC. The curve prototype of the hydrological regime is characterized by an intense period of high waters in winter (mainly from November to February), followed by a slight second regrowth at the end of April–May, which serves to slow down the languid drop little bit in flow experienced during spring. The flow reaches its most critical point in summer (from July to September). The average flow of the main river of the CMG (Lor) (Fig.1) is ~ 14 m3 s−1, with a monthly maximum and minimum in the series of 64 m3 s−1 and 1.7 m3 s−1, respectively. CMG is a rural area in clear demographic decline, with scarce urban centers (such as small villages) distributed throughout its territory. This decline caused many cultivated meadows to be abandoned and a process of plant succession began. Currently, CMG has a wide expanse of native forests (oak, chestnut, holm oak) dotted with shrubby undergrowth and meadows. Methodological framework Previous considerations Methodological workflow is structured in two stages. The first one focuses on defining the relative steepness index. Here, I employed several rivers from the lithotopo unit to which CMG belongs. The second stage looks toward calculating the stream gradient on the CMG rivers to identify the knickzones and analyzing their control factors and geomorphic characteristics. A geodatabase for the whole NW Iberian Peninsula was generated, so the study area and other zones used in this paper were clipped from the geodatabase. Data Digital elevation model A 5 × 5m digital elevation model (DEM) was done by creating a mosaic with 81 DEMs-sheets at 1:50,000 scale (50K) from the Spanish National Geographic Institute (IGN). Then, DEM-5 × 5m was resized by means of a bilinear procedure to new raster with pixel size 25 × 25m and 50 × 50m. The different spatial resolutions were used to determine the most appropriate DEM according to draw the drainage network and detect knickzones. Lithology The lithological data used correspond to the named MAGNA series of the Geological and Mining Institute of Spain (IGME) at 50 K. The original cartography contains 13 lithological units for the study area classified according to different variables (e.g., rock type, age, degree of permeability, and alterability). The data was supplied in a GIS vector layer. Table 1 Main hydrographic and topographical characteristics by basin of the study area Basin Area (km2) Elevation (m) Slope (o) Length main river (km) Length fluvial network (km) Min Max Range Mean Mean Lor 366.8 225.4 1639.5 1414.1 891.7 17.9 56.2 364.2 Soldón 90.9 254.2 1540.4 1286.2 879.7 20.1 23.7 85.3 Quiroga 80.4 236.1 1607.1 1371.0 795.0 24.3 22.6 81.3 Selmo (river head) 37.5 784.7 1637.7 853.0 1205.6 27.3 9.0 35.1 Parteme 16.5 233.3 1123.5 890.2 700.9 24.9 9.0 16.5 Castelo das Eiras 10.0 235.7 1174.1 938.4 700.5 27.0 10.3 12.2 Casas 9.4 284.2 1132.4 848.3 807.9 23.2 3.9 11.4 Environmental Earth Sciences (2023) 82:475 1 3 475 Page 4 of 19 Fluvial network Fluvial network was developed through GIS hydrological tools (ArcGIS PRO) following the common steps to get the drainage network and basins from a DEM. Initially, I used a DEM with cell size 25 × 25m. The reason for this resolution choice is because it fits well to the drawing of the river from national topographic map 25K. Besides, DEM-5 × 5m generates a fluvial network with excessive detail, while DEM50 × 50m is too rough. Briefly, those steps were: (1) filling of holes or small depressions to eliminate possible imperfections in the data. (2) Create a flow direction raster from each cell to its neighbor(s) with the steepest downward slope (method Deterministic 8, D8). (3) Creation of an accumulated flow raster in each cell following the previously calculated flow directions. (4) Using map algebra, a threshold value of 20,000 accumulated pixels was set as the start of the drainage network delineation, that is, the channel drawing will be made from those cells that have more than 20,000 cells flowing toward them. The designation of the threshold value does not follow any analytical method (e.g., Tarboton etal. 1991), but a visual comparison in which the density and traced of the drainage network calculated was like the drainage network of the 25K national topographic map. (5) Determine the contributing area over a set of cells from a raster (i.e., catchment). (6) Conversion of the accumulated flow to a linear vector layer and the basins to polygonal vector layer. The reaches artificialized by reservoirs were eliminated (source: IGN). Also stream reaches over Cenozoic sediments (source: IGME). This drainage network was used to identify the relative steepness index. Selection ofrivers todefine therelative steepness index Rivers selected to define the relative steepness index (first stage of the methodological process; see Sub-section.“Previous considerations”) had to meet two conditions: (i) belonging to the same lithotopo unit where is the CMG rivers, and (ii) have a “high probability” to be a bedrock river. First point was explained in the study area section. Regarding the second one, I followed a previous work in which the probability of bedrock rivers in the drainage network of the NW Iberian Peninsula was identified (García and PérezAlberti 2021). These authors applied a method supported by GIS to discriminate (at 102–103m scale analysis) between potentially bedrock rivers, and those that are potentially not, through the topo-geomorphological descriptors: type of lithology, surroundingslope, and fractures density. The results offer three levels of bedrock reach probability (high, medium, and low) for the whole Galician fluvial network. Finally, a total of 25 rivers with both conditions were selected to define the relative steepness index (Fig.2). The results achieved were used later to apply on the CMG rivers and identify the knickzones and analyze their geomorphic characteristics and control factors. The rivers selected sum 585km and combine whole rivers (i.e., from the head to the mouth) and continuous parts of rivers (i.e., from the head to reaches excluded according to criteria explained above), moving on a length range from ~ 8 to 72km and mean and standard deviation values of 23.3km and 15.2km, respectively. Identification ofknickzones I followed the method of Hayakawa and Oguchi (2006) to detected knickzones thought the KET (Knickzone Extraction Tool) GIS toolset for automatic extraction of knickzones from a DEM based on multi-scale stream gradients Veiga de Logares Rodil Eo Suarna Neira Queizón Rao Larxentes Ser Castelo Cervantes Cales Navia Narón Sarria Lor Lóuzara Cabe Saa Pequeno Quiroga Soldón Loureiro Ribeira Grande 05 km 10 Courel Mountains Geopark River studied Watershed outlet Couso Fig. 2 Rivers/reaches selected to define the relative steepness index. See Fig.1 for more information about the location. Note in this figure Courel Mountains Geopark does not include the surrounding area (see about this in Section: “Study area”) Environmental Earth Sciences (2023) 82:475 1 3 Page 5 of 19 475 (Zahra etal. 2017) (second stage of the methodological process; see Sub-section“Previous considerations”). The tool was run with python language from PyCharm. The procedure is quickly explained here because its details are described in the cited article. The calculation procedure can be summed up in three steps: (1) calculation of stream gradients with several measurement lengths to set a relative steepness index; (2) determine knickzones using a threshold value previously defined; and (3) extraction of the knickzones the general properties. Stream gradient (Gd; m m−1) at each reach is defined as (Eq.1): where e1 and e2 are the elevations (m) upstream and downstream from the midpoint (d/2) of a stream reach, respectively; d refers at the stream reach length (m). Using a regression line is possible to analyze the Gd–d relation (Eq.2): where a and b are coefficients of the linear regression. Because Gd decreases with increasing d, the negative slope of the regression line (Rd = –a; m−1) marks how the local gradient of the measurement point is steeper than its trend gradient (Hayakawa and Oguchi 2009). The Rd value applied along the riverbed draws the relative steepness. The Rd standard deviation was used to set an extraction threshold and locate those river segments as knickzones. For this study, I calculated the stream gradients along the streamline previously extracted of the DEM. The minimum distance between measurement point was 150m to ensure an interval longer than the diagonal length of the DEM resolution (for example, for the lowest resolution DEM, 50 × 50m = 70.7m diagonal length of the cell size) and reduce the presence of negative values between two consecutive measurement points. Gd values were estimated for divers d from 150 to 10,000m. A first set of values was divided at intervals of 150m (from 150 to 1500), which yields 10 types of segments with different length. The second one, it was divided at intervals of 300m from 1500 to 3000, i.e., 5 types of segments with different length. The last set of values moves from 3000 to 10,000 at intervals of 1,000m, i.e., 7 types of segments with different length. (1) G d= ( e2−e1 ) d, (2) Gd = a ⋅ d + b, Results Calculation andanalysis ofstream gradients toidentify knickzones DEM cell size adjusts The calculation of Gd was performed by extracting values from the DEM with cell size of 5 × 5m, 25 × 25m and 50 × 50m for measurement points with d = 150, d = 300, d = 450, d = 600 (Fig.3). The presence of negative values was checked as quality measurement to verify the most appropriate DEM resolution. Gd < 0 is because of adjacent valley-side slopes and/or the fact that all cells were forced to flow outward to avoid sinks in the DEM. This explains why some of the downstream measurement points could get a higher elevation value than the upstream one, which is unrealistic in terms of flow movement. According to the data in Table2, the sum of negative values of Gd increases considerably with DEM cell size but decreases with increasing d. Considering the total of Gd values obtained for the 25 rivers analyzed to define the relative steepness index and for DEM-5 × 5m, ~ 12% of the Gd values are negative for d = 150, ~ 5% for d = 300, and < 0.6% since d = 900. Based on these observations, d can be equal or larger than 150m, assuming some values (the negatives) must be removed (or averaged) for calculations of the relative steepness index. Using d = 300 the number of negative values decreases ostensibly. Figure3 shows an example of Gd calculated for the Lóuzara River (main Lor tributary; see Fig.1) with different DEM cell sizes and d = 300. Gd obtained from 50 × 50m DEM reveals a much broader threshold value than for the other two resolutions, with numerous negative cases (Fig.3A). The same pattern is repeated considering only the absolute values, but it is not a valid procedure because numerous unrealistic values are generated as potential knickzones (Fig.3B). In the case of the Lóuzara River (~ 32km), only five errors were detected due to Gd < 0, three with an elevation difference between measurement points low or very low, and two slightly higher (Fig.3C). Analysis ofstream gradient The relation between Gd and d displays a tendency from short d to mark the local slope to longer d to draw the trend gradient (Table3). Gd values calculated along each river show a large fluctuation with d ≤ 900 and smoothing with d > 900 (Fig.4). The first type reflects the local riverbed form and the second one the trend gradient. In other words, the local gradient is defined by a much stronger Environmental Earth Sciences (2023) 82:475 1 3 475 Page 6 of 19 0.14 0.02 0.04 0.06 0.08 0.10 0.00 -0.02 0.12 -1 G (m·m ) d 81624320 C Distance from riverhead (km) ErrorError -0.05 0.00 0.05 0.10 0.15 0.20 -1 G (m·m ) d -0.10 -0.20 A DEM 05x05m DEM 25x25m DEM 50x50m -0.15 0.10 0.15 0.20 0.05 0.00 -1 G (m·m ) d B Fig. 3 Stream gradients from measurement points (m) obtained each 300 m from three different DEM resolutions along Lóuzara River, where: streams gradients calculated directly from cell values extracted from the DEM, therefore negative values could be present (A); streams gradients calculated directly from cell values extracted of the DEM but as absolute values (i.e., without negative values) (B); errors detected after crossing streams gradient calculated with and without negative values (C) Environmental Earth Sciences (2023) 82:475 1 3 Page 7 of 19 475 rate of change than the gradient trend. The slope of linear regression from d = 300 to d = 3000 considering the Gd mean and standard deviation is −1.36·10–6 and 2.61·10–6, respectively. For local gradient, the slope is −3.68·10–6 (mean) and 4.87·10–6 (standard deviation), while for trend gradient, the slope is −1.03·10–6 (mean) and −2.24·10–6 (standard deviation). Internally, transition zones between gradients can also be differentiated. In the local gradient, with d ≤ 600, the rate of change is very steep, and it decreases until d = 900. In the trend gradient, since d = 1500, the rate of change is softer. Identification ofknickzones throughtherelative steepness index Following Hayakawa and Oguchi (2006) method, the standard deviation of Rd for each river was calculated. Its average for all rivers was 5.05 10–6, so river segments with values of Rd larger than this value were classified as knickzones (Fig.4). Height less than 10m should be removed to ensure an interval longer than the diagonal length of the DEM resolution (5 × 5m = 7.1m diagonal length of the cell size). Note that the second package of the KET GIS toolset (“Extraction of knickzones”), was applied using 300m like knickzone length threshold (d) and DEM-5 × 5m to minimize the number of negative values (see Table2). Knickzones were drawn under a flow accumulation threshold of 20,000, following the same criteria as the fluvial network done to define the relative steepness index. The use of a DEM with more detail (5 × 5m) than the one used to calculate the relative steepness index is also due to (i) CMG area and rivers are smaller; and (ii) presence of waterfalls, gorges and canyons everywhere, which requires a higher degree of finesse to detect them. Description ofknickzones fromCourel Mountains Geopark The number of the extracted knickzones rises to 325. The total length of the knickzones is 285km, which represents about half of the drainage network as knickzone (47%), and a frequency of knickzones of 0.467 km−1. The mean height, the length, and the gradient of all the knickzones are ~ 110m (min.: 10m; max.: 507m), ~ 880m (min.: 305m; max.: 3,463m), and 0.178m·m−1 (min.: Table 2 Sum of negative values of Gd for different DEM cell size and value of d (Navia River data) Cell size d = 150 d = 300 d = 450 d = 600 DEM05 97 65 51 39 DEM25 171 131 106 90 DEM50 205 190 174 169 Table 3 Summary statistics of stream gradients d = 150 d = 300 d = 450 d = 600 d = 750 d = 900 d = 1050 d = 1200 d = 1350 d = 1500 d = 1800 d = 2100 d = 2400 d = 2700 d = 3000 N3864 3828 3814 3778 3764 3728 3714 3678 3664 3628 3578 3528 3478 3428 3378 Mean 0.029 0.028 0.028 0.028 0.028 0.027 0.028 0.027 0.027 0.027 0.027 0.026 0.026 0.026 0.026 Stand. dev 0.036 0.032 0.031 0.030 0.030 0.029 0.029 0.028 0.028 0.027 0.027 0.026 0.025 0.025 0.024 Min −0.078 −0.032 −0.015 −0.013 −0.008 −0.007 −0.006 −0.006 −0.002 −0.002 −0.001 −0.001 0.000 0.000 0.001 Max 0.366 0.298 0.280 0.248 0.230 0.215 0.200 0.196 0.185 0.182 0.172 0.166 0.160 0.151 0.145 Range 0.444 0.330 0.295 0.261 0.238 0.222 0.206 0.202 0.187 0.184 0.173 0.167 0.160 0.151 0.144 Coeff. var 127.2 113.9 110.2 107.6 106.3 104.9 104.0 103.0 102.3 101.4 100.1 98.8 97.4 95.8 94.4 Environmental Earth Sciences (2023) 82:475 1 3 475 Page 8 of 19 0.007m·m−1; max.: 1.083m·m−1), respectively. Largest knickzones (from 3rd quartile of length and height) correspond to the ~ 5% of total knickzones. Smallest knickzones (until 1st quartile of length and height) correspond to the ~ 6%. Considering only one of the two conditions, the largest occupy 45% and the smallest 44%. In the next two subsections, I address in more detail the knickzones abundance and its forms. In the third, the -10 -5 0 5 10 20 15 25 Relative steepness index -1 R (m ) d -6 x10 Threshold for knickzone: -6 R = 5.05·10 d -7 R = -6.57·10 d -5 R = 2.40·10 d 0 16 24 32 Distance from riverhead (km) 8 400 1,000 1,400 Altitude (m) Knickzone B C 2 1 -0.02 0.14 0.00 0.02 0.04 0.06 0.10 0.08 0.12 -1 G (m·m ) d Local gradient Trend gradient 300 900 From 450 to 2,700 G with d (m) of: d 3,000 A y = -2E-05x + 0.1175 3,000300 0.14 0.12 0.10 0.08 0.06 0.04 -1 G (m·m ) d d (m) 1 1 2 y = 7E-07x + 0.0091 0.016 0.014 0.012 0.010 0.008 0.006 -1 G (m·m ) d 3,000 d (m) 300 2 Fig. 4 Example of workflow for the definition of knickzones (Lóuzara River). A represents the stream gradients for all values of d along the Lóuzara River. B, C shows the longitudinal profile and knickzone segments above a threshold of 5.05 × 10–6 Environmental Earth Sciences (2023) 82:475 1 3 Page 9 of 19 475 correlation between the knickzones and different parameters is analyzed. Knickzones abundance Knickzone abundance was analyzed crossing for each central point of the knickzones its value knickzone frequency with four watershed relevant variables organized by classes (elevation, distance from riverhead, gradient of knickzone, and drainage area) (see Fig.5). Knickzone frequency was obtained from GIS through a line density operation. The search radius to estimate the knickzone per square kilometer was 1000m (raster output cell size of 25 × 25m), equivalent to a circumference of ~ 0.3 km2 (circumference area = (0.3/Π)1/2 × 1000). The minimum value of knickzone frequency detected for the study area is 0.102 km−1, and maximum of 2.152 km−1. Mean value is 0.725 km−1 (median 0.708 km−1). Note that this frequency was calculated from the average of the raster pixels from the algebraic operation performed with GIS, and it is different from the frequency calculated directly by dividing knickzone length (km) by watershed area (km2) (0.467 km−1). The GIS procedure was necessary to obtain a particular value for each knickzone (i.e., x, y coordinate of midpoint of the line). Except for the elevation variable, the other three variables shown in Fig.5 reflect a decreasing trend from low values to the highest. This trend is not very intensive, but it is clear. Happened the same with the elevation, there is a soft growing trend. In general terms, apparently there is no distribution pattern and the knickzones are everywhere. 0.50 0.55 0.60 0.65 0.70 0.75 0.80 <400400-600 600-800800-1000 >1000 -1 Knickzone frequency (km ) A Altitude (m) Min.: 270 Median: 725 Mean: 731 Max.: 1,311 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 <0.5 0.5-3 3-6 6-9 9-12 12-15 >15 -1 Knickzone frequency (km ) B Distance from riverhead (km) Min.: 0.151 Median: 0.738 Mean: 1,941 Max.: 35,590 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 <0.2 0.2-0.4 0.4-0.6 0.6-0.8 >0.8 -1 Knickzone frequency (km ) -1 Gradient of knickzone (m·m ) C Min.: 0.007 Median: 0.114 Mean: 0.178 Max.: 1.083 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 <1 1-5 5-10 10-15 15-20 20-25 >25 -1 Knickzone frequency (km ) D 2 Drainage area (km ) Min.: 0.524 Median: 1.249 Mean: 4.186 Max.: 273 Fig. 5 Knickzone frequency for different classes of altitude (A), distance from riverhead (B), gradient of knickzone (C), and drainage area (D) Environmental Earth Sciences (2023) 82:475 1 3 475 Page 16 of 19 A B C Seco Stream waterfall Environmental Earth Sciences (2023) 82:475 1 3 Page 17 of 19 475 enhancement of the knickzones as singular elements could help to slowly reverse this situation for the CMG case. Conclusions Four main conclusions can be drawn from the findings of this research. First, CMG can be considered as the “knickzone land” (i.e., “waterfall land”) from Galicia. Close half of the drainage network was identified as knickzone (47%), which translates into a mean knickzone frequency of 0.467 km−1. In this study, a distribution pattern was not identified and the knickzones are everywhere. Second, the longitudinal profile forms that the river adopts in the form of waterfalls and rapids to overcome the heavy steepness is a direct consequence of its lithological characteristics and tectonic structures. Tectonic controls river trajectory, while lithological control plays a fundamental role as an influential factor in the appearance of characteristic forms and gives “personality” to CMG rivers. Third, knickzones form, analyzed from the height, length, and relative steepness variables, are not governed by a specific spatial distribution. For example, streams order 1 host 70% of the knickzones, and greater length and relative steepness. However, if the height is looked at, the mean value is higher in order 2 than in order 1 knickzones. Fourth, the visual expression in CMG is resolved with a landscape with marked scenic impact. Findings from this study can be used by managers to develop and/or improve strategies for conservation, valorisation, and how to approach the tourist who visits the CMG. Scientific tourism is able to offer a unique and educational travel experience, allowing participants to learn about bedrock rivers and knickzones in the biophysical and cultural context of CMG. This can also contribute to the local economy by supporting science-related businesses and research institutions. Future research on bedrock rivers planned by CMG Scientific Committee looks forward (i) to focus on more detailed morphologies from LiDAR data and topographical fieldwork (e.g., hydraulic geometric); (ii) to analyze the interaction of flow and sediment transport in the self-formed waterfalls (e.g., Scheingross etal. 2019); and (iii) measure the flow characteristics (such as water depth, discharge, specific stream power) to relationship them with forms and processes of the bedrock rivers). With all this, is hoped to know more about the effects of rock properties and environmental factors on the frequency and shape of knickzones. Acknowledgements The author thanks Augusto Pérez Alberti, Ramón Vila “Moncho”, Martín Alemparte, Alejandro Gómez Pazo, and Daniel Ballesteros for her diverse help. Also, to the canyoning colleagues and the reviewers for their comments and suggestions. This work was performed within the framework of the research projects promoted by Courel Mountains UNESCO Global Geopark and the Ribeira SacraCourel LocalActionGroup. Author contributions As it is an article with a single author, it was this author who carried out all the content of the research. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Data availability Not applicable. Declarations Conflict of interest The author has not disclosed any competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. References Alemparte M and Vila R (2018), "Guía de barrancos y cañones Montañas do Courel. Proyecto Geoparque Mundial de la UNESCO", pp. 135. Álvarez-Vázquez M, De Uña-Álvarez E (2021) An exploratory study to test sediments trapped by potholes in Bedrock Rivers as environmental indicators (NW Iberian Massif). Cuaternario y Geomorfología 35(1–2):69–88 Ashmore P (2015) Towards a sociogeomorphology of rivers. Geomorphology 251:149–156 Boothroyd A, McHenry M (2019) Old processes, new movements: the inclusion of geodiversity in biological and ecological discourse. Diversity 11(11):1–17 Brierley G, Fryirs K, Outhet D, Massey C (2002) Application of the River Styles framework as a basis for river management in New South Wales, Australia. Appl Geogr 22(1):91–122 Brierley G, Fryirs K, Cullum C, Tadaki M, Huang HQ, Blue B (2013) Reading the landscape: Integrating the theory and practice of Fig. 10 Diverse images about the bedrock rivers characteristics from CMG. Canyoning map to downfall the Ferreiriño Stream (Quiroga River tributary) and Seco Stream (Ferreiriño Stream tributary); courtesy from Alemparte and Vila (2018) (A). Quarrying erosion processes with different detail scale, showing how form and processes are the same, only change the zoom (B). Two examples of waterfalls(©Ramón Vila) (C): O Pombar waterfall (20m height) in Quiroga River (left), where a quartzite band offer major resistance to the erosion than slate or schist, and waterfall in Seco Stream with ~ 35m height (right), where the erosive action triggered a setback in this reach embedded in quartzite and slate (see image A to locate the waterfall) ◂ Environmental Earth Sciences (2023) 82:475 1 3 475 Page 18 of 19 geomorphology to develop place-based understandings of river systems. Progress Phys Geogr Earth Environ 37(5):601–621 Bierman PR, Montgomery DR (2013) Key concepts in geomorphology, p. 525. https:// www. amazon. com/ KeyConce ptsGeomo rphol ogyPaulBierm an/ dp/ 14292 38607 Byun J, Paik K (2021) The development process of the Korean coastal mountain range: examination from spatial distribution of knickzones. Progress Phys Geogr Earth Environ 45(4):541–563 Cantreul V, Bielders C, Calsamiglia A, Degré A (2018) How pixel size affects a sediment connectivity index in central Belgium. Earth Surf Proc Land 43(4):884–893 Chen A, Lu Y, Ng YCY (2015) Role of Tourism Earth-science in Tourism Development. In: The principles of geotourism. pp. 25–38. Springer Berlin Heidelberg, Berlin Church M, Ferguson RI (2015) Morphodynamics: Rivers beyond steady state. Water Resour Res 51(4):1883–1897 Comas X, Wright W, Hynek SA, Fletcher RC, Brantley SL (2019) Understanding fracture distribution and its relation to knickpoint evolution in the Rio Icacos watershed (Luquillo Critical Zone Observatory, Puerto Rico) using landscape-scale hydrogeophysics. Earth Surf Proc Land 44(4):877–885 Crofts R (2014) Promoting geodiversity: learning lessons from biodiversity. Proc Geol Assoc 125(3):263–266 De Uña-Álvarez E, Vidal-Romaní JR, Rodríguez R (2009) Erosive forms in river systems. In: Advances in studies on desertification. Murcia, pp. 465–468 De Uña-Álvarez E, Álvarez-Vázquez M and Rodríguez R (2014), "Tipología de formas graníticas en el tramo medio del río Miño (Ourense, Galicia, NW del Macizo Ibérico)", In Avances de la Geomorfología en España 2012–2014. Cáceres, pp. 434–437. de Vicente G, Vegas R (2009) Large-scale distributed deformation controlled topography along the western Africa-Eurasia limit: Tectonic constraints. Tectonophysics 474(1–2):124–143 DiBiase RA, Rossi MW, Neely AB (2018) Fracture density and grain size controls on the relief structure of bedrock landscapes. Geology 46(5):399–402 García H (2014) Geomorfología fluvial en sistemas atlánticos: metodología de caracterización, clasificación y restauración para los ríos de Galicia. Thesis at: University of Santiago de Compostela. Santiago de Compostela García H, Pérez-Alberti A (2021) Aproximación a la identificación y caracterización de ríos en roca a escala regional mediante variables topo-geomorfológicas (Galicia, Noroeste de la Península Ibérica). Cuadernos De Geografía 107:217–242 García JH, Ollero A, Ibisate A, Fuller IC, Death RG, Piégay H (2021) Promoting fluvial geomorphology to “live with rivers” in the Anthropocene Era. Geomorphology 380:107649 Garzón G, Ortega JA and Garrote J (2008), "Morfología de perfiles de ríos en roca. Control tectónico y significado evolutivo en el Bajo Guadiana", Geogaceta, pp. 63–66. Garzón-Heydt G (2010), "Significado de los ríos en roca y la importancia de su preservación", In Patrimonio geológico: los ríos en roca de la Península Ibérica. Madrid, pp. 17–36. Gray M (2008) Geodiversity: developing the paradigm. Proceedings of the Geologists' Association. Vol. 119, pp. 287–298 Hankock GS, Anderson RS, Whipple KX (1998) Beyond power: bedrock river incision process and form. In: Rivers over rock: fluvial processes in bedrock channels, Vol. 107, pp. 35–60. American Geophysical Union, Washington Hayakawa YS, Oguchi T (2006) DEM-based identification of fluvial knickzones and its application to Japanese mountain rivers. Geomorphology 78(1):90–106 Hayakawa YS, Oguchi T (2009) GIS analysis of fluvial knickzone distribution in Japanese mountain watersheds. Geomorphology 111(1):27–37 Henriques MH, dos Reis RP, Brilha J, Mota T (2011) Geoconservation as an Emerging Geoscience. Geoheritage 3(2):117–128 Horacio J, Montgomery D, Ollero A, Ibisate A, Pérez-Alberti A (2018) Application of lithotopo units for automatic classification of rivers: concept, development and validation. Ecol Indicators 84:459–469 Horacio J, Muñoz-Narciso E, Sierra-Pernas JM, Canosa F, PérezAlberti A (2019) Geo-Singularity of the Valley-Fault of Teixidelo and Candidacy to Geopark of Cape Ortegal (NW Iberian Peninsula): Preliminary Assessment of Challenges and Perspectives. Geoheritage 11(3):1043–1056 Hose TA, Marković SB, Komac B, Zorn M (2011) Geotourism—a short introduction. Acta Geogr Slov 51(2):339–342 Ibisate A, Ollero A, Díaz E (2011) Influence of catchment processes on fluvial morphology and river habitats. Limnetica 30(2):169–182 Lima AG, Binda AL (2013) Lithologic and structural controls on fluvial knickzones in basalts of the Paraná Basin, Brazil. J South Am Earth Sci 48:262–270 May C, Roering J, Snow K, Griswold K, Gresswell R (2017) The waterfall paradox: How knickpoints disconnect hillslope and channel processes, isolating salmonid populations in ideal habitats. Geomorphology 277:228–236 Mikhailenko AV, Nazarenko OV, Ruban DA, Zayats PP (2017) Aesthetics-based classification of geological structures in outcrops for geotourism purposes: a tentative proposal. Geologos 23(1):45–52 Miller VC (1953) A quantitative geomorphic study of drainage-basin characteristics in the Clinch Mountain area, Virginia and Tennessee. Thesis at: Off. Nav. Res. (U.S.), Geogr. Branch Montgomery DR (1996) Influence of geological processes on ecological systems. In: The Rain Forests of Home: Portrait of a North American Bioregion, pp. 43–68. IslandPress Muthusamy M, Casado MR, Butler D, Leinster P (2021) Understanding the effects of Digital Elevation Model resolution in urban fluvial flood modelling. J Hydrol 596:126088 Ollero A (2017), "Hidrogeomorfología y geodiversidad: el patrimonio fluvial" Ayuntamiento de Zaragoza. Ortega JA (2007) El estudio de la morfología de los ríos en roca. Implicaciones hidrológicas y evolutivas en dos barrancos españoles. Bol Geol Min 118(4):803–812 Ortega JA (2010), "Morfologías en los ríos en roca. Variaciones y tipologías", In Patrimonio geológico: los ríos en roca de la Península Ibérica. Madrid, pp. 44–78. Ortega JA and Durán JJ (2010), "Patrimonio geológico: los ríos en roca de la Península Ibérica" Madrid, pp. 504. Instituto Geológico y Minero de España. Ortega-Becerril JA, Jorge-Coronado A, Garzón G, Wohl EE (2017) Sobrarbe Geopark: an Example of Highly Diverse Bedrock Rivers. Geoheritage 9:533–548 Ortega-Becerril JA, Polo I, Belmonte A (2019) Waterfalls as Geological Value for Geotourism: The Case of Ordesa and Monte Perdido National Park. Geoheritage 11:1199–1219 Pérez-Alberti A (1982), "Hidrografía", In Xeografía de Galicia. Vol. Tomo I - O Medio, pp. 97–110. Ed. Sálvora. Pérez-Alberti A (1985), "Un exemplo de estudio integral do medio:a cunca do río Miño", I Cuaderno de Xeografía. Vol. SociedadeGalega de Xeografía, pp. 11–32. Pérez-Alberti A (1993), "Xeografía de Galicia: Xeomorfoloxía" Santiago de Compostela Vol. 3 Xunta de Galicia. Pérez-Alberti A (2000), "La red fluvial", In Geografía de Galicia. , pp. 193–240. Faro de Vigo. Ruban DA (2016) Representation of geologic time in the global geopark network: a web-page study. Tourism Manag Perspect. 20:204–208 Scheingross JS, Lamb MP, Fuller BM (2019) Self-formed bedrock waterfalls. Nature 567:229–241 Environmental Earth Sciences (2023) 82:475 1 3 Page 19 of 19 475 Scott DN, Wohl EE (2019) Bedrock fracture influences on geomorphic process and form across process domains and scales. Earth Surf Proc Land 44(1):27–45 Selby MJ (1993) Hillslope materials and processes. Oxford Sharma A, Tiwari KN, Bhadoria PBS (2011) Determining the optimum cell size of digital elevation model for hydrologic application. J Earth Syst Sci 120:573–582 Simon A, Castro J and Rinaldi M (2016) Channel form and adjustment: characterization, measurement, interpretation and analysis. In: Tools in fluvial geomorphology, pp. 237–259. Wiley Blackwell Tarboton DG, Bras RL, Rodriguez-Iturbe I (1991) On the extraction of channel networks from digital elevation data. Hydrol Process 5:81–100 Tinkler KJ (2004) Knickpoint. In: Encyclopedia of geomorphology. pp. 595–596. Routledge Ltd. Taylor & Francis Group, London and New York Wohl EE (1998) Bedrock channel morphology in relation to erosional processs. In: River over rock: fluvial processes in bedrochk channels, vol. 107, pp. 133–151. American Geophysical Union Wohl E (2015) Particle dynamics: The continuum of bedrock to alluvial river segments. Geomorphology 241:192–208 Wright M (2020) “A typology of hydraulic barriers to salmon migration in a bedrock river. Thesis at: Simon Fraser University. Canada Zahra T, Paudel U, Hayakawa YS, Oguchi T (2017) Knickzone extraction tool (KET)—a new ArcGIS toolset for automatic extraction of knickzones from a DEM based on multi-scale stream gradients. Open Geosci 9(1):73–88 Zgłobicki W, Baran-Zgłobicka B (2013) Geomorphological Heritage as a Tourist Attraction. A Case Study in Lubelskie Province, SE Poland. Geoheritage 5(2):137–149 Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.