Generation of Natural Runoff Monthly Series at Ungauged Sites Using a Regional Regressive Model

Many hydrologic applications require reliable estimates of runoff in river basins to face the widespread lack of data, both in time and in space. A regional method for the reconstruction of monthly runoff series is here developed and applied to Sicily (Italy). A simple modeling structure is adopted, consisting of a regression-based rainfall–runoff model with four model parameters, calibrated through a two-step procedure. Monthly runoff estimates are based on precipitation, temperature, and exploiting the autocorrelation with runoff at the previous month. Model parameters are assessed by specific regional equations as a function of easily measurable physical and climate basin descriptors. The first calibration step is aimed at the identification of a set of parameters optimizing model performances at the level of single basin. Such “optimal” sets are used at the second step, part of a regional regression analysis, to establish the regional equations for model parameters assessment as a function of basin attributes. All the gauged watersheds across the region have been analyzed, selecting 53 basins for model calibration and using the other six basins exclusively for validation. Performances, quantitatively evaluated by different statistical indexes, demonstrate relevant model ability in reproducing the observed hydrological time-series at both the monthly and coarser time resolutions. The methodology, which is easily transferable to other arid and semi-arid areas, provides a reliable tool for filling/reconstructing runoff time series at any gauged or ungauged basin of a region.


Introduction
Rivers offer several fundamental goods and services to human society and ecosystems.The knowledge of streamflow generated from river basins at different time scales is required for many hydrological applications, including preliminary studies concerning water supply, hydropower, recreation activities, irrigation, and watershed management.Time series recorded at different sites of a region contain crucial information about the spatial, seasonal, and interannual variability of streamflow, which are all basic elements of water resource systems planning and management.
Recent years have been witnessing an increasing qualitative-quantitative degradation of river water due to the coupled effect of human activities (e.g., land use changes or overexploitation of water resources) and climate change.New disciplines, such as socio-hydrology (e.g., [1,2]) or others related to issues of the sustainability and conservation of ecosystems (e.g., [3][4][5][6]), are being developed with the aim of studying the complex mutual interactions between water and society.New efficacious tools are Water 2016, 8, 209 3 of 24 provided by Cutore et al. [47].In their study two regionalization procedures were compared through application to a Sicilian river basin.
In this paper, a regional regressive model is developed and applied to the entire island of Sicily (Italy).The model is based on a simple, non-linear, regression equation whose parameters can be derived by opportune regional equations as a function of some climatic and physical attributes of the basins.The study involves almost all the gauged basins of Sicily, and encompasses a final validation of six representative basins, not considered during the model calibration.
A first calibration phase involves the definition of a simple model structure, valid across the different basins of the region, and the assessment of "optimal" model parameter sets for all the calibration basins, based on the best agreement between observed and simulated runoff series.The successive phase is aimed to define and calibrate regional regression equations relating such optimal parameters to some significant basin attributes, identified among a set of previously selected possible predictors.
Natural runoff time series are often requested in water resource assessment, water quality, and river ecology studies.This paper proposes a pragmatic solution for the estimation of monthly runoff at ungauged sites, analyzing its applicability and limitations.Model accuracy and robustness are verified through different performance indexes and specific analyses at monthly and coarser time resolutions (i.e., seasonal and annual).The methodology can be easily transferred to other arid and semiarid Mediterranean areas.At the same time, it provides a prompt and reliable tool for the evaluation of the available water resources at any basin of Sicily, a region generally characterized by water scarcity and where similar tools are particularly needed and urgent.The model could also be used in river ecology studies, for example with regard to regulated river basins and with the aim of estimating to what extent the hydrological regime has been modified compared to the natural one.

Study Area
The selected study area is the region of Sicily (Italy), which is the largest island in the Mediterranean Sea, with an area of over 25,700 km 2 .Despite the presence of several singularities, the region can be summarily distinguished into three different topographic zones: the northern coast area, mostly mountainous (principal mountain ranges: the Madonie, the Nebrodi, and the Peloritani); the southern coast area, essentially hilly; and the eastern coast area, dominated by the Hyblaean Mountains and the largest active volcano in Europe, i.e., Mount Etna.
The island is characterized by a typical Mediterranean climate, which, according to the Köppen-Geiger climate classification, belongs to the warm temperate class (group: Csa-Hot-summer Mediterranean climate), with, usually, a hot and dry summer and a maximum of precipitation during the colder wintry months.Nevertheless, considerable differences in the thermo-pluviometric regime can be observed among different sub-areas within the region, especially if rainfall and temperature are analyzed at monthly or lower time scales (i.e., daily); in general, the northern coast area is rainier and colder than the southern and eastern coastal areas.
Sicily is drained by several rivers, especially in the central area.The climate and the morphology in the region concur in determining a prevalence of small river basins, usually characterized by non-perennial hydrological regimes.Rivers in the northern coast zone are numerous, short, and ephemeral, while some rivers in the southern cost zone drain wider basins and are perennial.The eastern zone is characterized by some rivers that are relevant in terms of annual discharge.A considerable portion of the region preserves natural conditions in terms of land use: northern Sicily is characterized by the presence of some wide natural areas (e.g., Ficuzza Nature Reserve or Nebrodi Regional Park), where the main forests of the region are settled, while the non-cultivated natural areas in southeastern Sicily are mainly populated by vegetation typical of arid and semi-arid Mediterranean regions.

Data Sources and Regional Database
The OA-ARRA (Osservatorio delle Acque-Agenzia Regionale per i Rifiuti e le Acque) regional archive provides long historical data series for almost 100 Sicilian river basins, collected by a dense network of gauging stations across the territory (almost 90 hydrometric and 250 meteo-climatic gauging stations).Initially only the basins with at least five continuous years of simultaneous runoff-rainfall-temperature monthly observations have been considered for this study; this reduces the number of examined watersheds to a total of 71 gauged basins (i.e., 39 main basins and 32 sub-basins), as reported in Figure 1a.The assigned identification number (ID), the names of the watershed and the relative hydrometric gauging station, the basin typology (main or sub-basin), and the observation period for each basin are reported in Table 1.The location of all the OA-ARRA meteo-climatic stations are depicted in Figure 1b, along with reporting information about the station typology (Temperature, Rain-Temperature, or Rainfall gauging stations). of gauging stations across the territory (almost 90 hydrometric and 250 meteo-climatic gauging stations).
Initially only the basins with at least five continuous years of simultaneous runoff-rainfall-temperature monthly observations have been considered for this study; this reduces the number of examined watersheds to a total of 71 gauged basins (i.e., 39 main basins and 32 sub-basins), as reported in Figure 1a.The assigned identification number (ID), the names of the watershed and the relative hydrometric gauging station, the basin typology (main or sub-basin), and the observation period for each basin are reported in Table 1.The location of all the OA-ARRA meteo-climatic stations are depicted in Figure 1b   For each basin, historical monthly series of runoff depth at the basin outlet and simultaneous monthly areal precipitation have been retrieved from the official OA-ARRA annual reports (i.e., Annali Idrologici).The corresponding series of mean areal temperature have been estimated starting from the data collected by all the OA-ARRA thermometric stations (Figure 1b) and using the same spatial interpolation techniques used in Pumo et al. [48]; more specifically, algorithms based on a station weighting technique [49,50], which also take into account the elevation, have been used.
The proposed methodology requires a preliminary identification of significant basin descriptors, which could be suitable for the multi-regressive analysis of the second calibration phase.An appropriate basin descriptor selection set may be important in regionalization method performance depending on the adopted regionalization approaches [51].A preliminary analysis has been carried out for a subset of 20 basins, analyzing several characteristics related to the size, morphology, soils, vegetation, land use, and climate of the basins, and chosen mainly on the basis of their availability and ease of retrieval.In particular, a set of over 30 different basin attributes has been derived from the OA-ARRA regional database or by properly clipping some thematic maps available for the area of interest, such as: a 100-m resolution DEM (Digital Elevation Model) of Sicily (Figure 2a); the Corinne Land Cover Map of Sicily of 2006 [52]; the Soil Map of Italy and Soil Data Base of Italy's Soil Information System (SISI, [53]); the regional Land Use Map of Sicily (at 1:250,000 scale); the Regional Map of CN (Figure 2b), created in 2004 during the editing of the regional hydrological system plan PAI (Piano di Assetto Idrogeologico) and reporting the Curve Number, i.e., the parameter adopted in the SCS (Soil Conservation Service) Runoff Curve Number method [54].Such basin attributes, first examined for cross-correlations, were then individually analyzed in terms of their ability to explain the spatial variability of the basin's hydrologic response; more specifically, the Pearson correlation coefficient has been computed for all the possible combinations of basin attributes (independent variable) and three statistical values (mean, standard deviation, asymmetry) relative to the monthly observed runoff (dependent variable).The least significant explanatory variables have been dropped, selecting only basin attributes with sufficiently high coefficients (over 0.3, significant at 0.05 level) for at least one of the three examined statistics.
Seven basin descriptors have been identified at the end of the selection process and derived for all the other basins (Table 2): the mean annual precipitation (MAP in mm/year); a seasonal precipitation index (p) given by the ratio between the mean precipitation of August (i.e., usually the least rainy month in Sicily) and the MAP; the drainage area (A in km 2 ); the mean areal Curve Number (CN); the main stream length (L in km); the average altitude (H in m a.s.l.); and the basin relief (z in m), given by the difference between the maximum and minimum elevation of the basin.In particular, parameters representative of the pluviometric regime (MAP and p) have been directly derived from the OA-ARRA database; the topographic features (A, L, H, z) have been derived by the DEM (Figure 2a), while the mean areal Curve Number (CN) was from the Regional Map of CN (Figure 2b).
The basin descriptor set obtained is consistent with other works concerning the regionalization of hydrologic models in Sicily.For example: in Cutore et al. [47], the same attributes H and L were considered together with the permeable soil percentage, which could be considered as a proxy of the variable CN; Viola et al. [55] considered A, MAP, H, CN, and the permeable soil percentage of each basin (not considered here); Cannarozzo et al. [56] used the variables A, L, H, and MAP.
In this study, Sicily as a whole has been initially assumed to be a homogeneous hydro-climatic region.As will be discussed in the following, this kind of approach has successively highlighted the need to consider a more efficient subdivision of Sicily into smaller homogeneous regions.For this purpose, three different sub-zones (A, B ,and C) within Sicily have been considered (Figure 3), using the same subdivisions suggested by Cannarozzo et al. [57] and adopted in Viola et al. [55].Such basin attributes, first examined for cross-correlations, were then individually analyzed in terms of their ability to explain the spatial variability of the basin's hydrologic response; more specifically, the Pearson correlation coefficient has been computed for all the possible combinations of basin attributes (independent variable) and three statistical values (mean, standard deviation, asymmetry) relative to the monthly observed runoff (dependent variable).The least significant explanatory variables have been dropped, selecting only basin attributes with sufficiently high coefficients (over 0.3, significant at 0.05 level) for at least one of the three examined statistics.
Seven basin descriptors have been identified at the end of the selection process and derived for all the other basins (Table 2): the mean annual precipitation (MAP in mm/year); a seasonal precipitation index (p) given by the ratio between the mean precipitation of August (i.e., usually the least rainy month in Sicily) and the MAP; the drainage area (A in km 2 ); the mean areal Curve Number (CN); the main stream length (L in km); the average altitude (H in m a.s.l.); and the basin relief (z in m), given by the difference between the maximum and minimum elevation of the basin.In particular, parameters representative of the pluviometric regime (MAP and p) have been directly derived from the OA-ARRA database; the topographic features (A, L, H, z) have been derived by the DEM (Figure 2a), while the mean areal Curve Number (CN) was from the Regional Map of CN (Figure 2b).
The basin descriptor set obtained is consistent with other works concerning the regionalization of hydrologic models in Sicily.For example: in Cutore et al. [47], the same attributes H and L were considered together with the permeable soil percentage, which could be considered as a proxy of the variable CN; Viola et al. [55] considered A, MAP, H, CN, and the permeable soil percentage of each basin (not considered here); Cannarozzo et al. [56] used the variables A, L, H, and MAP.
In this study, Sicily as a whole has been initially assumed to be a homogeneous hydro-climatic region.As will be discussed in the following, this kind of approach has successively highlighted the need to consider a more efficient subdivision of Sicily into smaller homogeneous regions.For this purpose, three different sub-zones (A, B ,and C) within Sicily have been considered (Figure 3), using the same subdivisions suggested by Cannarozzo et al. [57] and adopted in Viola et al. [55].
Among the 71 available gauged basins, a total of 18 basins have been excluded from calibration due to different factors.More specifically, 12 basins have been cautiously neglected: some because their runoff series could be affected by regulation due to the upstream presence of artificial reservoirs or diversion dams and, consequently, they cannot be considered representative of natural runoff conditions; and others because some missing data and outliers values have been detected in their hydro-climatic data series.Moreover, a total of six representative basins have been explicitly hidden for validation purposes, selecting a couple of basins, one small (i.e., A < 60 km 2 ) and the other larger (i.e., A > 250 km 2 ), for each of the three considered sub-zones (Figure 4).

Model Description and Assumptions
The proposed methodology is based on a two-step approach, which refers to two different phases of the calibration procedure.The first step includes the definition of a simple model structure for the creation of regression rainfall-runoff models performing well across all the basins within the study area.
The rainfall-runoff transformation process is usually characterized by complex nonlinear relations, for which the consideration of other interrelated variables could prove to be extremely useful.When snowmelt contributions can be neglected (e.g., arid and semi-arid regions), precipitation occurring during a certain month represents the main water input that, after a series of

Model Description and Assumptions
The proposed methodology is based on a two-step approach, which refers to two different phases of the calibration procedure.The first step includes the definition of a simple model structure for the creation of regression rainfall-runoff models performing well across all the basins within the study area.
The rainfall-runoff transformation process is usually characterized by complex nonlinear relations, for which the consideration of other interrelated variables could prove to be extremely useful.When snowmelt contributions can be neglected (e.g., arid and semi-arid regions), precipitation occurring during a certain month represents the main water input that, after a series of different hydrological processes, is partially released during the same month as runoff; as a consequence, the runoff results are often significantly correlated with rainfall.In addition to the precipitation, another variable that plays an important role in runoff modeling is the potential evapotranspiration, which represents a control factor for water losses from the soil due to evaporation and transpiration.In practical applications, temperature is often used as a proxy for the potential evapotranspiration [38,47,58], since this is rarely directly measured.Finally, the water accumulated in the soil (both saturated and unsaturated) could influence the slow response of the basins, providing further contributions to runoff generation due to the mechanisms of subsurface and groundwater runoff.This could explain the significant autocorrelation often observed in the runoff time series of basins in arid and semi-arid regions.
A backward analysis of the entire available database has confirmed the conceptual scheme mentioned above: runoff was highly correlated with precipitation and, to a lesser extent, with temperature for all the basins; moreover, in most of the basins, the runoff process was significantly autoregressive.Thus, the following regression-based rainfall-runoff model, suitable for many non-humid regions and already used in a similar way by Cutore et al. [47], has been assumed: where k refers to the analyzed basin, Q k (t) and Q k (t ´1) are the monthly natural streamflow values (mm) at the months t and t ´1, respectively, P k (t) and T k (t) are the precipitation (mm) and the mean temperature ( ˝C), respectively, at month t, and, finally, a i,k (with i = 1, 2, 3, 4) are four site-specific model parameters.The first calibration step requires us to individually calibrate the models, whose general form is given by Equation ( 1), for all the gauged basins across the region, finding, for each basin, the set of optimal parameters a i,k associated to the best model reproduction of the observed runoff series.
The second step consists in a regression analysis aimed at identifying a set of four possible regional relationships among such optimal parameters and some, or all, of the previously selected basin descriptors according to the following generic expression: where f i (with i = 1, 2, 3, 4) refers to the generic analytic function describing the regression model found for parameter a i , for which m regional regressive parameters (i.e., b i,g with i = 1, 2, 3, 4 and g = 1, . . .,m) have to be evaluated.The four models of Equation ( 2) may have a different structure and, consequently, a different number of regional parameters to be estimated for each of them.In particular, the various basin descriptors could have a different weight in explaining the a i optimal parameters' variability, and the identification, with consequent exclusion, of those poorly significant could be extremely useful in reducing the number of model parameters [23,24,38]; equifinality and overparameterization are, in fact, common problems for such kind of models and could reduce model robustness, introducing possible sources of model uncertainty.Initially, an attempt to fit a unique model for the entire region was carried out, trying to identify a unique set of equations characterizing the general form of Equation (2); nevertheless, as expected, better results were obtained by separately analyzing the three homogeneous sub-zones of Figure 3 and determining a different set of equations for each of them.

Regression Rainfall-Runoff Models for the Calibration Basins: Step 1A
The two sequential phases of the calibration procedure are schematically represented in Figure 5.The first step is aimed at identifying the set of parameters a i,k maximizing the agreement between observed and simulated monthly streamflow series at each calibration basin.The simulated-observed agreement is measured by the Nash-Sutcliffe Efficiency, NSE [59], defined as: where Q emp,t and Q sim,t are the observed (empirical) and the simulated (theoretical) streamflow at the time t, respectively, while Q emp is the average observed streamflow over the entire sample of observations whose size is N.The NSE values can range between ´8 and 1 (perfect agreement), with positive values generally denoting "behavioral" modeling.Moriasi et al. [60] introduced a performance rating that, for watershed simulations at the monthly time scale, classifies model performances as "satisfactory" for NSE values between 0.50 and 0.65, "good" for NSE between 0.65 and 0.75, and "very good" for NSE over 0.75.
time t, respectively, while emp Q is the average observed streamflow over the entire sample of observations whose size is N.The NSE values can range between −∞ and 1 (perfect agreement), with positive values generally denoting "behavioral" modeling.Moriasi et al. [60] introduced a performance rating that, for watershed simulations at the monthly time scale, classifies model performances as "satisfactory" for NSE values between 0.50 and 0.65, "good" for NSE between 0.65 and 0.75, and "very good" for NSE over 0.75.The calibration process has been performed using a specific software tool in Matlab (based on the functions optimset and lsqcurvefit) for unconstrained nonlinear optimization that provides, for the k-th basin, the set of parameters ai,k corresponding to the highest NSE index (i.e., best agreement).The resulting index can be thought of as the maximum efficiency that the model of Equation ( 1) can reach at the level of a single basin.The parameters obtained at this calibration phase (named Step 1A) for the 53 different basins are shown in Figure 6, where the optimal values of a1,k, a2,k, a3,k, and a4,k are represented in four different panels, using different markers for each sub-zone.Despite the use of an unconstrained optimization procedure, the resulting parameters vary among the different basins in The calibration process has been performed using a specific software tool in Matlab (based on the functions optimset and lsqcurvefit) for unconstrained nonlinear optimization that provides, for the k-th basin, the set of parameters a i,k corresponding to the highest NSE index (i.e., best agreement).The resulting index can be thought of as the maximum efficiency that the model of Equation ( 1) can reach at the level of a single basin.The parameters obtained at this calibration phase (named Step 1A) for the 53 different basins are shown in Figure 6, where the optimal values of a 1,k , a 2,k , a 3,k , and a 4,k are represented in four different panels, using different markers for each sub-zone.Despite the use of an unconstrained optimization procedure, the resulting parameters vary among the different basins in very narrow ranges, with positive values for a 1 , a 2 , and a 4 and negative values for a 3 .Given the model structure used (see Equation ( 1)), this implies that streamflow, at a certain month, is positively correlated to rainfall in the same month and streamflow in the previous month, while it was negatively correlated to the mean temperature in the same month, as expected.
Water 2016, 8, 209 very narrow ranges, with positive values for a1, a2, and a4 and negative values for a3.Given the model structure used (see Equation ( 1)), this implies that streamflow, at a certain month, is positively correlated to rainfall in the same month and streamflow in the previous month, while it was negatively correlated to the mean temperature in the same month, as expected.

Regional Regression Analysis: Step 1B and Step 2
The second step (Step 2) aims to establish the functional relationships characterizing Equation (2) and relating the set of optimal model parameters to some basin attributes among those previously selected (i.e., MAPk, pk, Ak, Lk, Hk, zk, CNk).
A first attempt of regionalization has been carried out considering Sicily as a unique homogeneous region and using all seven basin descriptors as predictors in linear regression equations; nevertheless, none of the tested regressive models has provided satisfying coefficients of determination (highest R 2 equal to 0.26).A general and significant improvement of the same models has been observed considering the region's subdivision into the sub-zones of Figure 3, even if the value of R 2 associated with the regression models for the parameter a1,k remained rather low, especially for sub-zone A (R 2 = 0.33).For this reason, the value of a1 has been assumed as constant for each sub-zone and equal to the mean of all the values of a1,k obtained at Step 1A for each sub-zone.The optimizing procedure of Step 1 has been successively iterated with constant a1 for each sub-zone, finding a new set of optimal values for a2,k, a3,k, and a4,k for each basin (circles in Figure 7).The parameters assessed during this calibration phase (named Step 1B) are not significantly dissimilar from those obtained at Step 1A and exhibit very similar values in terms of mean, standard deviation, maximum, and minimum over the 53 considered basins.
At Step 2, different linear and nonlinear multiple regressions among the new parameter sets (Step 1B) and the different basin descriptors are analyzed with the objective of finding a unique (i.e., the same structure in the different sub-zones), simple (i.e., few parameters), and efficient (i.e., high R 2 ) regression model for each parameter of Equation ( 1).An important aspect in multi-regressive analysis is the evaluation of the relative importance of each regressor (i.e., basin descriptors) in explaining the dependent variables' (i.e., coefficients ai,k) variance, which allows for identifying a set of only significant covariates, reducing model parameters and, thus, possible sources of uncertainty.To this end, for each parameter ai,k and for two possible regression model forms (linear or log-linear), the following stepwise procedure has been adopted: (1) several models are set by adding, in successive steps, a certain predictor to all the possible subsets of the remaining

Regional Regression Analysis: Step 1B and Step 2
The second step (Step 2) aims to establish the functional relationships characterizing Equation ( 2) and relating the set of optimal model parameters to some basin attributes among those previously selected (i.e., MAP k , A first attempt of regionalization has been carried out considering Sicily as a unique homogeneous region and using all seven basin descriptors as predictors in linear regression equations; nevertheless, none of the tested regressive models has provided satisfying coefficients of determination (highest R 2 equal to 0.26).A general and significant improvement of the same models has been observed considering the region's subdivision into the sub-zones of Figure 3, even if the value of R 2 associated with the regression models for the parameter a 1,k remained rather low, especially for sub-zone A (R 2 = 0.33).For this reason, the value of a 1 has been assumed as constant for each sub-zone and equal to the mean of all the values of a 1,k obtained at Step 1A for each sub-zone.The optimizing procedure of Step 1 has been successively iterated with constant a 1 for each sub-zone, finding a new set of optimal values for a 2,k , a 3,k , and a 4,k for each basin (circles in Figure 7).The parameters assessed during this calibration phase (named Step 1B) are not significantly dissimilar from those obtained at Step 1A and exhibit very similar values in terms of mean, standard deviation, maximum, and minimum over the 53 considered basins.
At Step 2, different linear and nonlinear multiple regressions among the new parameter sets (Step 1B) and the different basin descriptors are analyzed with the objective of finding a unique (i.e., the same structure in the different sub-zones), simple (i.e., few parameters), and efficient (i.e., high R 2 ) regression model for each parameter of Equation ( 1).An important aspect in multi-regressive analysis is the evaluation of the relative importance of each regressor (i.e., basin descriptors) in explaining the dependent variables' (i.e., coefficients a i,k ) variance, which allows for identifying a set of only significant covariates, reducing model parameters and, thus, possible sources of uncertainty.To this end, for each parameter a i,k and for two possible regression model forms (linear or log-linear), the following stepwise procedure has been adopted: (1) several models are set by adding, in successive steps, a certain predictor to all the possible subsets of the remaining predictors; (2) each model is fitted for each sub-zone by the ordinary least squares method; (3) the mean value of the coefficient of determination R 2 over the three sub-zones is associated to each model; (4) the relative weight of that predictor is given by the average R 2 over all the tested models; (5) the procedure is iterated for all seven basin descriptors; (6) descriptors explaining less than 5% of the dependent variables variance are rejected, while the remaining ones are considered.At the end of this procedure, the following linear regressive models have been obtained: where the symbol Z refers to the sub-zone (Z = A, B, or C) of the k-th basin, while b i,j,Z (with i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5) are 16 regional regression parameters, estimated for each sub-zone by the least squares method and reported in Table 3.With this structure, the three parameters a 2,k , a 3,k , and a 4,k of the rainfall-runoff model of Equation ( 1) are all dependent on the basin mean altitude and mean areal curve number (H k and CN k , respectively) and other two basin descriptors; more specifically, the exponents for the climatic variables P k (t) and T k (t) in Equation (1) (i.e., a 2,k and a 3,k , respectively) are both also related to the basin pluviometric regime descriptors MAP k and p k , while the exponent (i.e., a 4,k ) for the previous month's runoff, Q k (t ´1), is derived as a function of the basin drainage area (A k ) and the relief (z k ), two basin attributes that influence the streamflow transferring process toward the basin outlet and whose variability among the different basins could influence the lag-1-month autocorrelation in the corresponding streamflow series.
Water 2016, 8, 209 12 of 24 predictors; (2) each model is fitted for each sub-zone by the ordinary least squares method; (3) the mean value of the coefficient of determination R 2 over the three sub-zones is associated to each model; (4) the relative weight of that predictor is given by the average R 2 over all the tested models; (5) the procedure is iterated for all seven basin descriptors; (6) descriptors explaining less than 5% of the dependent variables variance are rejected, while the remaining ones are considered.At the end of this procedure, the following linear regressive models have been obtained: where the symbol Z refers to the sub-zone (Z = A, B, or C) of the k-th basin, while bi,j,Z (with i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5) are 16 regional regression parameters, estimated for each sub-zone by the least squares method and reported in Table 3.With this structure, the three parameters a2,k, a3,k, and a4,k of the rainfall-runoff model of Equation ( 1) are all dependent on the basin mean altitude and mean areal curve number (Hk and CNk, respectively) and other two basin descriptors; more specifically, the exponents for the climatic variables Pk(t) and Tk(t) in Equation (1) (i.e., a2,k and a3,k, respectively) are both also related to the basin pluviometric regime descriptors MAPk and pk, while the exponent (i.e., a4,k) for the previous month's runoff, Qk(t − 1), is derived as a function of the basin drainage area (Ak) and the relief (zk), two basin attributes that influence the streamflow transferring process toward the basin outlet and whose variability among the different basins could influence the lag-1-month autocorrelation in the corresponding streamflow series.The regionalized values of a 2,k , a 3,k , and a 4,k for all the calibration basins (model parameters at Step 2), obtained by Equation ( 4) with the regression parameters of Table 3, are also reported in Figure 7 and compared with the parameter values resulting after Step 1B.The parameters at Step 1B (optimal) and Step 2 (regionalized) have almost the same mean value over the 53 basins, with differences below 3%, while the variance of the regionalized parameters is reduced by about 16% for a 2,Z and a 3,Z and by 34% for a 4,Z .The most significant variation between the parameters before and after regionalization can be noticed for parameter a 3,Z , with a maximum at the basin with ID 59 (sub-zone C) where the regionalized value results are one order of magnitude higher than the corresponding optimal parameter at Step 1B.
Table 3. Regional regression parameters b i,j,Z for the estimation of the parameters a i,k by Equation (4).

Model Performances for the Calibration Basins
The model in Equation ( 1), with the four parameters individually calibrated for each site (Step 1A, Figure 6), has shown, over the 53 calibration basins, a high capacity for reproducing observed streamflows, with NSE values ranging from 0.64 to 0.93 and a mean of 0.77 (Table 4a).The NSE values for such basins are represented in the histogram of Figure 7 (upper-left panel) and in the map of Figure 8a.From the latter as well as the results synthesized in Table 4a, it can be observed that the model in sub-zone C provides, on average, higher performance than results in the other two sub-zones.
The assumption that a 1 is constant for each sub-zone has provided new sets of optimal values for a 2,k , a 3,k , and a 4,k (Step 1B) whose corresponding performance, in terms of NSE, is reported in the same histogram of Figure 7 and Table 4a.The results are quite similar to those previously discussed (Step 1A), denoting a not significant decrease in the model's performances; on average, a reduction in the order of 2% for the NSE indexes can be noticed.A significant NSE reduction (from 0.73 to 0.44) has been detected only for the ID 66 basin (sub-zone C), which is characterized by the lowest value of MAP and the second highest value of p.As expected, the use of the regional relationships in Equation ( 4) for the a i,k parameter assessment (Step 2), reduces model performance with respect to the adoption of optimal (Steps 1A and B) parameter sets (histogram of Figure 7).Actually, this performance reduction is a price that has to be paid in order to provide a procedure able to assess runoff at ungauged sites.Nevertheless, the final model performances (after the Step 2) (Table 4a) remain in the range of acceptable values in terms of NSE, i.e., performance equal or more than satisfactory, according to the performance rating of Moriasi et al. [60], for all the calibration basins; the only two exceptions are the abovementioned ID 66 basin (sub-zone C), which had the lowest starting NSE value (i.e., NSE = 0.44 at Step 1B), and ID 22 basin (sub-zone A), where the resulting values are, however, positive and equal to 0.33 and 0.35, respectively.For this last basin, model performance reduction after the regionalization procedure is rather consistent, since, after Step 1B, the NSE was quite high (0.77).This could be attributable to the fact that the parameter values derived for this basin after Step 2 are markedly different from the corresponding values obtained at Step 1B, especially with regard to the parameter a 4 , whose regionalized value results were halved with respect to the previously obtained optimal value (from 0.43 to 0.23).the final model performances (after the Step 2) (Table 4a) remain in the range of acceptable values in terms of NSE, i.e., performance equal or more than satisfactory, according to the performance rating of Moriasi et al. [60], for all the calibration basins; the only two exceptions are the abovementioned ID 66 basin (sub-zone C), which had the lowest starting NSE value (i.e., NSE = 0.44 at Step 1B), and ID 22 basin (sub-zone A), where the resulting values are, however, positive and equal to 0.33 and 0.35, respectively.For this last basin, model performance reduction after the regionalization procedure is rather consistent, since, after Step 1B, the NSE was quite high (0.77).This could be attributable to the fact that the parameter values derived for this basin after Step 2 are markedly different from the corresponding values obtained at Step 1B, especially with regard to the parameter a4, whose regionalized value results were halved with respect to the previously obtained optimal value (from 0.43 to 0.23).The percent reduction of the NSE after the regionalization, with respect to the NSE computed at the Step 1B, ranges from 0.17% to 57%, with the mean NSE (=0.69) over all the calibration basins reduced by almost 9% and denoting overall good performances.The efficiency reductions after the regionalization with respect to Step 1A are slightly more marked, with NSE reduced by 11% on average.This comparison is also emphasized in Figure 8, where the NSE obtained with the regionalized parameters (Figure 8b) and with the parameter sets of Step 1A (Figure 8a) are represented by the same color bar.Also, after the regionalized procedure, the model exhibits, on average, a higher accuracy for the basins in sub-zone C, with a mean NSE of 0.70 (Table 4a).The ID 47 basin (sub-zone C) has shown the highest NSE at every calibration phase (Table 4a), with a performance reduction after the regionalization (Step 2) equal to about 3%.The percent reduction of the NSE after the regionalization, with respect to the NSE computed at the Step 1B, ranges from 0.17% to 57%, with the mean NSE (=0.69) over all the calibration basins reduced by almost 9% and denoting overall good performances.The efficiency reductions after the regionalization with respect to Step 1A are slightly more marked, with NSE reduced by 11% on average.This comparison is also emphasized in Figure 8, where the NSE obtained with the regionalized parameters (Figure 8b) and with the parameter sets of Step 1A (Figure 8a) are represented by the same color bar.Also, after the regionalized procedure, the model exhibits, on average, a higher accuracy for the basins in sub-zone C, with a mean NSE of 0.70 (Table 4a).The ID 47 basin (sub-zone C) has shown the highest NSE at every calibration phase (Table 4a), with a performance reduction after the regionalization (Step 2) equal to about 3%.
Model performances at Step 2 have also been assessed using three further statistical criteria measuring the agreement between observed and simulated monthly streamflow series: the mean error, ME (in mm/month); the dimensionless mean error, ME/ Q emp , given by the ratio between the ME and the average observed streamflow (Q emp ); and the root mean square error, RMSE (in mm/month).For each index, the mean, the standard deviation, the minimum, the maximum, and the best values over all the calibration basins are synthesized in Table 4b.The results of this analysis have confirmed the outcomes relative to the performance index used for calibration (i.e., NSE), confirming the validity of the adopted calibration procedure.The model has globally provided a satisfying accuracy in terms of all the analyzed indexes and can be considered unbiased, as demonstrated by the relatively low values of ME.The values of ME, in absolute value, are lower than 3 mm/month for about 80% of the analyzed basins, while the mean ME/ Q emp over the calibration basins is almost null.In terms of RMSE, the indexes have been of the same order of magnitude as those found through a similar modeling approach by Cutore et al. [47] for the Simeto river sub-basins (Sicily, Italy), with the worst performance (about 45 mm/month) obtained at basins ID 61, 67, and 68, which are all basins corresponding to NSE values lower than the mean value (0.69, 0.63, and 0.58, respectively).
The basin with the highest NSE (i.e., ID 47) also provided satisfactory performance in terms of ME, ME/ Q emp , and RMSE (0.428, 0.042, and 5.05 mm/month, respectively).The basins with the lowest NSE (i.e., ID 22 and 66), showed the worst performance with respect to the other considered statistical criteria (ME = ´6.67 and ´2.54 mm/month, ME/ Q emp = ´0.49and ´0.69, RMSE = 15.37 and 9.97 mm/month, respectively).Only one basin (ID 40; sub-zone C) associated with a satisfying NSE (0.63), even if it was below the mean of the other basins, has provided relevant values for ME (11 mm/month) and ME/ Q emp (0.35), denoting a weak performance that is, however, comparable to the worst performance obtained in different studies (e.g., [47]).
A comparison between monthly flow duration curves (FDCs) based on observed and simulated runoff values is reported in Figure 9 for different representative basins.Four examples, referring to basins with performance indexes approximately equal to the averages over all the calibration basins, are reported: basin ID 18 can be considered representative of basins with mean NSE, basin ID 42 is representative of basins with mean ME, basin ID 25 is representative of basins with mean ME/ Q emp , while basin ID 59 is representative of basins with mean RMSE.Moreover, since the FDCs typically depend on the period considered, a comparison of the basins with the shortest (ID 45) and the longest (ID 23) sample size is also reported.Although some important details of the variations in flows can be obscured when the FDCs are computed using monthly data, rather than daily or finer resolution data, this analysis has been useful in confirming model accuracy.For all the examined basins, in fact, the model has satisfactorily reproduced the magnitudes associated with the various observed monthly runoff values, from the lowest and more frequent values to the highest and rarer ones, proving to be effective also in the assessment of the probability that a certain runoff value will be equaled or exceeded.
Despite the fact that the validation basins are characterized by marked differences in terms of rainfall-runoff transformation, with mean annual runoff coefficients ranging from 0.14 (ID 43) to 0.34 (ID 49), the model has shown, over the six basins, an equal ability to capture the different basins' hydrological response, with performance that can be classified as "good" in all the basins and, for two cases (i.e., basins ID 24 and 46), even as "very good" (i.e., NSE > 0.75).
Figure 10 depicts a comparison between observed and simulated monthly specific streamflow series (mm/month) for the six validation basins, also reporting the corresponding precipitation series.Simulated series are quite close to the observed series, reproducing well most of the peaks and null values.Despite the marked differences, in terms of both observed rainfall and streamflow, that can be noticed among the basins, the model captures rather accurately the monthly streamflow variability for both the smaller (left panels) and larger (right panels) basins in the three sub-zones.For example, the ID 12 basin is characterized by a seasonal streamflow regime, with about five months per year almost dry and frequent winter peaks with streamflow of about 100 mm/month, while the ID 24 basin is characterized by more regular behavior in the observed streamflow, with, on average, less rainfall and streamflow, and with only five months out of eight years having streamflow on the order of 100 mm/month; for both basins the model showed similar performance.Despite the fact that the validation basins are characterized by marked differences in terms of rainfall-runoff transformation, with mean annual runoff coefficients ranging from 0.14 (ID 43) to 0.34 (ID 49), the model has shown, over the six basins, an equal ability to capture the different basins' hydrological response, with performance that can be classified as "good" in all the basins and, for two cases (i.e., basins ID 24 and 46), even as "very good" (i.e., NSE > 0.75).
Figure 10 depicts a comparison between observed and simulated monthly specific streamflow series (mm/month) for the six validation basins, also reporting the corresponding precipitation series.Simulated series are quite close to the observed series, reproducing well most of the peaks and null values.Despite the marked differences, in terms of both observed rainfall and streamflow, that can be noticed among the basins, the model captures rather accurately the monthly streamflow variability for both the smaller (left panels) and larger (right panels) basins in the three sub-zones.For example, the ID 12 basin is characterized by a seasonal streamflow regime, with about five months per year almost dry and frequent winter peaks with streamflow of about 100 mm/month, while the ID 24 basin is characterized by more regular behavior in the observed streamflow, with, on average, less rainfall and streamflow, and with only five months out of eight years having streamflow on the order of 100 mm/month; for both basins the model showed similar performance.

Model Performances at Different Aggregated Time Scales
Model performances have been further evaluated at different temporal aggregations, also considering the seasonal and the annual time scales.This analysis has been performed with regard to both the 53 calibration basins and the six validation basins, considered as a unique sample.More specifically, simulated monthly streamflow has been aggregated at the annual scale and also at the seasonal scale, considering the year as divided into two seasons: dry season, from April to September, and wet season, the remaining six months of the year.
Figure 11 compares all the estimates of the monthly, annual, and seasonal streamflow obtained from the regional model for all the basins with the corresponding observed streamflow.The high predictive ability of the model at the different aggregation scales is demonstrated by the fact that, for all the plots, most of the points are rather close to the perfect agreement lines (also reported in all the graphs), with high values of the coefficient of determination R 2 .These values are greater than 0.92 at the monthly level and 0.90 at the annual level.At the seasonal level, the model reproduces the dry season streamflow with an R 2 of 0.95 and the wet season values with an R 2 of 0.90.Although the cloud of points appears to be more disperse than in the other plots, the high value of R 2 obtained at the monthly scale can be explained by the presence of a considerable number of observed values that are identically reproduced by the model (i.e., 22% of the null values and about 5% of the not-null streamflow).

Model Performances at Different Aggregated Time Scales
Model performances have been further evaluated at different temporal aggregations, also considering the seasonal and the annual time scales.This analysis has been performed with regard to both the 53 calibration basins and the six validation basins, considered as a unique sample.More specifically, simulated monthly streamflow has been aggregated at the annual scale and also at the seasonal scale, considering the year as divided into two seasons: dry season, from April to September, and wet season, the remaining six months of the year.The high predictive ability of the model at the different aggregation scales is demonstrated by the fact that, for all the plots, most of the points are rather close to the perfect agreement lines (also reported in all the graphs), with high values of the coefficient of determination R 2 .These values are greater than 0.92 at the monthly level and 0.90 at the annual level.At the seasonal level, the model reproduces the dry season streamflow with an R 2 of 0.95 and the wet season values with an R 2 of 0.90.Although the cloud of points appears to be more disperse than in the other plots, the high value of R 2 obtained at the monthly scale can be explained by the presence of a considerable number of observed values that are identically reproduced by the model (i.e., 22% of the null values and about 5% of the not-null streamflow).

Conclusions
This study applies a regional regressive model for the reconstruction of natural runoff series at the monthly or coarser time resolutions, intensively analyzing its performance.The importance of developing appropriate and easily transferable rainfall-runoff models able to reproduce reliable and long runoff series from climatic data is related to the fact that, in many practical cases, a runoff dataset is not available or insufficient in size, while long data series of climatic variables are much more often available, and, then, can be exploited to fill gaps in runoff records.A practical means of estimating natural runoff dynamics at ungauged stream sections or for rivers subject to regulation is fundamental for a variety of modern hydrological applications that involve a number of sectors (e.g., engineering, economic, environmental, etc.).
The idea proposed in this paper is to implement a calibration procedure for rainfall-runoff models based on the knowledge of "soft" basins information, usually available, such as rainfall, temperature, morphology, or land use.The model has been calibrated, through a simple two-step procedure, for the island of Sicily using a considerable number of basins and it has been validated in six basins representative of different basin sizes and different climatic subareas within the examined region.
This application may be considered as a benchmark for similar hydrological and climatic areas.Through a detailed description of the adopted methodology, some practical aspects that may be of interest for the model's application to other regions have been explored.For instance, in this study, model performance has been revealed to be higher if the region is considered as divided into three different sub-zones.Moreover, six different basin attributes have been demonstrated to be functional in the definition of opportune regional equations for the estimation of the rainfall-runoff model parameters.The same methodology, applied to arid and semi-arid regions different from the study area of this work, could converge towards a different subdivision of the area of interest and/or different regional equations, characterized by different basin descriptors; nevertheless, the adopted structure for the rainfall-runoff models, which has shown extremely high accuracy at all Sicilian basins before the regionalization step, is expected to be equally valid.On the contrary, the application to humid regions should also encompass a redefinition of a more suitable conceptual scheme for the regression model, even if the proposed methodological approach could equally be followed.
For all the basins analyzed in this study, the regionalized model has shown satisfying accuracy, measured by different performance indexes, in reproducing the observed runoff at different time scales from monthly to aggregated seasonal and annual scales.Considering the Nash-Sutcliffe Efficiency (NSE) for the monthly series reproduction, it was rather weak for only two cases, but nevertheless positive (>0.33), denoting a "behavioral" modeling, while the NSE indexes for the other basins were all over 0.50, with "satisfactory" performance for 21% of cases, "good" for 53%, and "very good" for the remaining 26%, according to the performance criteria given by Moriasi et al. [60].At coarser time resolutions, model performance did not show significant variations, with a very low percentage of error in the reproduction of both the seasonal and the annual mean runoff in any basin and an ability to also capture the seasonal and interannual runoff variability.The robustness of the model is strongly confirmed by the results in validation, in which a model reproduction of the historical runoff series at all the tested time scales gave results comparable to those obtained in calibration.
The obtained results have demonstrated how the model discussed in this study can be seen as a prompt and suitable option for the estimation of natural runoff time series at ungauged sites.It is worth emphasizing that an important application of the model is currently under development.In particular, a specific plug-in is being implemented within open-source GIS software (i.e., Quantum GIS 2.10) with the aim of automatizing data retrieval and processing procedures.After a preliminary preparation of an opportune Database Management System (DBMS) supporting the plug-in, the regional model, entirely re-coded in Python scripting language, will be able to rapidly assess natural runoff series (at monthly or higher time scales) for any basin of Sicily and any desired time window, both selected by an opportunely implemented Graphical User Interface (GUI).The plug-in will be based on a Data Processing Module (DPM), which, using spatial analysis techniques, will first proceed with the watershed classification (i.e., into sub-zones) and delineation, derivation of the monthly time series of mean areal precipitation and temperature, extraction of the six basin descriptors needed for the model and, then, derivation of the monthly runoff time series.

Figure 1 .
Figure 1.Study area (Sicily, Italy).(a) Location of the analyzed river basins with position of the relative hydrometric gauging stations.The hydrographic network is also highlighted; (b) Location and typology of the meteo-climatic stations.

Figure 1 .
Figure 1.Study area (Sicily, Italy).(a) Location of the analyzed river basins with position of the relative hydrometric gauging stations.The hydrographic network is also highlighted; (b) Location and typology of the meteo-climatic stations.

Figure 2 .
Figure 2. (a) Elevations (m a.s.l.) from the considered Digital Elevation Model (DEM) of Sicily; (b) spatial distribution of the CN parameter (raster map by Regione Sicilia).

Figure 2 .
Figure 2. (a) Elevations (m a.s.l.) from the considered Digital Elevation Model (DEM) of Sicily; (b) spatial distribution of the CN parameter (raster map by Regione Sicilia).

Figure 5 .
Figure 5. Schematic representation of the different phases of model calibration and validation.

Figure 5 .
Figure 5. Schematic representation of the different phases of model calibration and validation.

Figure 6 .
Figure 6.Optimal calibration parameters (a1, a2, a3, and a4) after Step 1A for all the calibration basins.Different markers are used for each sub-zone.

Figure 6 .
Figure 6.Optimal calibration parameters (a 1 , a 2 , a 3 , and a 4 ) after Step 1A for all the calibration basins.Different markers are used for each sub-zone.

Figure 7 .
Figure 7. Parameters and model performances at the different calibration phases: optimal calibration parameters (a2, a3, and a4) after Step 1B and regionalized parameters after Step 2 for all the calibration basins; in the upper-left histogram, the resulting NSE values after each calibration phase (Step 1A, Step 1B, and Step 2).

Figure 7 .
Figure 7. Parameters and model performances at the different calibration phases: optimal calibration parameters (a 2 , a 3 , and a 4 ) after Step 1B and regionalized parameters after Step 2 for all the calibration basins; in the upper-left histogram, the resulting NSE values after each calibration phase (Step 1A, Step 1B, and Step 2).

Figure 10 .
Figure 10.Observed (marks) and simulated (dashed lines) monthly streamflow series (mm/month) for the validation basins, with histograms denoting the underlying precipitation series.Only the first eight years are reported.Upper panels refer to the basins of sub-zone A, middle panels to the basins in sub-zone B, while lower panels refer to sub-zone C. Left panels refer to the smaller basins (A < 60 km 2 ), while right panels refer to the largest ones (A > 250 km 2 ).

Figure 10 .
Figure 10.Observed (marks) and simulated (dashed lines) monthly streamflow series (mm/month) for the validation basins, with histograms denoting the underlying precipitation series.Only the first eight years are reported.Upper panels to the basins of sub-zone A, middle panels to the basins in sub-zone B, while lower panels refer to sub-zone C. Left panels refer to the smaller basins (A < 60 km 2 ), while right panels refer to the largest ones (A > 250 km 2 ).
In the top left panel of the figure, a total of 12,312 theoretical monthly streamflow estimates are plotted against the corresponding empirical values, while the other scatter diagrams on the top refer to the annual streamflow values (middle panel, 1026 values) and the mean annual totals of streamflow for each basin (right panel, 59 values).In the bottom panels of Figure 11, the empirical and theoretical dry season (left panel), wet season (middle panel) streamflow values and the mean seasonal totals of streamflow for each basin (right panel) are similarly compared.

Figure 11
compares all the estimates of the monthly, annual, and seasonal streamflow obtained from the regional model for all the basins with the corresponding observed streamflow.In the top left panel of the figure, a total of 12,312 theoretical monthly streamflow estimates are plotted against the corresponding empirical values, while the other scatter diagrams on the top refer to the annual streamflow values (middle panel, 1026 values) and the mean annual totals of streamflow for each basin (right panel, 59 values).In the bottom panels of Figure 11, the empirical and theoretical dry season (left panel), wet season (middle panel) streamflow values and the mean seasonal totals of streamflow for each basin (right panel) are similarly compared.

Figure 11 .
Figure 11.Model performance at the monthly, seasonal, and annual scale for all the calibration and validation basins.In the right two panels, the empirical mean annual (top) and seasonal (bottom) streamflow totals for each basin are plotted against the corresponding theoretical values.Dry season is the six months from April to September, while the Wet season is the remaining part of the year.Dashed lines indicate the perfect agreement between observed and simulated values.

Figure 11 .
Figure 11.Model performance at the monthly, seasonal, and annual scale for all the calibration and validation basins.In the right two panels, the empirical mean annual (top) and seasonal (bottom) streamflow totals for each basin are plotted against the corresponding theoretical values.Dry season is the six months from April to September, while the Wet season is the remaining part of the year.Dashed lines indicate the perfect agreement between observed and simulated values.

Table 1 .
Watersheds and hydrometric gauge station names (St.Name); identification code (ID); typology (Type): main basin (MB) or sub-basin (SB) and sub-zone (A, B, or C); observation period: first and last year of observation and total sample size in years (yr).

Table 2 .
Basin descriptors: MAP = mean annual precipitation; p = seasonal precipitation index; A = drainage area; L = main river length; H = mean basin altitude; z = basin relief; CN = Curve Number.

Table 4 .
(a) Analysis of NSE indexes (mean for sub-zone, and mean, standard deviation, minimum, and maximum over all the calibration basins) at the different calibration stages (Step 1A, Step 1B, Step 2); (b) Analysis (mean, standard deviation, minimum, and maximum) of other performance indexes (ME = mean error; ME/Q emp = dimensionless mean error; RMSQ = root mean square error) over the calibration basins after the regionalization (Step 2).The basin ID codes corresponding to the best, minimum, and maximum performance indexes are also reported within brackets.

Table 4 .
(a) Analysis of NSE indexes (mean for sub-zone, and mean, standard deviation, minimum, and maximum over all the calibration basins) at the different calibration stages (Step 1A, Step 1B, Step 2); (b) Analysis (mean, standard deviation, minimum, and maximum) of other performance indexes (ME = mean error; ME/ emp Q = dimensionless mean error; RMSQ = root mean square error) over the calibration basins after the regionalization (Step 2).The basin ID codes corresponding to the best, minimum, and maximum performance indexes are also reported within brackets.(a)Mean NSE for Sub-Zone NSE (53 Calibration Basins)