Investigating the Sediment Yield Predictability in Some Italian Rivers by Means of Hydro-Geomorphometric Variables

In the present work, preliminary results are reported from an ongoing research study aimed at developing an improved prediction model to estimate the sediment yield in Italian ungauged river basins. The statistical correlations between a set of hydro-geomorphometric parameters and suspended sediment yield (SSY) data from 30 Italian rivers were investigated. The main question is whether such variables are helpful to explain the behavior of fluvial systems in the sediment delivery process. To this aim, a broad set of variables, simply derived from digital cartographic sources and available data records, was utilized in order to take into account all the possible features and processes having some influence on sediment production and conveyance. A stepwise regression analysis pointed out that, among all possibilities, the catchment elevation range (Hr), the density of stream hierarchical anomaly (Da), and the stream channel slope ratio (∆Ss) are significantly linked to the SSY. The derived linear regression model equation was proven to be satisfactory (r2-adjusted = 0.72; F-significance = 5.7 × 10−8; ME = 0.61), however, the percentage standard error (40%) implies that the model is still affected by some uncertainties. These can be justified, on one hand, by the wide variance and, on the other hand, by the quality of the observed SSY data. Reducing these uncertainties will be the effort in the follow-up of the research.


Introduction
The assessment of river sediment yield (SY) is a focal issue in geomorphological research.This is because of the off-site effects produced by sediment routing in the fluvial system, e.g., reservoir siltation, degradation of aquatic habitats, overflooding, coastal erosion/accretion etc.In the matter of land planning and natural hazard management, being able to quantify the magnitude of these processes is of primary importance, given their repercussions on human activities.Unfortunately, monitoring data records are often lacking because of the absence of an adequate gauging network.
Consequently, theoretical models are needed to fill this gap.For years, this topic has been addressed by many authors who proposed different solutions for sediment yield prediction, from simple parametric equations to refined physically-based models.These latter are commonly utilized by hydraulic engineers mainly for design purposes.Eurosem [1], Système Hydrologique Européen (SHE) [2], Medrush [3], Watem/Sedem [4] are some physically-based models that have been developed in Europe in the last decades (see also [5,6] for the state-of-the-art about prediction models).An important limit of such models is that they are heavily data-requiring and time-consuming, so then they are generally suitable for fields and small catchments.
The quantitative geomorphology introduced by Horton [7], Strahler [8][9][10], Schumm [11], and others has allowed for parameterizing the dimensional, morphological, and topological characteristics of river networks and catchments, providing a quantitative description of the physiography of drainage basins.These parameters can theoretically represent meaningful predictors in river dynamics modeling, as largely proven in the worldwide literature since the end of the 1950s [12][13][14][15][16][17][18][19][20][21][22][23][24].Therefore, an alternative parametric-geomorphologic approach has been introduced in the sediment yield prediction issue, via the application of statistical correlation analysis and the development of simple or multiple linear regression models.Differently from the physical approach based on fluid and sediment dynamics, such methodology looks notably profitable and parsimonious with regards to its applicability at any scale.Moreover, the development of modern image processing technologies and geographical information systems (GIS) has simplified the computation of geomorphic parameters, along with the geographic-spatial handling of elevation, soil, geology, land cover, and climatic-hydrologic data, and a new discipline, geomorphometry, was born.In this field, new promising techniques, such as the DEMs of Difference (DoD) or Structure for Motion (SfM), based on high-resolution digital elevation models (DEMs) from LiDAR, have been investigated in the latest research.Interesting achievements have been obtained in the view of improved prediction models of sediment erosion and transport at hillslope and subcatchment scale [25][26][27].
With regard to the statistical approach, several regression model equations have been proposed in Italy, during the past years [28][29][30][31][32][33][34][35][36].Nevertheless, no one of such models can thoroughly account for SY variability in Italian watercourses, to date, since they are affected by many uncertainties deriving from how most of these models have been developed and presented.One reason is represented by the limited number of observations considered so far, mostly referred to restricted regional areas, as in [28,31,[34][35][36].Moreover, no significance or model robustness tests were reported in most of the mentioned works.Gazzolo & Bassi [29] carried out an extended investigation on 53 drainage basins located over the entire Italian territory from the Alps to Sicily, but the analytical expression they found, based in large part on observations prior to 1950s, was not supported by a suitable model performance evaluation.Therefore, extending the investigation to a larger number of catchments and supporting the analysis by suitable statistical diagnostic tests could allow the achievement of better results.
In the present work, the first phase of a study program concerning 45 Italian drainage basins of varied extensions, distributed along the Apennine ridge, is presented.The predictability of sediment yield by means of hydro-geomorphometric variables is investigated by examining a first tranche of 30 catchments located in southern Italy and Sicily.New-generated parameters were included in the analysis, with the aim to possibly improve the previously developed prediction models.Ultimately, an extended dataset comprising a wide spectrum of variables was gathered, with the aim to consider the influence of main catchment physiographic components that, in different ways and to varying degrees, can affect the sediment yield.A statistical analysis supported by confidence and performance tests was then performed to assess the reliability of the considered independent variables in explaining the different behavior of the selected catchments with regard to sediment delivery.

Input Data
To the aim of the present study, a preliminary selection of catchments has been made on the basis of the availability of sediment yield data records covering an adequately long timespan.In Figure 1 and Table 1, the locations and characteristics of the selected catchments are shown.The sediment load measurements were carried out, approximately during the period 1950-1990, by different territorial offices of the former Italian National Hydrographic and Tide Monitoring Service (SIMN), whose data are provided in annual reports publicly available online.After the year 1990, the measurements began to be run discontinuously, hampering gathering of a meaningful dataset.In 2002, the jurisdiction on hydrologic monitoring was officially transferred to the regional authorities, and several gauging stations ceased to operate.The amount of sediment transported to the outlet of selected catchments is, here, expressed through the area-specific sediment yield (SSY), as already done by previous authors, in order to compare catchments with different area extents.SSY data are referring to the suspension load, thus, neither the bedload nor the solute fraction are considered here.This is because the suspended load is more easily measured and, consequently, long monitoring data records can be found, while measurements of bed and solution loads are definitely more rare, and only covering short time periods.Nevertheless, there is a general agreement that the suspended sediment fraction can account for up to 90% or more (depending on bedrock type) of the total river load [18,37].

Input Data
To the aim of the present study, a preliminary selection of catchments has been made on the basis of the availability of sediment yield data records covering an adequately long timespan.In Figure 1 and Table 1, the locations and characteristics of the selected catchments are shown.The sediment load measurements were carried out, approximately during the period 1950-1990, by different territorial offices of the former Italian National Hydrographic and Tide Monitoring Service (SIMN), whose data are provided in annual reports publicly available online.After the year 1990, the measurements began to be run discontinuously, hampering gathering of a meaningful dataset.In 2002, the jurisdiction on hydrologic monitoring was officially transferred to the regional authorities, and several gauging stations ceased to operate.The amount of sediment transported to the outlet of selected catchments is, here, expressed through the area-specific sediment yield (SSY), as already done by previous authors, in order to compare catchments with different area extents.SSY data are referring to the suspension load, thus, neither the bedload nor the solute fraction are considered here.This is because the suspended load is more easily measured and, consequently, long monitoring data records can be found, while measurements of bed and solution loads are definitely more rare, and only covering short time periods.Nevertheless, there is a general agreement that the suspended sediment fraction can account for up to 90% or more (depending on bedrock type) of the total river load [18,37].In addition to the SSY data, rainfall and overland runoff data records were also gathered, paying attention to the chronological correspondence between the sediment yield and the hydrologic data series.As with the SSY data, rainfall and overland runoff data were collected from the SIMN annual reports as well.
Other input data were the fluvial network lines and catchment divides.These were derived, respectively, from the 1:25,000 scale official topographic maps edited by the Italian Army's Geographic Institute (IGMI) and the 20 m × 20 m resolution digital elevation model (DEM) of Italy.Further basic information data regarded soil loss and geolithological features of the selected catchments.Geolithologic descriptions were derived from the Geolithological Map of Italy, publicly available on the Italian National Geoportal [38].Soil loss data were extracted from the 100 m × 100 m resolution Map of Soil Loss by Water Erosion in the European Union [39,40].

Data Processing
For each of the selected catchments, different groups of parameters have been computed (see Tables A1-A7 in Appendix A).One group comprises the geomorphometric parameters, describing the catchment dimensions and morphology along with the stream network texture and topological organization.Another set of variables helps to describe the sediment potential and the catchment ability in sediment delivery.Finally, some hydrologic data and derived indices have been taken into account, in order to consider the driving forces acting in the processes of sediment particles detachment and transfer through the catchment.These forces are basically produced by the energy input given by falling and running water.

Geomorphometric Analysis
Geomorphometric parameters were derived from cartographic sources or DEM by measuring basic catchment features, such as perimeter, surface area, length, and relief, as far as the total length of stream channels forming the fluvial network and ordered in a hierarchic tree-shaped structure.
The calculation of parameters has been performed in QGIS environment by means of the QMorphoStream toolset [41], specifically designed and developed with the aim to integrate, in a single logical frame, all the steps of the geomorphic quantitative analysis.The toolset integrates original algorithms with both QGIS built-in functions and processing tools available in the open-source GIS community.The parameters computed by the QMorphoStream tools are listed and described in Tables A1-A4 (Appendix A).For a more detailed description and explanation, see the cited references.
In addition to the well-known parameters described in the literature [7][8][9][10][11][42][43][44], a series of new geomorphometric variables have been introduced, with particular reference to the stream channel slope (see Table A1 in Appendix A).In this regard, new parameters were derived in various ways, based on the average slopes of streams grouped by Strahler's orders and computed by means of the specific QMorphoStream processing tool through the intersection of the stream network layer with the 20 × 20 m resolution DEM.The stream slope data can be averaged for the whole concerned catchment as simple mean (Ss ave ) or as weighted mean, on the basis of the Strahler's order (Ss O ), or of the frequency of streams in a given order (Ss N ).Moreover, considering the roles played, respectively, by the low-order streams of the upper (mountain) portion of the catchment, characterized by steeper bedslopes more favorable to sediment transport, and by the high-order streams in the lower valley treats, whose lower slopes foster the sediment deposition, two different averaging parameters, respectively, Ss up and Ss down , have been developed.The belonging to one of the two categories (up or down valley) is established on if the average slope of the i-th order streams is greater or lower than the average stream slope (see Table A1 of the Appendix A for explanation).This distinction was made in order to differentiate the catchments with large alluvial plains from the prevailingly mountainous catchments with limited alluvial plain, assuming that the first would undergo a relative reduction in specific sediment delivery compared to the latter.To better express this circumstance, a single index ∆S s , referred to hereafter as "stream channel slope ratio", was also introduced, given by the Ss down /Ss up ratio.The rationale for this new parameter is that, with average upvalley stream gradient being equal, if average downvalley stream gradient decreases, ∆S s decreases too, then, catchments with large alluvial plain will most likely show lower ∆S s .On the other hand, for different watersheds having quite similar values of average downvalley stream slopes, a decrease in ∆S s indicates, on average, higher values of upvalley stream slopes.
After all, the physical meaning of the stream slope parameters is analogous to the H y index, that is, the integral of the hypsometric curve representing the distribution of elevations within the catchment [9] (Figure 2).The shape of the hypsometric curve, that can be also parameterized in terms of statistical moments or density functions [45], as well as the value of the hypsometric integral, are indicative of the morphological-evolutionary stage of the catchment, implying the state of equilibrium/disequilibrium between erosion, transport, and sedimentation rates [9,46].
On the other hand, the stream gradient can also account for the bedrock lithology.According to Hack [47], the stream slope S is strictly connected with the stream length L through a couple of coefficients k and n, changing in function of bedrock lithology: This correlation ultimately means that the stream longitudinal profile conveys the link between catchment constitutive materials and landforms, as a result of fluvial processes [18].

Sediment Supply and Connectivity
Recently proposed indices [48] referred to as "sediment supply potential" (SSP) and "simplified connectivity index" (SCI) have been also computed, on a geographical basis, by taking into account the inverse relative distance, measured along the stream network, of each soil loss unit area from the river outlet (see Table A5 of Appendix A).In particular, the SSP is the weighted mean of soil loss rates per unit area while the SCI index is devised to quantify, in a simplified mode, the catchment efficiency in delivering the soil-loss amounts to the outlet.In addition to these indices, the Lem index of bedrock erodibility has also been utilized.This index is based on a score-rating of the proneness to erosion (Es), owned by the different geolithological units recognized in the selected catchments [36].The erodibility scores, here utilized, are reported in Table A6 of Appendix A. Furthermore, considering the functional relationship between bedrock lithology and landscape, as theorized by Hack [47], and Cooke & Doornkamp [18], before cited, a new parameter (LHR) was also introduced, given by the ratio between the Lem index and the hypsometric integral Hy, in order to couple the bedrock erodibility with the catchment morphology.

Hydrologic Parameters
The mean yearly rainfall (Ptot), the total catchment overland flow and the total overland rainfall depth data were finally considered.Both Ptot and the total overland rainfall depth can be seen as input quantities, in the catchment energy balance, while the corresponding output is represented by the total overland flow at the outlet.Moreover, the flow-coefficient (overland flow/rainfall depth ratio) also appeared to be useful to the present analysis, since it can be considered as a measure of the catchment efficiency in transforming precipitations in overland water flow, that is, the medium by which the sediments move downvalley.In addition, the average rainfall erosivity (R), as defined in the Universal Soil Loss Equation (USLE) [49,50], was computed basing on annual, maximum daily, and maximum hourly rainfall data records, according to the correlation model proposed by Diodato [51].Additional indices were also computed by multiplying both Ptot and R by their own standard deviation i, with the aim to take into account their variability during the observed time periods [32,33].All the above introduced parameters are described in Table A7 of Appendix A.

Sediment Supply and Connectivity
Recently proposed indices [48] referred to as "sediment supply potential" (SSP) and "simplified connectivity index" (SCI) have been also computed, on a geographical basis, by taking into account the inverse relative distance, measured along the stream network, of each soil loss unit area from the river outlet (see Table A5 of Appendix A).In particular, the SSP is the weighted mean of soil loss rates per unit area while the SCI index is devised to quantify, in a simplified mode, the catchment efficiency in delivering the soil-loss amounts to the outlet.In addition to these indices, the L em index of bedrock erodibility has also been utilized.This index is based on a score-rating of the proneness to erosion (E s ), owned by the different geolithological units recognized in the selected catchments [36].The erodibility scores, here utilized, are reported in Table A6 of Appendix A. Furthermore, considering the functional relationship between bedrock lithology and landscape, as theorized by Hack [47], and Cooke & Doornkamp [18], before cited, a new parameter (LHR) was also introduced, given by the ratio between the L em index and the hypsometric integral H y , in order to couple the bedrock erodibility with the catchment morphology.

Hydrologic Parameters
The mean yearly rainfall (P tot ), the total catchment overland flow and the total overland rainfall depth data were finally considered.Both P tot and the total overland rainfall depth can be seen as input quantities, in the catchment energy balance, while the corresponding output is represented by the total overland flow at the outlet.Moreover, the flow-coefficient (overland flow/rainfall depth ratio) also appeared to be useful to the present analysis, since it can be considered as a measure of the catchment efficiency in transforming precipitations in overland water flow, that is, the medium by which the sediments move downvalley.In addition, the average rainfall erosivity (R), as defined in the Universal Soil Loss Equation (USLE) [49,50], was computed basing on annual, maximum daily, and maximum hourly rainfall data records, according to the correlation model proposed by Diodato [51].Additional indices were also computed by multiplying both P tot and R by their own standard deviation σ i , with the aim to take into account their variability during the observed time periods [32,33].All the above introduced parameters are described in Table A7 of Appendix A.

Results and Discussion
After computing all the parameters described above, an extended set of 48 variables was ultimately gathered, whose efficiency and significance in predicting the SSY was analyzed through statistical techniques.Basically, a stepwise linear regression between the dependent (the SSY data) and the independent variables (the parameters listed in Tables A1-A7 in the Appendix A) from the whole sample of 30 catchments was carried out for selecting the most significant variables to enter in the prediction model.
The selection of variables was here made on the basis of their F-significance shown at each regression step.The stepwise procedure was stopped at the third predictor, since it was observed that further inputs would not have significantly improved the model performance while reducing, at the same time, the model significance.
From the stepwise regression, the catchment elevation range (H r ), the density of hierarchic anomaly (D a ) and the stream channel slope ratio (∆S s ) appeared, in order, as the most eligible predictors for SSY.
A multiple linear regression model was thus derived, having the form: It is to be noted that the selected predictor variables are different from those commonly pointed out as strongly significant, such as the drainage density or the total rainfall.On the other hand, it is obvious to presume that all the computed parameters and variables comprised in the performed analysis hold some link with and physical influence on the sediment yield process, but it could be arbitrary to establish, a priori, which are the most significant.In this regard, the catchment elevation range (H r ), pointed out as mostly significant from the statistical analysis, has its own reason, from the physical point of view, since the higher it is, the higher is the potential energy driving the erosional and transport processes.On the other hand, an indicator of network disorder, such as the density of hierarchic anomaly (D a ), which is strongly influenced by the number of anomalous low-order streams, may be associated to the probability of overflooding events (e.g., first-order, non-perennial streams directly draining into main-valley river course, fed by heavy shower and downloading huge amounts of sediments hampering the normal flow) which may justify the negative sign in the relationship with SSY.Lastly, the negative relationship of the stream channel slope ratio (∆S s ), for which an increase in SSY corresponds to a reduction in ∆S s , appears more complex to explain.Basing on how the ∆S s ratio was defined (see Table A1 in Appendix A), a reasonable explanation to this behavior outcomes from the trend analysis of stream slopes in the examined catchments.In fact, in the observed sample, the average slope of upvalley, low-order streams (Ss up ) grows more than the average slope of downvalley, high-order streams (Ss down ), as shown in Figure 3.This means that the sediment delivery capacity from hillslope streams grows in parallel, hence, the observed SSY increases accordingly.
Discussing in detail the significance and correlation degree of each of the 48 computed parameters with SSY is far from the aim of the present work.One remark can be made about the very poor correlation shown by the simplified connectivity index (SCI) and the sediment supply potenial (SSP) derived from soil loss data.These latter are referred to the period after 1990, thus, they could reflect land-use conditions likely different from those existing in the observation period here concerned.This could also explain the better correlation shown by the L em index that, being referred to inherent characteristics of bedrock material, may be considered time-independent, and more suitable for long-term assessments.
To verify the model performance and accuracy of the regression Equation (2), the ANOVA test was performed.The results are shown in Table 2.As can be seen, the statistics look very good with respect to regression significance and determination.In this regard, it should be noted that the application of two of the most significant regression equation among those previously developed for Italian rivers, respectively, Equation ( 9) by Ciccacci et al. [32] and Equation ( 2) by the same authors [33], both utilizing the drainage density and the index of hierarchic anomaly, and even used until recently [52], provided a significantly lower accuracy compared to the present model.Namely, the mentioned regression models, recalibrated on the data set here considered, showed adjusted-r 2 values of 0.132 and 0.205, respectively, with model efficiencies ME [53] less than zero, meaning that the observed mean is a better predictor than the models.The efficiency of the regression model ( 2) here developed, on the contrary, is quite satisfactory.Nevertheless, at the same time, the standard error (40% of mean observed SSY) points out some uncertainty.This would suggest that the chosen variables or a single equation cannot fully explain the variance in the SSY data distribution.Actually, in the latter, this is very high, changing from about 10 Mg•km −2 •year −1 to more than 2000 Mg•km −2 •year −1 ; consequently, a single equation can very hardly predict such wide variations from a catchment to another.As already suggested by Anderson [12], one solution to make the statistical procedure more reliable can be given by introducing the principle of similarity between catchments, i.e., discriminating the catchments on the basis of definite characteristics so as to individuate different model groups which will incorporate different variables.To do that, the catchment elevation range (Hr) can reasonably represent a distinctive feature for catchments grouping, given its highest correlation coefficient with SSY (r = 0.69) in the parameter dataset (see Table 3).As can be seen, the statistics look very good with respect to regression significance and determination.In this regard, it should be noted that the application of two of the most significant regression equation among those previously developed for Italian rivers, respectively, Equation ( 9) by Ciccacci et al. [32] and Equation ( 2) by the same authors [33], both utilizing the drainage density and the index of hierarchic anomaly, and even used until recently [52], provided a significantly lower accuracy compared to the present model.Namely, the mentioned regression models, recalibrated on the data set here considered, showed adjusted-r 2 values of 0.132 and 0.205, respectively, with model efficiencies ME [53] less than zero, meaning that the observed mean is a better predictor than the models.The efficiency of the regression model (2) here developed, on the contrary, is quite satisfactory.Nevertheless, at the same time, the standard error (40% of mean observed SSY) points out some uncertainty.This would suggest that the chosen variables or a single equation cannot fully explain the variance in the SSY data distribution.Actually, in the latter, this is very high, changing from about 10 Mg•km −2 •year −1 to more than 2000 Mg•km −2 •year −1 ; consequently, a single equation can very hardly predict such wide variations from a catchment to another.As already suggested by Anderson [12], one solution to make the statistical procedure more reliable can be given by introducing the principle of similarity between catchments, i.e., discriminating the catchments on the basis of definite characteristics so as to individuate different model groups which will incorporate different variables.To do that, the catchment elevation range (H r ) can reasonably represent a distinctive feature for catchments grouping, given its highest correlation coefficient with SSY (r = 0.69) in the parameter dataset (see Table 3).The Figure 4 shows the trend of H r values, plotted in ascending order.As can be seen, from the distribution of H r data, three likely homogeneous groups of catchments can be outlined: one on the left extreme of the diagram, formed by catchments with lowest H r values less than 700 m (from the Alaco to the Venosa catchment); another on the right extreme, formed by catchments with H r greater than 1300 m (from the S. Leonardo V. to the Sinni river catchment); and the third in the central section of the distribution, given by catchments with intermediate H r values (from the Casanova to the Agri G. river catchment).The Figure 4 shows the trend of Hr values, plotted in ascending order.As can be seen, from the distribution of Hr data, three likely homogeneous groups of catchments can be outlined: one on the left extreme of the diagram, formed by catchments with lowest Hr values less than 700 m (from the Alaco to the Venosa catchment); another on the right extreme, formed by catchments with Hr greater than 1300 m (from the S. Leonardo V. to the Sinni river catchment); and the third in the central section of the distribution, given by catchments with intermediate Hr values (from the Casanova to the Agri G. river catchment).Then, the same stepwise procedure was followed for each of the three groups so identified.
It is notable that the resulting model equations contain different predictors, except for the variable Da, present in both Equations ( 3) and ( 5).This means that grouping the catchments in different combinations implies that different parameters or variables assume more significance at the expense of others, since other physical characteristics of the catchments assume more relevance.As shown in Table 4, statistical indices, standard errors, and model efficiency (ME) have been considerably improved, compared to the global model ( 2), but at the expense of the F-significance, probably because of the reduced number of observations per each group.The worst significance is shown by the model Equation ( 5), where only 5 observations are included, too few to consider the obtained model efficiency (ME = 1) reliable and realistically accurate.With these results in mind, the global model (2), in the present phase, seems to be preferable, as regards its overall performance; nevertheless, clustering the catchments in homogeneous groups can be an interesting aspect to investigate in further research, in view of increasing the number of observed catchments.Then, the same stepwise procedure was followed for each of the three groups so identified.The resulting model equations, distinct for lowest, intermediate, and highest H r class group respectively, are the following: It is notable that the resulting model equations contain different predictors, except for the variable D a , present in both Equations ( 3) and ( 5).This means that grouping the catchments in different combinations implies that different parameters or variables assume more significance at the expense of others, since other physical characteristics of the catchments assume more relevance.As shown in Table 4, statistical indices, standard errors, and model efficiency (ME) have been considerably improved, compared to the global model ( 2), but at the expense of the F-significance, probably because of the reduced number of observations per each group.The worst significance is shown by the model Equation (5), where only 5 observations are included, too few to consider the obtained model efficiency (ME = 1) reliable and realistically accurate.With these results in mind, the global model (2), in the present phase, seems to be preferable, as regards its overall performance; nevertheless, clustering the catchments in homogeneous groups can be an interesting aspect to investigate in further research, in view of increasing the number of observed catchments.In Figure 5, the observed SSY values, taken from the input dataset, are plotted against the SSY values calculated applying the global model ( 2) to each catchment.Despite a certain scatter around the theoretical line of perfect agreement, the model allows a good correlation between observed and calculated SSY values; an ideal regression line would show a determination coefficient r 2 equal to 0.71.In Figure 5, the observed SSY values, taken from the input dataset, are plotted against the SSY values calculated applying the global model ( 2) to each catchment.Despite a certain scatter around the theoretical line of perfect agreement, the model allows a good correlation between observed and calculated SSY values; an ideal regression line would show a determination coefficient r 2 equal to 0.71.The model robustness was also tested by means of the jackknife technique [54].Thirty submodels having the same form of the full model (2) were derived in successive steps by leaving out one catchment and repeating the same procedure for each of the 30 catchments.Each of the submodels was then applied to the excluded catchment, and the obtained SSY data were compared with the results from the full model.The comparison is given in the diagram in Figure 6, showing an optimal agreement between the test models and the full model.In this case, an ideal regression line of the two datasets would show a determination coefficient r 2 equal to 0.94.The model robustness was also tested by means of the jackknife technique [54].Thirty submodels having the same form of the full model (2) were derived in successive steps by leaving out one catchment and repeating the same procedure for each of the 30 catchments.Each of the submodels was then applied to the excluded catchment, and the obtained SSY data were compared with the results from the full model.The comparison is given in the diagram in Figure 6, showing an optimal agreement between the test models and the full model.In this case, an ideal regression line of the two datasets would show a determination coefficient r 2 equal to 0.94.As a general comment to the presented results, a reasonable adaptation of the SSY data to the derived regression formulae has been observed, despite the presence of prediction errors deviating the results from a full agreement, that can be one of the main criticisms of the method here adopted.
However, in this regard, the inherent variability of the processes investigated and the measurement techniques for sediment transport rates are to be considered.As already pointed out, a wide inter-catchment variance is evident in the SSY data distribution; moreover, the mean annual SSY data themselves are characterized by a relatively high coefficient of variation, with CV values up to 2 (Table 1).In addition, no information is available with regard to the data accuracy, and possible discrepancies between measuring techniques and procedures followed at different gauging stations cannot be excluded.Then, a certain degree of uncertainty is reasonably to be expected, that ends up affecting the model performance.
More in general, the above-introduced uncertainties represent a well-known critical issue for researchers and practitioners dealing with estimation methods for sediment yield, as already synthesized in Delmas et al. [21].There is a huge amount of literature on this topic, and an extensive discussion is beyond the aim of the present work.It is just worth mentioning the focal issues already pointed out by Walling [55] about the problems of temporal and spatial lumping and of the "blackbox" nature of sediment delivery, including the intermittent, non-stationary character of sediment transport.In this regard, the author himself highlighted the need for more detailed empirical investigations.
To date, a number of field investigations and laboratory experiments are reported in the literature, dealing with the uncertainties inherent in the sediment load estimation and measurement.For example, comparisons have been attempted between results from different prediction formulae, commonly used in research and design practice, and in situ measures [56,57].Results obtained using both suspended sediment load and total sediment load prediction formulae [58] showed that the computed values can significantly deviate from the measured ones, with mean ratios between computed and measured values from 0.23 up to more than 5, and, in the worst cases, greater than 10.As a general comment to the presented results, a reasonable adaptation of the SSY data to the derived regression formulae has been observed, despite the presence of prediction errors deviating the results from a full agreement, that can be one of the main criticisms of the method here adopted.
However, in this regard, the inherent variability of the processes investigated and the measurement techniques for sediment transport rates are to be considered.As already pointed out, a wide inter-catchment variance is evident in the SSY data distribution; moreover, the mean annual SSY data themselves are characterized by a relatively high coefficient of variation, with CV values up to 2 (Table 1).In addition, no information is available with regard to the data accuracy, and possible discrepancies between measuring techniques and procedures followed at different gauging stations cannot be excluded.Then, a certain degree of uncertainty is reasonably to be expected, that ends up affecting the model performance.
More in general, the above-introduced uncertainties represent a well-known critical issue for researchers and practitioners dealing with estimation methods for sediment yield, as already synthesized in Delmas et al. [21].There is a huge amount of literature on this topic, and an extensive discussion is beyond the aim of the present work.It is just worth mentioning the focal issues already pointed out by Walling [55] about the problems of temporal and spatial lumping and of the "blackbox" nature of sediment delivery, including the intermittent, non-stationary character of sediment transport.In this regard, the author himself highlighted the need for more detailed empirical investigations.
To date, a number of field investigations and laboratory experiments are reported in the literature, dealing with the uncertainties inherent in the sediment load estimation and measurement.For example, comparisons have been attempted between results from different prediction formulae, commonly used in research and design practice, and in situ measures [56,57].Results obtained using both suspended sediment load and total sediment load prediction formulae [58] showed that the computed values can significantly deviate from the measured ones, with mean ratios between computed and measured values from 0.23 up to more than 5, and, in the worst cases, greater than 10.
Another comparison between suspended sediment loads, estimated using the transport-curve method, and in situ measures taken at 10 USGS stream gauging stations, showed mean annual errors up to 63% in absolute value, with a maximum annual error of 526% [59,60].
A recent experimental study was performed by Alhasan et al. [61], in which results of laboratory tests on dam breaching are compared with sediment transport rates predicted using state-of-the-art formulae.The average ratios between computed and experimentally measured specific sediment loads are comprised between 4.9 and 24.8, depending on the formula adopted.These results are indicative of how it is difficult to provide accurate estimates of sediment discharges, even in controlled experimental conditions.
Another recent experimental research study by Praskievicz [62] has shown the relative sensitivity of sediment yield to soil and watershed characteristics under different precipitation intensities.In this regard, measurements should be more frequent under flooding conditions, when the most part of sediment yield is expected to occur.
On these remarks, it can be concluded that the approximation in SSY prediction obtained using the regression model Equation ( 2), here presented, can be considered quite acceptable, even though more efforts could be made trying to achieve even better results.

Conclusions
The preliminary investigation on the predictability of SSY in Italian rivers by means of simple geomorphometric and hydrological parameters, here presented, has concerned a wide spectrum of variables descriptive of the physiographic characteristics of drainage basins.The main aim of the present study was to test the effectiveness of the geomorphologic-statistical approach, by considering a larger number of catchments encompassing a wider geographical extent than was made in previous research in Italy.The statistical analysis helped to recognize the most significant parameters to enter in the regression model based on the available sediment yield observations.The prediction model derived from the analysis carried out to date was satisfactory and encouraging for continued research, even though a certain degree of uncertainty still persisted.This can reasonably be ascribed, on one hand, to the wide variability inherent in the sediment transport process, depending on the extremely variable combination, in space and time, of various environmental conditions.This implies that a single model equation can hardly predict the whole range of data with the same precision.Conversely, the research of partial models per catchment groups, based on distinctive characteristics, did not adequately improve the results, at least in this phase of the investigation.On the other hand, a significant source of uncertainty is likely to be ascribed to the quality of observed data.No information is available about the accuracy in measurements and uneven techniques and procedures could have affected the consistency of the data records here utilized.That represents a thorny problem for the study and management of river basin systems, and points out the common issue of maintaining an efficient and standardized long-term hydrological monitoring network, mainly with regard to sediment yield measurements.In any case, the advantage given by its simple applicability makes such methodology worthwhile to be further improved, since geomorphometric parameters can be easily derived from a cartographic source or DEM, that is, a suitable alternative to costly field investigations.
In the follow-up of this research, more drainage basins will be included in the analysis, in order to deal with an even more representative number of observed cases, and to try to further reduce the residual uncertainty in the model performance.∑ imax r=i+2 Na i,r × f i,r , where Na i,r are the numbers of stream segments of ith order anomalously draining into segments of rth order and f i,r = 2 r−2 − 2 i−1 for i ≤ r − 2 , as defined in [33] Table A4.Parameters describing the catchment relief and shape., where • Se i : percent area of specific lithotype i with a given erodibility score

Figure 1 .
Figure 1.Geographic location of the selected river catchments.Red dots indicate SY gauging stations, numbered in alphabetical order (seeTable1).

Figure 1 .
Figure 1.Geographic location of the selected river catchments.Red dots indicate SY gauging stations, numbered in alphabetical order (seeTable1).

Figure 2 .
Figure 2. Example of hypsometric curves derived from QMorphoStream.(a) The S. Leonardo M. catchment shows the typical S-shaped curve and hypsometric integral of a basin at a "mature" stage of morphological evolution; (b) the upward concave curve and the low value of the integral of the Canale S. Maria catchment indicate a typical "senile" stage of morphological evolution.

Figure 2 .
Figure 2. Example of hypsometric curves derived from QMorphoStream.(a) The S. Leonardo M. catchment shows the typical S-shaped curve and hypsometric integral of a basin at a "mature" stage of morphological evolution; (b) the upward concave curve and the low value of the integral of the Canale S. Maria catchment indicate a typical "senile" stage of morphological evolution.

Figure 3 .
Figure 3. Upvalley stream slopes vs downvalley stream slopes of the selected catchments (see Table A1 in Appendix A for definition).

Figure 3 .
Figure 3. Upvalley stream slopes vs downvalley stream slopes of the selected catchments (see TableA1in Appendix A for definition).

Figure 4 .
Figure 4. Elevation range (Hr) data distribution in the selected catchments.The red dashed lines mark the hypothetical partition of catchments in three subgroups.

Figure 4 .
Figure 4. Elevation range (H r ) data distribution in the selected catchments.The red dashed lines mark the hypothetical partition of catchments in three subgroups.

Figure 6 .
Figure 6.Results of the jackknife test models plotted vs full model Equation (2).

•Figure 6 .
Figure 6.Results of the jackknife test models plotted vs full model Equation (2).

2 2 1
anomaly: I a = Na N1 , where N 1 is the number of 1st order streams D a km −Density of hierarchical anomaly: D a = Na A Table A3.Parameters describing the network texture.the drainage network, i.e., the sum of all the stream segments length F s km −Stream frequency: F s = N A , where N = total number of stream segments D d km −Drainage density: D d = Ltot A MoD d nd "Modified" drainage density: D d = Ltot √ A

•
orographic coefficient: H f = H mean • tana , where tana = Hmean A being the area above the contour level z and A the total basin area • y = z − Hmin Hrange R c nd Circularity ratio, i.e., the ratio of the basin area to the area of the circle whose circumference is equal to the basin perimeter R h nd Relief ratio: R h = Hr Le R e nd Elongation ratio: R e = Dc Le , where D c is the diameter of the circle equivalent to the basin area Table A5.Parameters describing the sediment supply potential.SL i : specific soil loss from the ith unit area (1 ha); • SL max : maximum specific soil loss in the catchment; • D i : distance of the ith unit area from the outlet; • d max is the maximum distance from the outlet; SSP Mg•ha −1 •year −1 Sediment Supply Potential: SSP = ∑ N

R 1 R = 1
(see[51]) MJ•mm•ha −1 •h −1 •year −Ny ∑ N i=1 EI 30−annual , where EI 30−annual = 12.142 (a × b × c) 0.644 and:•EI 30 is the rain erosivity, according to[49] • N is the number of years • a, b, and c are, respectively, the total annual rainfall, the annual maximum daily rainfall, and the annual maximum hourly rainfall, all expressed in centimeters P tot •σ mm Average yearly rainfall amount multiplied by its standard deviationR•σ MJ•mm•ha −1 •h −1 •year −1Average yearly rainfall erosivity multiplied by its standard deviationOverland flow mmDepth of water-layer having the same volume of water flowed at the catchment outlet in a given time interval and uniformly distributed over the catchment areaOverland rainfall depth mmDepth of water-layer whose volume is equal to the catchment rainfall amount in a given time interval, uniformly distributed over the catchment areaFlow-coefficient nd Ratio of the Overland flow to the Meteoric water inflow (commonly < 1)

Table 1 .
List of catchments and specific sediment yield (SSY) data with reference periods.CV is the coefficient of variation or relative standard deviation of mean annual SSY.

Table 3 .
Correlation matrix between SSY and most significantly correlated geomorphometric parameters (r > 0.1).For symbols, see Appendix A.

Table A6 .
[38]ibility scores assigned to main units reported in the Geolithological Map of Italy[38].

Table A7 .
Rainfall and overland flow parameters.