Resistance Analysis of Morphologies in Headwater Mountain Streams

: River ﬂow velocity is determined by the energy available for ﬂow motion and the energy fraction lost by ﬂow resistance. We compared the performance of different equations for the Darcy-Weisbach resistance coefﬁcient ( f ) and empirical equations to predict ﬂow velocity. The set of equations was tested using data from the Quinuas headwater mountain river in the Andean region. The data was collected in three Cascades , two Step-pools , and one Plane-bed covering a wide range of velocity magnitudes. The results reveal that nondimensional hydraulic geometry equations (NDHG) with a Nash-Sutcliffe efﬁciency index (EF) varying from 0.6–0.85 provide the most accurate velocity prediction. Furthermore, the study proposes a methodology applicable to all morphologies for deﬁning the NDHG parameters using easily measured ﬁeld data. The results show an improvement in predictability with EF values in the range of 0.81–0.86. Moreover, the methodology was tested against data from the literature, which was not divided into morphologies providing EF values of around 0.9. The authors encourage the application of the presented methodology to other reaches to obtain additional data about the NDHG parameters. Our ﬁndings suggest that those parameters could be related to reach characteristics (e.g., certain characteristic grain size), and in that case, the methodology could be useful in ungauged streams.


Introduction
Prediction of the mean velocity in a river is important from a scientific and practical point of view. Nondimensional hydraulic geometry equations (NDHG) are capable of directly estimating the mean velocity but have parameters that vary according to the river morphology [1]. Indirect mean velocity estimation using the Manning, Chézy, or Darcy-Weisbach equations is also common. Predictive empirical equations (PEEs) are focused mainly on the prediction of the Darcy-Weisbach (f ) dimensionless resistance parameter, which has a physical meaning. However, PEEs face several challenges in mountain rivers. First, each PEE is derived under different flow conditions and river morphologies, using different measuring techniques [2]. Second, mountain river characteristics such as steep slopes (bed slope (S 0 ) greater than 0.2%) [3], an average depth comparable to bed material size [4], and a coarse, poorly sorted clast [5] result in resistance patterns that differ from plane rivers. Third, mountain river morphologies such as Step-pools [6][7][8][9], Cascades [10,11], and Plane-beds [8] add complexity to resistance analysis as each morphology possesses different resistance characteristics. Authors have estimated PEE uncertainties of 30% [5] and 66% [8] in mountain rivers. Hence, an analysis of PEE for mountain streams under different flow, geographical, and dissipative conditions is needed to improve velocity estimation approaches.
For this study, experimental data were collected from a 1500 m longitudinal section in the headwaters of the Quinuas river (Figure 1), located in the province of Azuay, Ecuador. The river section is situated between 0 + 000 at 3664.4 masl and 1 + 431.13 at 3605.77 masl, with a mean slope of 4%. This reach was selected given the relatively large variation in river morphology consisting of three Cascades, two Step-pools, and one Plane-bed. Table 1 presents pictures and schemes of each morphology to illustrate its bed characteristics. The geometric characteristics of each morphology are presented in Table 2. The reach length was measured along the thalweg; the bed slope (S 0 ) for each morphology was obtained through linear regression of thalweg points.
The river cross-section (XS) in the reaches without abundant vegetation was measured with a differential GPS Trimble ® R6 instrument (Sunnyvale, CA 94085, USA); in reaches with abundant trees, a total station Sokkia ® 550 RX (Olathe, KS 66061, USA) was used. The average geometry of each reach was calculated from three XSs, as shown in Table 2, except for Cascade 3, where the reach was divided into five XSs. The wetted width (w) at each cross-section was estimated with a measuring tape, excluding the width of the boulders stretching above the water level. The average depth (d) at each XS was computed using the continuity equation assuming a rectangular XS (Equation (1)). These values were averaged per reach to yield a weighted average water depth [2].
where Q is the discharge (m 3 /s), U is the mean flow velocity (m/s), A is the XS area below the water surface (m 2 ), w is the wetted width at the water surface (m), and d is the mean water depth (m).
Water 2021, 13, x FOR PEER REVIEW 3 of 17

River Reach
For this study, experimental data were collected from a 1500 m longitudinal section in the headwaters of the Quinuas river (Figure 1), located in the province of Azuay, Ecuador. The river section is situated between 0 + 000 at 3664.4 masl and 1 + 431.13 at 3605.77 masl, with a mean slope of 4%. This reach was selected given the relatively large variation in river morphology consisting of three Cascades, two Step-pools, and one Plane-bed. Table  1 presents pictures and schemes of each morphology to illustrate its bed characteristics. The geometric characteristics of each morphology are presented in Table 2. The reach length was measured along the thalweg; the bed slope (S0) for each morphology was obtained through linear regression of thalweg points.

Cascade
Step-pool

Cascade
Step-pool

Cascade
Step-pool

Cascade
Step-pool Step-pool

Cascade
Step-pool

Cascade
Step-pool The roughness parameter was estimated using field measurements of discharge, velocity, and energy slope. The discharge was measured using the dilution-gauging method because in small streams, especially in low-flow conditions, measuring flow using the standard wading rod method is difficult and inaccurate. The dilution-gauging method is based on measuring the dilution of a known volume of conservative salt tracer [21]. NaCl was used for its low cost and wide use in small river studies [7,9,22,23]. Velocity was evaluated through the reach length and the subtraction of time-of-travel from conductance curves read upstream and downstream of the reach. Two HOBO U24-001 ® (U24-001 Bourne, MA 02532, U.S) freshwater conductivity data loggers with a resolution time of 1.0 s were placed upstream and downstream of the reach. The traveling time of each conductance curve was calculated using the harmonic method [11]. Energy slope (S F ) was estimated as the water surface slope (S W ) [8,24]. Pebbles were counted to estimate bed material distribution [25]. The number of sampled elements was 400 for each reach. Data were taken at flow magnitudes from 0.03 m 3 /s to almost 1 m 3 /s. During that flow range, there was in-bank flow. However, only for the highest flow reaches such as Plane-bed and Step-pool received bank full flow. The Cascades did not reach bank full flow because of its high slope. Table 2 depicts the range of Manning's coefficient measured at each studied reach. Hence, the data range presents a complete overview of different velocity magnitudes, which are closely related to the resistance characteristics. The roughness parameter used in this study is the Darcy-Weisbach resistance coefficient (f ); most literature from the last three decades has used this parameter given its physical interpretation and dimensionless units [4]. f is estimated with Equation (2).
where g is the gravitational acceleration (m/s 2 ), R H represents the hydraulic radius (m), S F is the energy slope, and U is the velocity (m/s).

Empirical Resistance Equations
The equations tested for velocity prediction are listed in Table 3. Most of the equations estimate the resistance parameter as (8/f ) 1/2 . The velocity was computed using Equation (2). Estimations of NDHG equations derive flow velocity through the following steps: (1) Estimation of dimensionless unitary flow; (2) estimation of dimensionless velocity with the NDHG equation; (3) velocity estimation with the dimensionless unitary velocity definition.

Statistical Performance Metrics
Generally, statistical metrics provide information about a single aspect or projection of the model error. Thus, it is advisable to use a combination of metrics to assess the overall model performance [26]. In this study, six metrics were selected: the root mean square error (RMSE), a qualitative methodology in which larger model errors have more weight than smaller ones [26,27]; the logarithm of RMSE to determine the model prediction capacity for low values [4]; the prediction errors (PE) counting the number of predicted values that are greater than twice or less than half as large as the observed values [4]; the average standard error of estimation (S X ), which gives the mean percentage of the model error relative to the observations [19]; the Nash-Sutcliffe efficiency index (EF), a metric widely used to determine the model goodness-of-fit with flexibility and reliability [28][29][30]; the mean average error (MAE) considered by Willmott and Matsuura [27] is a better indicator of the average error than RMSE. Ritter and Muñoz-Carpena [31] provided a table with a range of EF values, a useful tool for interpreting the score of the model goodness-of-fit.
RMSE and MAE are transformed into a relative version for some analyses. Hence, these two metrics are divided by the average of the observed values and multiplied by 100 to obtain a percentage.

Determination of NDHG Parameters
In this research, a methodology to determine the NDHG equation parameters is proposed. The Rickenmann and Recking [12] dimensionless velocity (Equation (4)) and unitary flow (Equation (5)) obtained by Nitsche et al. [11] in the dimensional analysis are used in this process. The form presented in Equation (3) was selected for the methodology development.
where q (m 2 /s) is the unitary flow (q = Q/w) and D 84 (m) corresponds to the 84th percentile of the grain-size distribution. Equation (3) was linearized through the application of logarithms: Equation (6) resembles the equation of a line: where the parameters m and a are related to Equation (6) and expressed as Equations (8) and (9): where m and a are values from the linear regression of U** and q**. However, there are two equations with three unknowns in the system to be solved. An additional equation, Equation (10), was obtained from Ferguson [4], who used the generalized power law and the bed shear stress.
where U * is the shear velocity (m/s) and D XX is a characteristic grain size (m) taken here as D 84 ; c and b are constant parameters.
In Equation (10), U * and d are replaced with their definitions. After a mathematical process, Equation (11) is obtained, which is the same obtained by Ferguson [4]. This equation contains parameter mo, defined in Equation (12).

Variance Decomposition Methodology (VDM)
VDM decomposes the variance of the total error, as shown in Equation (17). The variance of the total error increases as the model output increases. To obtain a constant variance independent of the model output magnitude, a Box-Cox transformation was applied.
Both observed and predicted values of velocity are transformed using Equation (16). λ is calibrated through the minimization of the variance of the error of the transformed predicted and measured velocities. The resulting lowest variance is taken as the total residual error variance in Equation (17). The observation error variance for velocity is based on an error of 5% obtained by Lee and Ferguson [18].
where Y is the model output variable, S is the standard deviation, Se 2 Y-Yo is the total residual error variance, Se 2 Y is the model error variance, and Se 2 Yo is the observation error variance.

Test with Data from Literature
An additional performance test was developed using data available from the literature. Jarret [15] provides data from 21 reaches in the Rocky Mountains of Colorado. However, data from two reaches could not be used because D 84 data was not available. Bathurst [16] presents data from 16 British rivers. Both data sets were joined comprising 121 measurements with flows ranging from 0.137 to 129 m 3 /s, for this data set the morphology of each reach was not specified.
The proposed methodology used 50% of the data randomly chosen to estimate a 1 , a 2 , and a 3 ( Figure 2) and predict U** with the remaining 50% of the data. Moreover, Zi12010 [14], the best fitting equation for Cascade and Step-pool, was used to predict the same data. In this equation, instead of using S 0 , S F was used since this parameter was provided in the dataset. (3) for 50% of the data provided in Jarret [15] and Bathurs [16] randomly chosen.

Best Fitting Equation
The performance of the empirical equations listed in Table 3 was compared with th measured velocities in the Cascades, Step-pools, and Plane-bed. Table 4 presents the NDHG equations with the best fitting properties for the Cascade, Step-pool, and Plane-bed rive reaches. For Cascades and Step-pools, the Zi12010 equation [14] is the best; for the Plane-bed the Co2009 equation [13] performs best, except for the SX metric for which th FeNHGE2007 [4] fits best. The metrics indicate that performance depends on the morphology. For Cascades, th Zi12010 prediction performs well according to Ritter and Muñoz-Carpena [31]. The rela tive RMSE and MAE are similar, 16% and 12%, respectively. For Step-pools, Zi12010 per forms acceptably. The difference between relative RMSE (23%) and MAE (17%) is th same as for Cascades, approximately 5%, indicating that there are no significant difference between residual magnitudes for these morphologies. For the Plane-bed, Co2009 demon strated unsatisfactory performance, with a higher difference between relative RMSE (30%  (3) for 50% of the data provided in Jarret [15] and Bathurst [16] randomly chosen.

Best Fitting Equation
The performance of the empirical equations listed in Table 3 was compared with the measured velocities in the Cascades, Step-pools, and Plane-bed. Table 4 presents the NDHG equations with the best fitting properties for the Cascade, Step-pool, and Plane-bed river reaches. For Cascades and Step-pools, the Zi12010 equation [14] is the best; for the Plane-bed, the Co2009 equation [13] performs best, except for the S X metric for which the FeNHGE2007 [4] fits best. The metrics indicate that performance depends on the morphology. For Cascades, the Zi12010 prediction performs well according to Ritter and Muñoz-Carpena [31]. The relative RMSE and MAE are similar, 16% and 12%, respectively. For Step-pools, Zi12010 performs acceptably. The difference between relative RMSE (23%) and MAE (17%) is the same as for Cascades, approximately 5%, indicating that there are no significant differences between residual magnitudes for these morphologies. For the Plane-bed, Co2009 demonstrated unsatisfactory performance, with a higher difference between relative RMSE (30%) and MAE (21%) than for Cascades and Step-pools, indicating higher residual values than for the other morphologies. The best equations for all morphologies do not have cases with predicted and observed values (PE) that differ by a factor greater than two or less than 0.5. The model error relative to the observed value (S X ) illustrates that for a Plane-bed the best equation is FeNHGE2007. The Co2009 equation omits S F , and FeNHGE2007 includes S F . The morphology with the lowest model error relative to the observed value is Cascade (17%); Step-pool and Plane-bed have similar values (approximately 25%).

Estimation of NDHG Parameters
Implementation of the proposed methodology to calculate the NDHG exponents first requires a check that log(q**)-log(U**) follows a linear trend. Figure 3 indicates a linear tendency for all of the morphologies, with a coefficient of determination (R 2 ) greater than 0.85 for all fittings, although the slope m and the independent factor a vary considerably between the studied morphologies.
The estimated coefficients of NDHG for each morphology are presented in Table 5. Parameter a 1 varies from 1.73-2.31; the value range of parameter a 2 depends on the morphology and is 0.75 for the Plane-bed, considerably higher for Cascades and Step-pools (0.48-0.57). The opposite is true for a 3 ; the Plane-bed value (0.12) is less than the values for Cascades and Step-pools (0.21-0.26).  Table 6 shows the performance metrics of the NDHG equations according to the proposed methodology (NDHGCA, NDHGSP, and NDHGPB). Comparison of Tables 4 and 6 reveals that the proposed approach produces improved quality metrics. The derived equations demonstrate good performance according to the EF metrics. The difference between relative RMSE and MAE has been reduced to 4.6% on average for all morphologies; the residual magnitudes are uniform. The model error relative to the observed values decreases significantly for Step-pool and Plane-bed morphologies; the improvement in the Cascade morphology is less evident. PE illustrates that the proposed NDHG equations did not produce any point with predicted and observed values that differ by a factor of less than 0.5 or greater than two.   Table 6 shows the performance metrics of the NDHG equations according to the proposed methodology (NDHGCA, NDHGSP, and NDHGPB). Comparison of Table 4 and

Variance Decomposition Methodology (VDM)
The calibration required for the Box-Cox transformation provides the following data: Cascade, λ = 0; Step-pool, λ = 1; Plane-bed, λ = 0. When λ = 0, the Box-Cox transformation is a logarithmic transformation BC(Y) = log(Y). The use of calibrated parameters allows the decomposition of the variance, as shown in Table 7. Table 7 reveals that most of the error variance is contained in the model output. There are slight differences in the variance of observation errors, but analysis can be conducted based on the ratio of the relative to total residual error variance. Hence, the Cascade morphology exhibited the largest observation error variance, followed by Plane-bed and Step-pool. The calibration data was used to calculate the band presented in Figure 4; for all three morphologies, approximately 70% of the data is inside the band. Table 7. Variance decomposition methodology for the studied morphologies.

Cascade
Step-Pool Plane-Bed   Figure 2 depicts the linear regression of U** and q** of randomly chosen 50% of literature data. The regression provides the information needed to estimate a1, a2, and a3. The obtained equation was called NHDGlit.   Figure 2 depicts the linear regression of U** and q** of randomly chosen 50% of literature data. The regression provides the information needed to estimate a 1 , a 2 , and a 3 . The obtained equation was called NHDGlit. Table 8 compares the performance of  NHDGlit and Zi12010. In this Table, three metrics have been used: RMSEa and MAEa are dimensionless versions of RMSE and MAE, defined as a percentage of the observations mean, and EF is the Nash-Sutcliffe efficiency index. RMSEa and MAEa depict a marked predictive superiority of NHDGlit against Zi12010. According to EF, NHDGlit has a very good performance rating, however, Zi12010 EF shows an unsatisfactory performance rating [31].

Characteristics of NDHG Equations
Bathurst [5] suggested that a resistance equation needs two parameters, one representing at-a-site resistance variation, and the other representing between-site resistance differences. At-a-site variations are usually related to the relative submergence d/D 84 (Maxwell and Papanicolaou, 2001). According to Bathurst [16], D 84 provides a 3-D image of the bed material disposition. However, Aberle and Smart [2] state that D 84 is not a good roughness height for Step-pools. Indeed, q* is considered a better at-a-site parameter, as any measurement error affects the observed and predicted values [4]. In this study, the same function is attributed to q**, as the equation structure is the same. The difference between q* and q** is that q** comprises S F in its denominator. The between-sites parameter, S F [14], represents the change in morphology [5]. Moreover, it is expected that the exponents in non-dimensional hydraulic equations change with morphology [1].
In NDHG equations, the dimensionless unitary flow is less sensitive to measurement errors [4], and its combination with S 0 explains most of the resistance variation [1,8]. These equations do not assume any distribution of velocity or resistance parameter [14]. The velocity distribution in mountain rivers has an S-shape [16]; in low-land rivers, the velocity distribution is semi-logarithmic. When Equation (1) is used to calculate velocity, there is an assumption of uniform flow; this is not the case for mountain rivers, which are characterized by changes in water depth and surface slope at each XS [16]. However, there is no better alternative for relating a resistance parameter with velocity. Dimensionless equations are preferred for the following reasons: (1) the exponents are also dimensionless; (2) in these equations, the common physics for all of the reaches are taken from empirical data [32,33].

NDHG Parameters
Previous studies proposed constant exponents for NDHG equations [4,13,14]. However, Nitsche et al. [11] identified variability in a 1 . Nitsche et al. [11] used data from six reaches with different morphologies and eight other Swiss mountain streams. In addition, a 1 depends on the concentration of boulders Γ; the term containing slope was not included in their equation and a 2 was fixed as 0.6. Nevertheless, we found a 1 , a 2 , and a 3 to have different values for each studied site. According to the proposed methodology, these parameters depend on the regression parameters of log(q**)-log(U**) and the energy slope. Hence, this methodology requires collecting field data at different flow magnitudes, which is not possible in all cases. The authors consider that regression parameters m and a may be related to bed material or profile characteristics. Lacking sufficient data, we could not conduct this analysis. However, it would be possible to find an expression for m and a with additional data.
NDHG parameters depict certain patterns as a function of the reach morphology. The order of magnitude of a 2 (~0.5) and a 3 (~0.2) for Cascades and Step-pools are the same as in the literature; this is not the case for the Plane-bed. The relation between a 2 and a 3 is key to the importance of at-a-site variations compared to between-site variations in resistance. In Cascades and Step-pools, the at-a-site variation parameter has an exponent (a 2 ) that is 2.4 times larger on average than the between-site variation parameter (a 3 ). In a Plane-bed, this difference increases to 6.02 (a 2 /a 3 ), illustrating that the effect of between-site variations in resistance is not relevant in a Plane-bed, which is logical, as there are no periodic bedforms in a Plane-bed. Ferguson [4] found an a 2 /a 3 ratio of 3:1 in pool-riffles, riffles, runs, and Step-pool reaches. This ratio is close to the value obtained for Cascade 2 in Table 5. Zimmermann [14] conducted flume experiments in self-formed Cascade streams with a resulting a 2 /a 3 ratio of 1.71, similar to the result for Cascade 1 in Table 5. Thus, the relations estimated with the proposed methodology are similar to those found in the literature for field and flume experiments.
The values for a 1 obtained with the proposed method are clearly higher than those presented in the literature. The a 1 value in our research varies from 1.83-2.87. Ferguson [4] provided a value in the range of 1-1.74. Zimmermann [14] obtained a value of 1.45. The difference in the values may be due to the data used to derive the equation. In this study, when a 1 , a 2 , and a 3 were obtained, each site was analyzed independently. However, in the other studies, data from different reaches were used to calibrate the equations. The separation of the reaches is due to the different at-a-site and between-site variations of resistance, which led us to believe that the parameters are dependent on the regression parameters and relate to the reach characteristics.
The values of m and a follow an evident pattern. Plane-bed values are higher than Cascade and Step-pool values. According to the results, this may be explained by smaller m and a values at higher resistance complexity. Table 7 represents the variance decomposition of the NDHG equations for the analyzed morphologies. The trend shows that the variance of the observation error is higher in Cascades and Plane-beds than in Step-pools. Cascade and Plane-bed conductivity sensors cannot be installed in places without turbulence; these morphologies have bed material of significant size, which contributes to turbulence. In contrast, Step-pool data is collected from pools with a nearly stationary flow and smaller bed material. Moreover, Cascades and Plane-beds have the same calibration parameter for Box-Cox transformation (λ = 0);

Variance Decomposition Methodology (VDM)
Step-pools have a calibration parameter of λ = 1. Table 7 clearly indicates that model error variance (Se 2 Y ) is the main component in the variance decomposition methodology. This term comprises different model output variances resulting from the model structure, the input, and the parameters. The model structure component is expected to be small because the NDHG equations represent the best equations for all examined morphologies. There was not a better equation structure for the performance of these equations. The inputs for this model are flow, gravity, wetted width, energy slope, and D 84 . The energy slope is approximated with the water level; the water level had an error of 1.5% of the standard deviation uncertainty. Wetted width had less than 0.14% of the standard deviation uncertainty. These values were computed through repeated measurements from different morphologies. Flow has an error of 5% according to Lee and Ferguson [18], who used tracers for flow and velocity measurement. The bed material 84th quartile (D 84 ) was obtained after sampling 400 particles at each reach to obtain the bed material distribution. Some studies sampled only 100 particles in studying poolriffles and boulder-cobble beds [5,16,34]. The number of samples increased to 300 particles when additional morphologies were examined [8]. It is evident that increasing the number of samples decreases the estimation error in the bed material distribution, given that the input parameters have a small influence on the model error variance. Hence, most of the model error variance corresponds to the parameters a 1 , a 2 , and a 3 . Table 4 presents the performance of different versions of NDHG equations with diverse parameter values. Figure 4 shows the band encompassing the mean error +/-standard deviation, which for a normal distribution includes 68.26% of the data. Given that 70% of the points fall inside the band for all morphologies, the normal distribution assumption is justified and the confidence interval of the model equals 70%.
The use of standard deviation of bed elevations (s) has been tested as characteristic roughness length in multiple studies, however, s provides good results only in some of them. Aberle and Smart [2] have successfully used s (s from 4.6 to14.6 mm) as a roughness parameter in a power equation in Step-pools to predict (8/f) 0.5 based on flume data. Lee and Ferguson [18] (s ranged from 0.068-0.257) use field and flume data from Step-pools to find an equation to predict (1/f ) 0.5 . A log law equation provided good results when the effective roughness was step D 84 . Nitsche [11] (s lies in the range of 0.07 to 0.47) found an NDHG equation whose dimensional macro roughness parameters to calculate velocity and dimensionless unitary flow was D 84 . In this study (s range 0.022 to 0.214) different equations were tested, some of them with different representations of effective roughness, but NDHG equations provided the best fitting, and a new methodology to determine its parameters has improved fitting performance. All the NDHG equations used in this research utilize D 84 to estimate velocity and dimensionless unitary flow, so this term seems to work well for these types of equations.

Conclusions
In this research, the flow resistance in three morphologies of a headwater mountain river was studied, in three Cascades, two Step-pools, and one Plane-bed. Each reach was divided into three to five XS where staff gauges were installed. Field measurements of water level and wetted width were collected at each XS. For each reach, flow and mean velocity were computed for different water level conditions. Different empirical equations for velocity prediction were tested using goodness-of-fit metrics. The equation with the best performance was calibrated to find expressions for its exponents. Moreover, variance decomposition methodology was used to estimate the uncertainty of the proposed method. As a final step, data from the literature was used to test the proposed methodology.
The findings clearly indicate that the best equations for the studied morphologies are NDHG equations; no other type of equation exhibited similar performance. A methodology to find the NDHG exponents was proposed using logarithmic regression, the bed shear stress, and the generalized power law. The resulting equations for the exponents in the NDHG equations have variations according to the type of reach and depend on regression parameters, namely the slope and the independent terms (m, a). The derived Step-pool and Cascade ratios (a 2 /a 3 ) are in accordance with data from the literature; for the Planebed, this ratio is larger due to the small influence of the between-site resistance variation component in this morphology. Besides, the proposed methodology was successfully used to predict data from the literature. The applicability of the proposed approach for estimating the exponents of NDHG equations can certainly be improved with additional data (experimental measures and other morphologies). As the m and a parameters may follow a certain pattern, this methodology is useful for ungauged streams. Funding: This research was developed within the framework of the project "Modelación hidrodinámica avanzada como herramienta para predicción de calidad de agua en ríos de alta pendiente" funded by the Research Directorate of the University of Cuenca (DIUC) under the XVth call for research proposals.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to a confidentiality agreement of the financing project "Modelación hidrodinámica avanzada como herramienta para predicción de calidad de agua en ríos de alta pendiente".