How Can Be Lotic Ecosystem Size More Precisely Estimated? Comparing Different Approximations in Pre-Pyrenean and Pyrenean Mountains

Rivers are among the most biodiverse and endangered ecosystems on earth. In Europe, concern over their conservation promoted the development of legal instruments for habitat and species conservation, the Habitats Directive, and water resource management, the Water Framework Directive. This legal protection demanded the estimate of river ecosystem surface for different purposes. Different approaches allow river surface to be measured at a low cost. Some accurate techniques like satellite images or LiDAR (Light Detection and Ranging) do not always work at a large scale or for streams and small rivers. We discuss here the use of the traditional hydraulics relationship between drainage area and bankfull width as a good approach to river surface estimation. We confirm that the use of this cheap and simple method could be a good approach to estimate river surface. However, we also proved that the development of regional curves, i.e., to establish the empirical relationship based on study area data, constitutes an essential improvement to estimation.

Rivers support some of the most biodiverse ecosystems in the world and provide essential ecosystem services to society [10], but, at the same time, freshwater ecosystems may be the most endangered ecosystems in the world [6], and rivers are between the most threatened [11,12] due to overexploitation, pollution, regulation, climate change, drainage-basin degradation, non-native species, and synergistic impacts [6,13]. As a consequence, this extensive habitat deterioration is causing a decline of biodiversity in freshwater ecosystems far greater than in the most affected terrestrial ecosystems [14].
Flow and its regime are considered as the main structuring factors of stream communities [15,16], and several studies have demonstrated a relationship between the hydrological regime and biological communities [17][18][19][20][21]. However, it is estimated that around 50% of the river volume is currently altered by either flow regulation or fragmentation, and it is expected to increase to 93% with the construction of new major hydropower dams [22].
Expanding hydropower is considered as one of the 12 pressing and emerging threats to freshwater biodiversity [11]. In Europe, most rivers are largely impacted and significantly modified by anthropogenic structures like dams, weirs, and barrages [23], negatively impacting the ecological processes and aquatic species [24]. The rivers also suffer from fragmentation by barriers to free flow [7,25], causing loss of habitat, altering the natural flow regimen, and reducing the connection between the main channel and floodplain due to flow regulation [23]. Hydrological variation caused by weirs and barrages is a substantial contributor to aquatic habitat alteration and community disturbance by altering the longitudinal continuity of lotic ecosystems, the lateral interaction with floodplains, the interchange with adjacent groundwaters, and the natural variation in their relative magnitude over time [26][27][28]. As a result, European rivers suffer from intense morphological habitat degradation, which is considered one of the most important river ecosystem alterations in Europe [29]. Pyrenean rivers are not an exception, as they are highly affected by reservoirs that impound about one-third of the mean annual runoff [30].
In this context, our societies have called for increased conservation of our ecosystems. Responding to these requirements, the Habitats Directive (HD), the main legislative work of the European Union's nature conservation policy, and the Water Framework Directive (WFD), the main legislative work regarding water resource management, have been developed to conserve our ecosystems. Both directives, from different perspectives and approximations, have common objectives and links [31]. These legal responsibilities demand the estimation of river size and its evolution at different scales. On one hand, the Habitats Directive requires surface estimation for the different running water habitats every six years [32,33]. On the other hand, the Water Framework Directive demands the size estimation of each WFD management unit [34], and by evaluating the hydromorphology of ecosystems the connectivity and continuity of these ecosystems are indirectly measured. Therefore, a method for river size estimation for all streams and rivers at a large scale is mandatory in order to achieve the legal requirements. However, the lack of a common method and some interpretation and implementation problems have led to calls for an established protocolized method for monitoring running water habitats.
River size is a key parameter related to several river ecosystem processes such as flux exchange with the atmosphere [35,36], species richness [37], and nutrient removal [38]. Many different measures and approximations have been used to estimate river size [39], e.g., length [2,37,40], drainage area [36,41], and river surface. At a large scale, easy measures like drainage area or length have been traditionally widely used [39]; however, new techniques like satellite images can be used to measure river surface [35,42]. As of yet, these estimation methods are limited by pixel resolution, allowing measurements for rivers wider than 30 m [35,42]. Other techniques, like LiDAR and photogrammetry, are limited at a large scale by their computing processing time [43,44]. Thus, several methods have been proposed to estimate river width and river surface [35,36], but none of them seems to be widely applicable with satisfactory results and with an optimal cost-benefit ratio.
When the river ecosystem surface is estimated, it is necessary to delimitate its lateral extent, which can refer to the water channel, the bankfull, the riparian zone, or even the floodplain. Estimation of the bankfull and water channel at a large scale has been widely developed [35,36,42,45]; however, riverine and floodplain surface estimation has only been developed at a local scale [46,47]. In this study, we considered the bankfull width as the lateral extent for the river ecosystem assessment due to two main reasons: (1) Ecosystem classifications usually differ riparian and floodplain ecosystem from running water ecosystems, e.g., habitat types of the HD [48]; (2) although the river channel and floodplain could be considered as components of a single dynamic system [27], the main material and energy of most European rivers are derived through local production and the riparian zone [49]. Therefore, bankfull width is indeed a good proxy for river lateral extent.
Bankfull width has been measured by different methods in stations or small sections including field measurements [50,51], high-resolution digital elevation models, and photogrammetry [43,52]. However, these measures only have been realized in stations or small sections. At a large scale, the most widely accepted estimation for bankfull channel geometry relates river width to discharge or drainage area in a power law model [42,53]. These power laws have been largely estimated in the USA at a regional scale and for the entire country [54][55][56]. However, few studies have estimated them in Europe and only at a regional or local scale [51]. A global-scale bankfull estimation method was recently developed but only for rivers wider than 30 m [42].
Therefore, the main aim of this study was to analyze the size of freshwater lotic ecosystems from different approximations and methodologies to determine which of the existing methodologies are the most accurate and adequate by identifying the advantages and disadvantages of these approximations and, as a consequence, quantifying the size of the freshwater lotic ecosystem to achieve the mandatory requirements of the HD and WFD by providing a proved tool that can be adapted to local contexts, enabling easy implementation.

Study Area and Data
This study was developed for Pre-Pyrenean and Pyrenean rivers, an ecosystem type defined in the methodologies for the monitoring of the conservation of habitat types in Spain [57]. This ecosystem type is defined based on physical, chemical, and biological variables and covers rivers along approximately 390 km of the southern slopes of the Pyrenees [57] (Figure 1). They are small rivers and streams located in altitudes between 700 and 900 m with a mean discharge between 1.37 and 4.93 m 3 · s −1 and a drainage area between 4.5 and 3800 km 2 . A total of 232 management water bodies defined according to the WFD are classified as this type of ecosystem. Bankfull widths for 76 river sections, or reaches, of 1-2 km were collected from the national cartography of legal boundaries delimitation [58]. This digital maps establish the bankfull channel limits based on historical and actual orthophotos, digital elevation models, and bankfull discharge [59]. However, they have only been developed for a few rivers. For each of the 76 water bodies where these guidelines have been established, we randomly selected a point and created a buffer of 1 km to delimitate and standardize a river section. For every section, we obtained the bankfull limits via the national guidelines of legal river boundaries ( Figure 1).

Methods: Size Estimation
The surface of each polygon delimitated by the bankfull channel banks was calculated ( Figure 1). This surface was considered as the real bankfull channel surface of each section (Sb). We extracted the centerline of each section using a skeletonization procedure of the bankfull channel bank lines. The length (L) of each centerline was calculated as an approximation of the river ecosystem size, a classical and probably the most widely used approximation of river size. The bankfull width of each section (Wb) was calculated as the median of the measurements of Wb every 50 m within each reach. The drainage area of each section (A) was calculated as the mean value of all sections. We used the drainage area raster developed by the Spanish government for its estimation [60]. All these operations were conducted using ArcMap 10.5 [61].
First, we estimated the power law that relates bankfull width and drainage area for our study area. We used the median value of the bankfull width (Wb) and the mean drainage area (A) of each section to estimate the α' and β' coefficients. To develop the power law for our study, we log-transformed both variables and estimated α' and β' through simple linear regression, i.e., Y = β 0 + β 1 X, where Y = log 10 (Wb), X = log 10 (A), and β 0 and β 1 are regression parameters (intercept and slope, respectively). Thus, the coefficient and the exponent of the power law relation (Equation (2)) are α' = 10 β0 and β' = β 1 , respectively.
To explore the relationship between bankfull channel width and drainage area, we used multiple linear regression models in combination with mean annual precipitation, channel slope (primary controls on discharge and channel morphology), and elevation as explanatory variables of channel width [56]. Mean annual precipitation in mm (MAP) data were collected from WordlCilm2 [74]. Mean altitude (MA) in m and mean channel slope (MS) were obtained using a digital elevation model of 5 m pixels. All variables, except elevation, were log-transformed, yielding a multiplicative effects model of the form [58]. We ranked the models using the corrected Akaike information criterion (AICc) and selected the model with the lowest AICc (∆AICc < 2).
To estimate the bankfull channel surface for each water body, the average bankfull width (Wbe) was calculated for each section using our regional power law, considering the mean drainage area of the section (A). Then, we estimated the bankfull surface (Sbe) for a water body as the product of the centerline length (L) and Wbe, considering the section as a parallelogram. This operation was repeated for all the power laws found in the bibliography to estimate bankfull surface (Wbe) for each one.

Methods: Statistical Analysis
We used multiple linear regression to explore the differences between the real bankfull surface (Sb) of each section and the estimated surface (Sbe) according to our regional power law. We included the length of the centerline (L) and median width (Wb) of the section, as well as the sinuosity index (SI) [2] (primary controls on size, width and length, and channel complexity, sinuosity index). We ranked the models using the corrected Akaike information criterion (AICc) and selected the model with the lowest AICc (∆AICc < 2). All data and results are shown in Supplementary Materials 1.
Finally, we compared the differences between the real bankfull surface (Sb) and the estimated surfaces (Sbe) via linear regressions. The Sb and Sbe values were previously logtransformed to achieve residual normality, linearity, and homoscedasticity. The resultant linear models had the expression log 10 (Sb) = m + n log 10 (Sbe), where Sb and Sbe are the real and the estimated bankfull surface, and m and n are the model parameters (intercept and slope, respectively). If we remove the logarithms, the resulting expression is Sb = 10 m Sbe n . Thus, when m is nearer to 0, and n to 1, the estimated values of Sb are more similar to the real Sbe.
All statistical analyses were carried out by using R version 4.0.2 [75].

Results
The relationship between the logarithms of Wb and A showed a significant increase in log 10 (Wb) with the increase of log 10 (A) (Figure 2) (F 1,74 = 62.89, p < 0.001, r 2 adjusted = 0.45). With the linear model coefficients, we calculated α and β values of Equation (2) for our regional data. The regional power law relationship is as follows: The multiple linear regression including in the relationship log 10 (MAP), log 10 (MS), and MA showed an improvement in the relation between log 10 (Wb) and log 10 (A) when log 10 (MS) was included in the model (Table 1), although the model with only log 10 (A) as an explanatory variable was roughly equivalent to the best model, and it was the only significant explanatory variable for all models (Table 1). Table 1. Comparison of Akaike information criterion (AICc) deltas and weights for linear models testing the effects of drainage area (A), mean channel slope (MS), mean annual precipitation (MP), and mean altitude (MA) and response (bankfull Width, Wb). Shown are the number of parameters (K), the AICc, the ∆AICc, AICc weight (wAICc), and the evidence ratios (ERs). The best model and models roughly equivalent to the best (∆AICc ≈ 2 or lower) are highlighted in bold. The only significant explanatory variable in all models is log 10 (A) (p < 0.005). The sample size is 76 for all models. The multiple linear regressions developed to assess the differences between the estimated bankfull channel surface (Sbe) and real bankfull surface (Sb) showed that size and river complexity influence the estimation error. As shown in Table 2, the most parsimonious model only included Wb and SI (F 2,73 = 94.35, p < 0.001, r 2 adjusted = 0.72). The difference between Sb and Sbe was positively correlated with Wb and SI. The surface was overestimated (Sbe > Sb or Sb − Sbe < 0) in sections with small Wb and underestimated (Sbe < Sb or Sb − Sbe > 0) where Wb was high (Figure 3a). SI showed a slightly positive correlation with the estimation error. The surface of complex rivers (SI > 1.5) was more frequently overestimated than the surface of less complex rivers (SI ≈ 1) (Figure 3b). Table 2. Comparison of AICc deltas and weights for linear models testing the effects of bankfull width (Wb), sinuosity index (SI), river length (L), and bankfull channel Surface (Sb) and response (differences between real and estimated bankfull channel surface, Sb -Sbe). Shown are the number of parameters (K), the AICc, the ∆AICc, AICc weight (wAICc), and the evidence ratios (ERs). The best model and models roughly equivalent to the best (∆AICc ≈ 2 or lower) are highlighted in bold. Lastly, we compared the results of the bankfull channel surface estimation obtained by our regional curve with the results obtained by the other power laws. Figure 4 shows the regional curve as one of the most similar to the real values. Approximately 50% of the linear models between log 10 (Sb) and log 10 (Sbe) had a higher r 2 adjusted than the model in which Sbe was estimated by our regional curve (F 1,74 = 51.24, p < 0.001, r 2 adjusted = 0.40) (the highest r 2 adjusted is 0.41). However, the mean differences between Sb and Sbe estimated by the regional curve were the third-smallest.

Discussion
As expected, and confirming previous similar studies, the drainage area is intensively related to width, and our results have shown that the estimated bankfull channel surface is directly related to the real bankfull surface, and this estimation significantly improves when the mean channel slope is included in the model. Moreover, the development of a regional power law based on real data improved the estimation. However, the surface was overestimated and underestimated due to river size and complexity. On one hand, upstream sections were normally overestimated, whereas downstream sections were underestimated (Figure 3a). On the other hand, sections highly sinuous were more frequently underestimated (Figure 3a).
We demonstrate the utility of power laws relating drainage area and bankfull width as a fast assessment of river ecosystem size for streams and small rivers at a large scale, confirming this method as a tool of easy implementation. Other studies indicated the importance of developing power law relationships between bankfull channel width and drainage areas for preliminary channel design for restoration projects [73] or modeling habitat availability and suitability for fish [76]. According to our results, these estimations are robust and can be consolidated as an optimal approximation of river size with a low costbenefit ratio, at least until LiDAR and remote sensing can be optimized for these purposes. We found that the best obtained results corresponded to the power law established with our data within the regional context of the Pre-Pyrenean and Pyrenean range. Thanks to these results, we can state that using regional power law models is the most optimal approximation, and in case comparing with other regions or models is required, models can be easily compared by recalculating the values. As outlined, using a regional power law provides the best results. However, in case it is not possible to use a regional power law due to the lack of digital terrain models or real data, it is also possible to use a power law model from similar ecoregions.
Our reported value for α could be considered as high compared with the values found in the bibliography, the 13th largest α detected/historically detected, whereas the β value is near the median value of all β values registered, 0.36 (see Data Availability Statement). High values of α, such as those estimated in our study area, are related to high pluviometry [56]. The mean pluviometry for all ecosystem sizes is 1273.84 ± 213.84 [57], a value that could be considered as relatively high and could be related to the fact that the Pyrenees is a mountain range that links the Atlantic Ocean with the Mediterranean Sea with a clear west-east gradient of humidity and pluviometry. The β value is lower the less the width grows with the watershed area, but this relationship can be affected by the degree of alteration of the river basin. It has been proved that human impacts can reduce β values [56,69,70]. On the other hand, high β values have also been related to climatic conditions, e.g., arid zones [56]. In our case, the high estimated β value is coherent with other studies [56,70] because most of our studied rivers can be considered as altered by dams, hydropower stations, etc. [30]. Moreover, the hydrological alteration and flow regulation caused by these perturbances can influence the estimated regional power law, increasing the β value [70].
To improve our model, we included in the power law other factors that are related to discharge (mean annual precipitation) and channel morphology (mean channel slope and mean altitude). We found only slight improvement when we included mean channel slope in the model (∆AICc < 2). This result is coherent since the chosen ecosystem type was defined partly based on physical variables such as discharge, altitude, or slope [57]. Therefore, it is logical that the relationship was not significantly improved based on these variables because they are intrinsically included.
As previously mentioned, the regional power law provides the most optimal results. However, this is an imperfect approximation, and the comparison between the real and the estimated bankfull channel area (Figure 3) shows an underestimation in wider sections that are located in downstream reaches and an overestimation in small lotic ecosystems ( Figure 3). Degradation of the riverbank often causes a widening processes [77], which could explain the observed underestimation, i.e., the downstream reaches are wider than expected. Most of the areas surrounding rivers in the low parts of our study area are occupied by cultures, whereas the high areas are mostly forest in a mountainous area where rivers flow through narrow valleys where widening is not possible and overestimation of Wbe is explained. According to these results, altitudinal gradient and orography could explain these over-and underestimation results. By including mean channel slope (MS) in our model, we obtained a good model (ER = 1) for explaining river width (Table 2). This is coherent with the mentioned under-and overestimation problems because upstream rivers have a higher MS than lower reaches. The surface was slightly underestimated in more complex sections with higher SI. Complex rivers had a higher heterogeneity of width [78], which was not captured in the estimations producing underestimations.
As is shown in Figure 4, the variation in the empirical values of the power laws triggered a huge variability in surface estimation. This variation has been reported in several studies developed in the USA [54][55][56], and it has been demonstrated that different factors influence this relationship, e.g., human alteration [69,70], riparian vegetation [69], riverbed material [56], and climate and topographic factors [56]. Therefore, it is important to develop a power law at a regional scale for river size estimation when possible. However, all estimated bankfull channel surfaces are significantly related to the real surface and could be useful for comparison of the river size.
These results confirm that regional power law models are superior to other models that have been developed for specific hydrogeological contexts. Of course, even in a regional context, we can expect variability and errors in the estimation as we have observed. For this reason, the use of regional power law models is conditioned by the characteristics of the region; therefore, these regions must be internally homogenous with similar climate, geology, pluviometry, and, of course, orography. However, this fast approach to river size might be replaced in the future by more accurate new techniques such as LiDAR or satellite images that currently present limitations [44]. Indeed, recent studies have demonstrated the precision of satellite images to estimate bankfull and water channel at a large scale for rivers wider than 30 m [35,42] or accurate measurements derived from LiDAR [79]. However, these methods still present limitations that restrain their use at a large scale for narrower rivers or sections. As a consequence, using regional power laws could be considered as an optimal temporary solution until these more accurate techniques can be applied to most of the lotic ecosystems. With the application of this method, we can address the lack of surface data from lotic ecosystems since the Habitats Directive was approved 30 years ago. Classical methods like river length also should be calculated because of historical reasons and due to the rivers mainly comprising a lineal ecosystem where the most important processes and patterns take place in longitudinal axes [80][81][82]. Acknowledgments: This manuscript would not be possible without the support and the motivation of Rafael Hidalgo, Olga Lamas Murúa, and Elena Bermejo Bermejo. This manuscript was born within the context of the Project "Establishment of a national monitoring system of the Conservation Status of the Habitat of Community Interest" 3081147 and supported by the Government of Spain and the Ministry for the Ecological Transition and Demographic Challenge.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.