Conclusions References

Abstract. A Regional Water Resources study was performed at basins within and draining to the Basque Country Region (N of Spain), with a total area of approximately 8500 km2. The objective was to obtain daily and monthly long-term discharges in 567 points, most of them ungauged, with basin areas ranging from 0.25 to 1850 km2. In order to extrapolate the calibrations at gauged points to the ungauged ones, a distributed and conceptually based model called TETIS was used. In TETIS the runoff production is modelled using five linked tanks at the each cell with different outflow relationships at each tank, which represents the main hydrological processes as snowmelt, evapotranspiration, overland flow, interflow and base flow. The routing along the channels' network couples its geomorphologic characteristics with the kinematic wave approach. The parameter estimation methodology tries to distinguish between the effective parameter used in the model at the cell scale, and the watershed characteristic estimated from the available information, being the best estimation without losing its physical meaning. The relationship between them can be considered as a correction function or, in its simple form, a correction factor. The correction factor can take into account the model input errors, the temporal and spatial scale effects and the watershed characteristics. Therefore, it is reasonable to assume the correction factor is the same for each parameter to all cells within the watershed. This approach reduces drastically the number of parameter to be calibrated, because only the common correction factors are calibrated instead of parameter maps (number of parameters times the number of cells). In this way, the calibration can be performed using automatic methodologies. In this work, the Shuffled Complex Evolution – University of Arizona, SCE-UA algorithm was used. The available recent year's data was used to calibrate the model in 20 of the most representative flow gauge stations in 18 basins with a Nash-Sutcliffe index higher than 0.6 (10 higher than 0.8). The calibrated correction factors at each basin were similar but not equal. The validation process (in time and space) was performed using the remaining data in all flow gauge stations (62), with 42 basins with a Nash-Sutcliffe index higher than 0.5 (25 higher than 0.7). Deficient calibration and validations were always related with flow gauge stations very close to the karstic springs. These results confirmed that it was feasible and efficient to use the SCE-UA algorithm for the automatic calibration of distributed conceptual models and the calibrated model could be used at ungauged basins. Finally, meteorological information from the past 50 years at a daily scale was used to generate a daily discharges series at 567 selected points.

Abstract.A Regional Water Resources study was performed at basins within and draining to the Basque Country Region (N of Spain), with a total area of approximately 8500 km 2 .The objective was to obtain daily and monthly long-term discharges in 567 points, most of them ungauged, with basin areas ranging from 0.25 to 1850 km 2 .In order to extrapolate the calibrations at gauged points to the ungauged ones, a distributed and conceptually based model called TETIS was used.In TETIS the runoff production is modelled using five linked tanks at the each cell with different outflow relationships at each tank, which represents the main hydrological processes as snowmelt, evapotranspiration, overland flow, interflow and base flow.The routing along the channels' network couples its geomorphologic characteristics with the kinematic wave approach.The parameter estimation methodology tries to distinguish between the effective parameter used in the model at the cell scale, and the watershed characteristic estimated from the available information, being the best estimation without losing its physical meaning.The relationship between them can be considered as a correction function or, in its simple form, a correction factor.The correction factor can take into account the model input errors, the temporal and spatial scale effects and the watershed characteristics.Therefore, it is reasonable to assume the correction factor is the same for each parameter to all cells within the watershed.This approach reduces drastically the number of parameter to be calibrated, because only the common correction factors are calibrated instead of pa-Correspondence to: J. J. Vélez (jjvelezu@unal.edu.co)rameter maps (number of parameters times the number of cells).In this way, the calibration can be performed using automatic methodologies.In this work, the Shuffled Complex Evolution -University of Arizona, SCE-UA algorithm was used.The available recent year's data was used to calibrate the model in 20 of the most representative flow gauge stations in 18 basins with a Nash-Sutcliffe index higher than 0.6 (10 higher than 0.8).The calibrated correction factors at each basin were similar but not equal.The validation process (in time and space) was performed using the remaining data in all flow gauge stations (62), with 42 basins with a Nash-Sutcliffe index higher than 0.5 (25 higher than 0.7).Deficient calibration and validations were always related with flow gauge stations very close to the karstic springs.These results confirmed that it was feasible and efficient to use the SCE-UA algorithm for the automatic calibration of distributed conceptual models and the calibrated model could be used at ungauged basins.Finally, meteorological information from the past 50 years at a daily scale was used to generate a daily discharges series at 567 selected points.

Introduction: problem framework
During last decade hydrological institutions and research groups have focused in the Prediction on Ungauged Basins (PUB) as a research topic which must be studied, analyzed and restructured if necessary in this decade (Sivapalan et al., 2003).McDonnell et al. (2007) mention that to make progress in watershed hydrology with coherence it is necessary to characterize the landscape heterogeneity exploring the principles that might underlie the heterogeneity and Published by Copernicus Publications on behalf of the European Geosciences Union.complexity.Nevertheless, the ability to generalize these findings to ungauged regions remains unsolved.However, in this paper the extrapolation of the parameters to ungauged basins with a hydrological distributed model is successfully exposed.
Knowing water budgets for the basin and its subbasins allows water-resource managers to identify areas of greatest concern in the basin and make more informed decisions about conservation and remediation.A Regional Water Resources study was performed at basins within and draining to the Basque Country Region (N of Spain), with a total area of approximately 8500 km 2 .The objective was to obtain daily and monthly long-term discharges in 567 points (initially 123 points, but at the end of the project more simulation points were demanded), obviously most of them ungauged, with basin areas ranging from 0.25 to 1850 km 2 .The study area is officially divided into 21 hydrological systems with a total of 41 basins (Fig. 1), only 17 of them with at least one flow gauge station.Northern basins drain to the Bay of Biscay in the Atlantic Ocean with a humid climate, and southern basins drain to the Ebro River (which flows to the Mediterranean Sea) with a dry continental climate.
In this study, distributed modelling was proposed mainly: -Due to the high number of simulation points and its vague definition at the beginning of the project.
-In order to extrapolate the calibrations at gauged points to the ungauged ones, without parameters regionalization but reproducing the natural spatial variability of the Hydrological Cycle.
-Because spatial information is available to be used in this type of models.
Distributed models divide the catchment into a number of smaller areas (grid elements or subcatchments), which are assumed to be uniform with respect to their hydrologic parameters.There has been an "exponentially growing" recent interest in distributed hydrologic modelling that has been fuelled by growing availability of GIS-related information (Schaake, 2003).
Temporal and spatial heterogeneity can be important in many hydrologic applications.The main processes of the hydrological cycle are variable in time and space.The temporal variation is observed in input temporal data as precipitation, potential evapotranspiration, temperature, etc.The spatial variation is mainly detected in soil properties, soil moisture, land use, etc., which are represented as maps.In some cases, simple aggregation of temporal variability can misrepresent important physical processes.Different hydrological processes as runoff production and channel network routing are affected by this temporal and spatial variability which must be considered in the selected model (Wood et al., 1988;Gan and Biftu, 1996).Heuvelmans et al. (2004) studied the spatial variability of soil parameters using regionalisation schemes and concluded that clustering of parameter sets gives a more accurate result than the single parameter approach and is, therefore, the preferred technique for use in the parameterisation of ungauged sub-catchments as part of the simulation of a large river basin.
The scale effect by spatial and temporal aggregation of nonlinear processes is complex because the mean process response is not the result of the true mean parameter; therefore it is necessary to introduce effective parameters in the model in order to overcome this problem.The effective parameters usually are: physically senseless, non-stationary (they are not valid if used at different ranges obtained during calibration), and highly uncertain.The proposed solution is to use distributed modelling with a small temporal discretization.
The TETIS model is a freeware which has been developed at DIHMA-UPV as a continuous and event model and corresponds to a distributed hydrological conceptual model, which has been chosen because of their very good performance in different basins and climates; Vélez (2001), Vélez et al. (2002aVélez et al. ( , b, 2007) ) and Francés et al. (2007).The TETIS model uses the kinematic wave methodology coupled to the basin geomorphologic characteristics in order to route the flow along the channel network.This procedure is known, according to Vélez (2001), as the "Geomorphologic Kinematic Wave" or GKW.
It is well known in hydrology that the quality in the output data of conceptual rainfall-runoff models depends on: the quality of input data, the model structure and the calibration process (Sorooshian et al., 1993;Lidén andHarlin, 2000 andMadsen, 2000).According to Grayson et al. (1992), conceptual models are better than physically based models since they are faster computationally and have less number of parameters when the suitable scale is adopted.Likewise, the Hydrol.Earth Syst.Sci., 13,[229][230][231][232][233][234][235][236][237][238][239][240][241][242][243][244][245][246]2009 www.hydrol-earth-syst-sci.net/13/229/2009/ conceptual models have a great capacity to compensate for errors and they are difficult to calibrate because they have been created for specific conditions.Sorooshian et al. (1993) indicate that conceptual rainfall-runoff models are difficult to calibrate using automatic methodologies.The successful application of a conceptual rainfall-runoff model depends on how well the model is calibrated.Beven (1989) has presented the limitations of current rainfall-runoff models and argued that the possible way forward must be based on a realistic assessment of predictive uncertainty.For this reason, Beven (1989) has pointed out that many distributed models are just lumped parameter models with a finer mesh, and a closer correspondence between model equations and field processes.This author also indicates that long-term records are required, which not always are available.In addition, the calibration must be performed under a physical point of view in order to interpret properly the parameters and to give them a correct interpretation to the spatial variability.
According to Beven (2000), the calibration of a model takes the same characteristics of an adjustment as a multiple regression, where the ideal parameters will be such that minimize the residual errors, but if there are residual errors it implies that uncertainty exists in the calibrated model.An important limitation of conceptual models, according to Gan and Burges (1990), is that they are not interpolative and there is no guarantee that the model predicts exactly a value during validation when it is used beyond the range of calibration.Likewise, it is recommended that a model should be tried in a range drastically different from that they used during the calibration.
The challenge is to estimate the best parameters set of distributed conceptual models.Due to the inability to accurately measure distributed physical properties of environmental systems, calibration against observed data is typically performed, which is most often achieved with limited rainfall-runoff data.The equifinality noted by Beven (1989) indicates that given the complexity of such models, many different combinations of parameter values may simulate the discharge equally well.These parameter sets may be located throughout many areas of the parameter space (Duan et al., 1992;Beven, 1989).This uncertainty of the appropriate parameter values yields predictive uncertainty, as has been demonstrated through applications of the Generalized Likelihood Uncertainty Estimation methodology (Beven and Binley, 1992;Freer et al., 1996;Beven, 2000).Lidén and Harlin (2000) mentioned that the risk of obtaining a bad model performance, in a conceptual rainfall-runoff model, increases as the initial state moisture condition diminishes.
According to Refsgaard and Knudsen (1996), the main criteria for the evaluation of the performance of hydrological models include the joint visualization of the observed and simulated hydrographs, the general water balance, the efficiency coefficients and the dispersion graphs.
Efficiency indexes must be treated with care in the models, since reliability and validity are not equivalent and the first one does not imply the second one, increasing the size of a sample does not make the estimation more precise with regard to the simulation, it only diminishes the standard deviation.It is also true that a model cannot be estimated reliably based on a small sample, and even if a longer sample is obtained, it does not guarantee a more precise simulation (Duckstein et al., 1985).
The validation process is in charge of demonstrating that physical dominant processes in the basin have been simulated appropriately in a specific site.After validation the model is capable of performing predictions that satisfy the precision criteria previously established (Klemeš, 1988;Refsgaard and Knudsen, 1996;Senarath et al., 2000;Andersen et al., 2001).A fundamental principle in the validation process is that the model must be validated for the same type of applications for which it was originally developed (Klemeš, 1986).In addition, the information used during validation must be different enough from that used during calibration.
The purpose of this paper is to show a Regional Water Resources application using distributed conceptual modelling.The initial parameter maps are estimated using available information and then these maps can be globally corrected without losing their internal variability.This splitparameter structure allows for an extrapolating parameter from gauged sites to ungauged points using advantages of distributed models.
In the next section a brief description of the TETIS model is performed.Then the split-parameter structure and automatic calibration methodology are presented where the initial estimation of spatial parameters is briefly described.The altitude effect on the rainfall is presented and calibration results are exposed jointly with validation and simulation.Finally, some discussion issues and main conclusions obtained during the Water Resources Regional study are exposed.

Model description
The TETIS model is a conceptual distributed model where each grid cell represents a tank model with six tanks connected among them.A conceptual scheme of the vertical movement of the water at each cell can be observed in the Fig. 2. The relationships among different tanks are different for every case, but always simple relationships were used.A more detailed description of the model appears in Vélez (2001), Vélez et al. (2002a, b) and Francés et al. (2007).
The first tank corresponds to the snow package, which can either exist or not.The snow melting process used in the model is the classic degree-day method, which is widely used in the scientific literature (Singh, 1989).For the Basque Country Region case study this tank was not included.The second tank corresponds to the static storage, where the only Runoff production in TETIS model and conceptual scheme of vertical movement at cell scale (Francés et al., 2007).flow exiting is the evapotranspiration.This tank also represents the initial abstractions and pond surfaces.The surface storage is the third tank, where the available water that is not infiltrated can be drained superficially as overland flow.The soil infiltration capacity has been associated with the soil saturated hydraulic conductivity.The fourth tank represents the gravitational storage; the percolation process is modelled according to both soil saturation conditions and the transport capacity, in vertical sense; the remaining water is available to feed the interflow.The fifth tank corresponds to the aquifer, where the vertical flow represents the system groundwater outflow and the horizontal flow is the base flow.Last tank represent the channel at the cell, where each cell is connected to the downstream cell according to the drainage network.Indeed, it is a three-dimensional model.In the Fig. 3 the behavior of the horizontal flow in TETIS model is observed.All cells drain towards the downstream cell until they reach the channel.Once the channel is reached the flow routing is performed according to the GKW methodology.
The time series required during the model execution are discharge, rainfall, evapotranspiration, snow water equivalent and temperature in case that the snow exists.Cartographic information uses raster format maps.Digital Terrain Model (DTM) and soil properties (available water and saturated hydraulic conductivities) are obtained based on soils studies, land use, geological maps, edaphologic information, hydrogeological data and other environmental topics that could be interesting and are available at the study area.(Vélez, 2001).
The infiltration model and the flow channel routing model proposed in TETIS include a few correction factors which correct globally for the different soil properties maps instead of each cell value of the calibration maps, thus reducing drastically the number of factors to be calibrated.This strategy allows for a fast and agile modification in different hydrological processes.These correction factors can be found using an automatic calibration.In general, if the TETIS model has been calibrated adequately these values must not change along the basin, thereby allowing extrapolation to ungauged subbasins.
TETIS model incorporates a correction factor (rainfall interpolation factor) which includes the variation of precipitation with the altitude.The rainfall interpolation factor has been called β and is given in mm/m.In those cases where it is not possible to calibrate the model due to a lack of water in the balance, then the rainfall interpolation factor can be included as an additional variable to be calibrated.
The routing along the channel network was carried out using the GKW, where nine geomorphologic parameters are required, which can be obtained from potential laws.The coefficients and the exponents of these potential relationships can be obtained with a geomorphologic regional study for hydrological homogeneous zones.The Basque Country Region uses the geomorphologic parameters recommended in scientific literature (Leopold et al., 1964;Vélez, 2001), since a geomorphologic regional study for hydrological homogeneous zones was not carried out.
All storages have relevancy inside the initial soil moisture conditions (state variables): snow pack supplied as snow water equivalent in mm.The static storage is given as a percentage of the maximum capacity of the static storage.The initial soil moisture in the surface storage, the gravitational storage and the aquifer are given as a water column in mm.Finally, the initial condition in the channel is supplied to the model as a percentage of the bankfull section.These initial conditions variables are global, e.g. they have the same value for all cells and must be adjusted or calibrated by means of automatic or trial and error procedures.In the case study if the temporal scale was a daily scale, then it is possible to use Hydrol.Earth Syst.Sci., 13,[229][230][231][232][233][234][235][236][237][238][239][240][241][242][243][244][245][246]2009 www.hydrol-earth-syst-sci.net/13/229/2009/ a previous period of time to calibrate in order to generate initial soil moisture conditions as real as possible.Normally, a range between one and two months has been used as previous period of time to calibrate to obtain an initial soil moisture condition, also known as a "warm up" period.
The TETIS model uses the inverse distance method to interpolate spatially temporal series of rainfall, evapotranspiration, temperature and the snow water equivalent.Due to the existence of numerous stations, it is necessary to indicate in the model the number of nearest stations to each cell, to carry out the spatial interpolation.

Split-parameter structure and automatic calibration
The split-parameter structure refers to the way the model estimate the parameters maps.The parameter estimation methodology tries to distinguish between the effective parameter used in the model at cell scale, and the watershed characteristic estimated from the available information, being the best estimation without losing its physical meaning.The relationship between them can be considered as a correction function or, in its simple form, a correction factor (Francés et al., 2007).
Initially, the maps are estimated a priori using environmental and available information (Puricelli, 2003).Then, correction factors are used to modify globally the previously estimated maps.In this way, the spatial variability captured in the initial estimated maps is kept and a global change in magnitude of parameter maps is performed with correction fac-tors (Francés et al., 2007).This approach reduces drastically the number of parameter to be calibrated, because only the common correction factors are calibrated instead of parameter maps (number of parameters times the number of cells).The correction factor can take into account the model input errors, the temporal and spatial scale effects and the watershed characteristics.Therefore, it is reasonable to assume that the correction factor is the same for each parameter to all cells within the watershed.Finally, the calibration of the correction factors is performed if there are available observed episodes.
According to the hydraulic behavior at cell scale, the required maps and correction factors are related in Table 1, where it can be found a brief description of the parameter structure and the information used during the initial estimation of the parameter maps.The required maps are the maximum static storage or available water H u , the vegetation cover index map λ v , the saturated hydraulic conductivity of soil k s , the overland flow velocity v h , and the saturated hydraulic conductivity of subsoil k p .The flow velocity v(t) is calculated according to the GKW procedure and correction factors R 1 , R 2 , . . ., R 9 must be calibrated.
The next sections are dedicated to outline briefly the initial estimation of parameter maps and the calibration strategy.

The initial estimation of spatial parameters
The initial estimation of spatial parameters deals with the prior estimation of three soil properties, according to Puricelli (2003) Fig. 4. Spatial variability between cartographic units (Puricelli, 2003).
conductivities of soil and subsoil, respectively.This prior estimation is made of edaphologic and not strictly edaphologic information, based on available data.This implies the interpretation (or inference) of the land properties (Abbaspour and Moon, 1992), extending the analysis to ground mean values.This work requires the knowledge of factors affecting the selected variables and their impact.The statistical technique used during the map estimation is the multivariate weighted least square method, in which the three hydrological parameters were related to environmental variables, deduced from available information.
The proposed procedure was carried out in four stages: 1) Describing land properties showing their spatial variations.
2) Establishing a sampling criterion to obtain a discrete set of representative values at the work area.3) Defining the general tendency of expected values.4) Finally, obtaining the thematic maps of physical soil properties.
The cartographic representations of environmental variables, land characteristic and soil properties are related to elements of the landscape.Such elements are related to the spatial distribution of hydrological properties.Consequently, it is possible to define two types of variables: a) Main variables: directly related to the hydrologic behavior as soil and subsoil hydraulic conductivity and available water.
b) Environmental variables: those variables capable of explaining part of the spatial variability of main variables.
The objective of the initial estimation is to obtain as reliable as possible the estimation of three main variables (maximum available water and saturated hydraulic conductivities of soil and subsoil, respectively).Consequently, a reasonable way to analyze its spatial variability is based on its relation with environmental variables.In real life, the available environmental variables are geology, land use and edaphology, which use categorical scales and are plotted on a map by means of cartographic units.Thus, for example, in the case of litology, the categories correspond to the different types of rocks described in the geologic map.Each category has associated a set of values corresponding to different physical properties, which are based on continuous scales of values.This association is actually stored into a database, where each cartographic unit is a record and properties are the different fields, including modal values.It is considered that the spatial variability of soil properties is significantly smaller within each cartographic unit than among adjacent units.
In Fig. 4 these concepts are exposed graphically.A small region can be observed in the part (1), whose physical properties are expressed at two cartographic units.According to cross section A-A', the variability of a physical property (k s ) between the cartographic units can be observed in part (2).Thus, it can also be verified that it is possible to define a general tendency of the values throughout this profile according to the expected spatial variability.Finally, if a cartographic unit is isolated, as shown in part (3), then it can be verified that a variability among the unit exists, which oscillates around its modal value.
The distributed hydrological models force one to use different physical properties extracted from the cartographic units.This implies an interpretation and process of the categorical information in order to incorporate its variation.
The most significant environmental variables able to explain the spatial distribution and the variability of the saturated hydraulic conductivity are: topographic index, slope, curvature, cumulative area and geology (McKenzie and Austin, 1993;Boer et al., 1996).
The relationship between the vegetation and the landscape has been extracted from DTM (Bolstad et al., 1998).The vegetation is a variable influenced by different edaphic qualities including the land use and vegetation coverage.
Once we've obtained the thematic layers corresponding to the main and environmental variables, the sampling scheme was based on assigning the modal value of each main variable to a portion of land which was considered as a sampling site.This portion of land is obtained from the intersection between all the cartographic units.
The mean values of the main and environmental variables were assigned to the new cartographic units.In the case of the main variables, the assigned weight was the area at each soil cartographic unit, as shown in Fig. 5.In the case of continuous environmental variables, the weight was the average value of all pixels affected by the intersection.The result was a finite number of distributed values on the study area, which represent the values of the main variables based on the spatial distribution of the environmental variables.
The method to estimate the parameter maps was based on the statistical adjustment among the environmental variables and the main variables.A classical statistical methodology as multivariate weighted least-squares was carried out because it is simple and can be easily adapted at cell scale.Different  (Puricelli, 2003).methodologies such as krigring, incorporating external information, could be applied to obtain the spatial distribution of the modal values among cartographic units, but those methods were not adopted because they are time-consuming.
The scheme of multivariate weighted least-squares proposed by Montgomery and Runger (1996) allows one to combine quantitative variables with categorical data and their interactions.The objective function to minimize is: where y i is the main variable, x ij is the environmental variable j , corresponding to the main variable, ω i is the main variable weight, β 0 is the estimation value when all the environmental variables are null and β i are the environmental factors.
The dummy variables (also known as indicator variables), are used to incorporate qualitative information in leastsquares procedures (Montgomery and Runger, 1996).Thus, a new dummy variable can be defined as the product of the other two dummy variables.
Using this procedure the quantitative information is maximized, categorical or continuous, being able to incorporate qualitative information and maintaining the coherence with the previous physical knowledge of the main variable.
Finally, the results at each cell are processed using GIS in order to obtain the final thematic layers, adding a spatial term to the original cartographic units according to: Since y * i is the estimated main variable value at cell i, y i is the main variable value after the multivariate weighted leastsquares at each cell in the cartographic unit.E (y i ) is the expected value of the main variable for each cartographic unit, and y c i is the modal value of the main variable, assigned to each cartographic unit.
In Table 1 are indicated the model parameter (soil properties at cell scale) and the source map used for their estimation.Interflow velocity uses the same initial estimation than infiltration capacity, which means than R6 must take into account also the anisotropy between the vertical and horizontal conductivities.Similar for the base flow velocity.With respect the groundwater outflow velocity, R7 depends on the basin water balance.A detailed description of each parameter estimation and specially the relationship among the main, environmental and dummys variables can be found in Puricelli (2003).
Finally, three maps of soil parameters are obtained and requested by model TETIS.Is shown in Fig. 6 the parameter map (H u ), obtained with the proposed initial estimation process, where the spatial variability among the Basque Country Region can be observed.
The GKW parameters can be found in Table 1 of Francés et al. (2007).

Calibration procedure
The calibration of most complex hydrological conceptual models has been performed using traditional manual calibration, using the trial and error as suitable methodology.In this type of calibration it is necessary to be an expert user in order to obtain results in a fast and reliable form.The automatic calibration procedures need an initial region of the parameters (physical and feasible values), which can be similar to neighbour basins or hydrologically and climatologically similar basins; Boyle et al. (2000).
The estimation of conceptual model parameters has proven to be a difficult task for several reasons, which include the existence of regions of attraction and multiple local optima.In the last few years numerous technologies have been developed for the automatic calibration.The methods of global search are the most popular, among them the genetic algorithms and the SCE-UA; Madsen (2000).Different authors have performed comparisons among different methods, among them the Simulated Annealing, Genetic Algorithm, MultiStart Simplex, MSX and the SCE-UA, emphasizing that the latter behaves better because the obtained parameters are more grouped and the number of evaluations of the objective function are lower (Duan et al., 1992;Sorooshian et al., 1993;Duan et al., 1994;Yapo et al., 1998 andEckhardt andArnold, 2001).Gan and Biftu (1996) mention that the SCE-UA is ideal to be operated by users who do not know the model and MSX needs an expert user and a process for stages, being more inefficient computationally whereas the SCE-UA obtains results in one execution.
And most recently, the Multiobjective Shuffled Complex Evolution Metropolis algorithm, proposed at the University of Arizona, is best suited for hydrologic model calibration applications that have small parameter sets and small model evaluation times (Tang et al., 2006 andVrugt et al., 2006).
The model performance criteria in different optimization methods take different execution times, which influences the measurement of their behavior (Thyer et al., 1999).Gan and Burges (1990) mention that the mass balance is the most efficient criterion; even tough iterative processes are used.They suggest that there is not an universal and valid scheme to calibrate a conceptual model following manual or automatic procedures.They recommend as optimization algorithms those based on direct search, considering them to be more robust.Numerous criteria can be found in the literature.
The TETIS model includes different objective function criteria, such as (1) the error in volume (BE), obtained according to the expression: where Vp is the total simulated volume and V o is the total observed volume.A positive value indicates an underestimation and negative values an overestimation, where zero is the expected value.
(2) The Root Mean Square Error, RMSE.This efficiency criterion is widely used: where Q i is the observed discharge at time i, Qi is the simulated value at time i and n is the total number of observations.This criterion looks for the minimum value of RMSE.
(2) The efficiency coefficient, E. It is estimated according to the expression: where, the mean value of the observed discharges is Q.This is also known as the efficiency index of Nash and Sutcliffe (1970).This criterion is commonly used because it involves standardization of the residual variance and its expected value does not change with the length of the record or the runoff magnitude (Kothyari and Singh, 1999).A perfect adjustment suggests a value equal to one; if the value is zero it indicates that the model is not better if it is compared by one variable model (for example, an average value) and negative values indicate that the model behaves worse (Beven, 2000).
Recently, Bárdossy (2007) show how parameter sets can be transferred from gauged basins to close ungauged catchments, if performance indexes and annual statistics in donor basins are good.Therefore, after the calibrating process at each basin, the correction factor sets could be used in ungauged catchments, if parameter map have been estimated previously.

Case study: basque country region
The Regional Water Resources study carried out at the Basque Country Region required an intensive information analysis and data process because collected data have different resolution, origin, size, quality and quantity.Additionally, other related topics to the case study as rainfall increment with altitude, calibration process, validation strategy and simulation are included in this section.
The spatial and temporal quality data analysis was performed in order to reduce data uncertainty.This analysis was carried out discarding any doubtful information and correcting temporal data if additional information was available.Finally, a field trip along the Basque Country region allowed improving the spatial pattern of soil properties.

Hydrologic information
The meteorological and hydrological information have been supplied by different entities, namely: the National Meteorological Institute, INM; the Basque Government, GV; the Navarra Government, GN; the Statutory Deputation of Guipuzcoa, DFG; the Statutory Deputation of Biscay, DFV and the Hydrographic Federation of the Ebro River, CHE.
In the Basque Country Region case study a daily series has been used.Analyzing data was not an easy task, mainly because they were supplied from different sources and some data had a different reference hour.This was the case for pluviographic data with a reference time at 12 a.m. and pluviometric data with a reference time at 9 a.m.For this reason it was necessary to process the information in order that the reference time was the same in all stations.Basically, the modifications consist of estimating the daily value as the sum of two thirds of the observed value in the day plus one third of the registered value the next day.
The discharge information was supplied by twelve different entities with 143 gauge stations data but just 113 gauge stations were selected, from which 71 stations were used during the calibration and validation process.It is important to mention some of those stations required to re-establish the natural regime and some flow gauges present a delay with regard to the rain, so it was necessary to realize the respective modification.
The rainfall information consists of a data set of 248 rain gauge stations with a total of 3767 years of daily data.Additionally, 88 climatological stations were used to estimate the PET.In each cell in the case study the information of the three nearest stations was used in order to estimate the spatial distribution of the input variables.

The cartographic information
The original spatial information must be digitized in vectorial format.The information analysis, estimation methods, spatial distribution of hydrological parameters and final results are processed using GIS tools and databases.
The DTM was supplied with a cell size of 25 m×25 m.A gross discretization may lead to poor simulation results whereas a very fine discretization would require more input data and significantly increased computation time and space during automatic calibration with little increase in accuracy.Therefore, the cell size selected was 500 m×500 m.This decision implies that the remaining information must be processed to the same resolution.The DTM was used to extract the basic information: flow direction map to establish connectivity among cells, the cumulative drainage area used to estimate flow velocity using GKW procedure, the slope map used to calculate the runoff velocity and flow velocity using GKW method.
In addition, based on available information (Table 1) the needed by TETIS model four parameter maps have been obtained: the three main soil characteristics (maximum available water and saturated hydraulic conductivities) and the vegetation index cover map as mentioned in previous section.

Rainfall increment with altitude
In those situations in which the precipitation data network does not have rain gauge stations located at high elevations, the main risk is to underestimate the precipitation value at high points.This effect could be possibly the greater source of error in the water balance (Daly et al., 1994).
In the Basque Country Region case study it was initially supposed that the correlation between rainfall and elevation does not exist (β=0.0), which contradicts the observed data on some rain gauge stations.The valley direction can determine the precipitation evolution in some zones of the Basque Country Region.In order to consider the variation of the precipitation with the altitude, it was necessary to estimate this gradient.Part of the supplied data consists of annual rainfall change rates with the elevation above sea level.These rates can be used to estimate a rainfall interpolation factor using the number of rainy days throughout the year.The obtained results were satisfactory once we clarified all the different types of influences on the precipitation.For example, in the Oria River basin the elevation and precipitation analysis were carried out.If all data are analyzed jointly, the influence between precipitation and elevation cannot be detected, as shown in the Fig. 7 (continuous line).Therefore, the basin has been differentiated in two zones according to the prevailing wind direction.The results of the lower basin zone are shown in Fig. 7 (circle points and dashed line), where the topographic influence is quantifiable around 90 mm, each 100 m.The constant of the equation indicates a coastal annual precipitation of approximately 1500 mm which is coherent with reality.In the upper basin sector (filled triangle and dashed dot line) a remarkable correlation is also made visible.It provides similar coefficients to the lower zone; 100 mm, each 100 m.
In the case of the Oria River basin the analysis allowed one to conclude that the precipitation gradient is around 100 mm for each 100 m of elevation increment due mainly to the wind effects.This effect was considered in the model using the rainfall interpolation factor during calibration, which establishes the rate of change of the precipitation with altitude, since the groundwater outflow is the only degree of freedom capable of adjusting the water balance, if we assumed no errors in the data.This is the reason that this study highlights  the necessity to count on rain gauge networks sufficiently extensive to diminish the use of this global factor.

Automatic calibration
The process of automatic calibration was carried out using the SCE-UA methodology.This procedure was performed independently at every basin or hydrological system, considering that nearby basins must have a hydrological similar behavior.Additionally, during the calibration process as a priori division was proposed, basically because the northern basins must have different behavior than southern basins, but must be similar inside each region.
The TETIS model allows one to choose the objective function during the calibration process, in the case study the RMSE has been selected.In addition, automatic calibration procedures require a search range and the initial values which can be observed in the Table 2. Nevertheless, in every particular case some of these values have been diminished or extended, as the hydrologist considers it suitable.This is quite important because the hydrologist experience and field knowledge are included the calibration process at this point.
Concerning to the equifinality problem it was solved using "hydrological sense" in: 1.The search range of each calibration variable (correction factors and β in the automatic calibration process, Francés et al., 2007) 2. A final manual correction of these variables, in most cases with the cost of reducing the calibration efficiency Therefore, the calibration process requires an iterative process, where correction factors are estimated automatically using SCE-UA, and then the analysis of the results is performed, checking the hydrological sense and the water balance in each hydrological system, in order to reduce the uncertainty associated with the parameter maps and model performance.The meaning of "hydrological sense" refers to the participation in the calibration process of local experts, the increase of our basins knowledge by several field trips and analyzing carefully the neighbor basins behavior, the basin water balance and the internal flow distribution into the different flow components.The initial soil moisture conditions are not included during calibration because a warm-up period of one or two months was included.
In Table 3 one can see the correction factors and the rainfall interpolation factor (β) obtained after applying the automatic calibration method SCE-UA.The correction factor variability among basins is not too strong except for correction factors R-3, R-5 and β.However, important differences can be observed between north and south regions, especially for the mentioned variables.Only two basins were calibrated with groundwater outflow; R-7 is greater than zero and therefore β must be zero; both basins located in the south zone.Differences in magnitude among correction factors are explained due to correction factor globally modifies the parameter map and it is in charge of different hydrological process.For example, R-3 and R-6 globally correct the parameter map (k s ), but R-3 is trying to modify the conductivities in order to improve the infiltration process and R-6 is multiplying the conductivities in looking for a good representation of the horizontal flow movement along the cell represented as interflow.In this way the order of magnitude for each correction factor can be explained.
There are two types of correction factor fluctuations among the 20 sets of correction factors shown in Table 3: 1.For correction factor R-7 (groundwater outflow) and for β; (precipitation increment with altitude) the differences from one basin to another are not random, because these two variables depend on the basin water balance of each basin.In fact, they can not be automatically transferred to other basins without a minimum of expert judgment.
2. For the rest of correction factors, the fluctuation is random (it must be) and reflects the correction factor uncertainty produced mainly by the basin differences of each flow gauge station used for calibration: -Different spatial information (in this case study not different source of information) and inputs (location, density,. . . ) and their uncertainty.
-Different hydrology at each region, or equivalently, different model conceptualization error.
The results of the water balance using the calibrated correction factors are shown in Table 4.In the case of basins located in the northern zone a high annual precipitation was observed, where interflow (IF) is the most important component, and overland flow (DR) and base flow (BF) show similar behavior.In southern basins the precipitation (PPT) and actual evapotranspiration (RET) are lower and the relevant Hydrol.Earth Syst.Sci., 13,[229][230][231][232][233][234][235][236][237][238][239][240][241][242][243][244][245][246]2009 www.hydrol-earth-syst-sci.net/13/229/2009/ flow is BF, where DR and IF are not as important.These aspects can be confirmed with land use and cover maps, where forests are located in the north basins and vegetation changes when you move southward.The relief and geomorphologic configuration also contributes to this effect.
In Fig. 8 (part a) one sees that consecutively the monthly index improves with regard to the daily index, which is favourable when a water resources study is performed.In Fig. 8 (part b), for the studied zones, the relation was presented between the area and the monthly Nash and Sutcliffe efficiency index.In conclusion, excellent calibrations were obtained in 10 cases (with E daily over 0.8), acceptable results at 9 basins (E daily between 0.6 and 0.8) and the one case left was moderate (E daily equal to 0.5); the monthly efficiency indexes show excellent results in most cases.In Fig. 9 the results are shown after calibration process at "A3Z1 Altzola" gauge station located in Deba River basin, where the error in volume was −1.99%.The performance efficiency indexes were 0.90 and 0.95 for daily and monthly data, respectively.This calibration results can be considered as excellent.

Validation process
Distributed models must be validated temporally, spatially and spatial-temporarily, according to the available data.The first case is the validation using the same gauge station used during the calibration but with a different period of time.Spatial validation is performed using the same period of time used during the calibration but in another subbasin, usually located upstream.And the spatial-temporal validation is performed using a different period of time and a different gauge station.All these validation approaches were performed in the Basque Country Region; Table 5 contains the  summary results, and the validation process was carried out using available data, therefore the validation period length varies from a few months until several years.In those basins with short data length calibration was not carried out and the validation process was performed using the correction factors of the closest more similar basin (i.e.Bidasoa used correction factors from Oiartzun).Northern basins have shown a better performance than southern basins on a daily scale.But the southern basins showed a better efficiency index at the monthly scale.In general, the volume error was greater at the small basins.It was not possible to include all detailed validation results in this paper but they can be found at INTECSA-INARSA-UPV (2004).Table 5 and Fig. 10b show the basin outlet discharge model efficiency (an inverse measure of uncertainty) at gauge stations not used for calibration and for same and different periods to the calibration ones.This model efficiency is representative what will happen simulating the discharges at truly ungauged points.The model efficiency includes the uncertainty and/or errors in the inputs (precipitation at gauge stations and PET), initial parameters estimation, model conceptualization and parameters calibration.I.e., the objective is to asses the uncertainty of the state variable of interest (which is the interest of the end user), not only to estimate the parameter uncertainty.
The worst efficiencies in these not calibrated points are related in all cases with karstic springs close to the flow gauge stations and/or very small basins with a few number of simulated cells.The first one is clearly a model conceptualization error (TETIS is not intended to simulate or predict concentrated base flows).The second case is clearly the initial parameter estimation uncertainty: when the basin size is smaller than the scale of the spatial information used for this estimation, it is a "lottery" (the cell parameter estimation uncertainty) to "hit" the actual value (always unknown).
In general, the quality of the southern basins information is poorer than for the northern ones (lower spatial resolutions, lower density of input gauges and lower quality of flow gauges).But this is not reflected in the efficiency results, mainly because the other sources of uncertainty have bigger effect on the discharge simulation efficiency.
Figure 10b (and also the Fig. 8b for the calibration) shows a dependence between the variability of the model efficiency and the basin area.It can be shown clearer if we compute the coefficient of variation of monthly E for four basin area ranges and plot them versus the intermediate basin area value.
These two dependences can be explained mainly by the initial spatial parameter estimation uncertainty at cell scale.The effect of this uncertainty is maximum for basin area equal to one cell (0.025 km 2 ) and its effect reduces with basin area because the discharge at the outlet is the sum of all cell runoffs within its basin.And by the Central Limit theorem, sums in statistics always reduce the new variability.
If the uncertainty of the simulated discharges at ungauged points is measured by the monthly value of the Nash-Sutcliffe efficiency index (monthly E) then, the regression equations for the mean monthly E and its coefficient of variation are very accurate indicators of this model uncertainty for our case study.It must be taken into account that it includes all sources of uncertainty and errors (as explained previously) and can be used as an example for other researches and engineering works.bad quality information, a scale problem at small basins due to cell size, or strong karstic presence.
Figure 10 shows the behavior of the monthly efficiency index at the basins where the validation has been performed.Those negative indexes that indicate that validation was not satisfactory were not included in the graph.The results indicate that the monthly indexes improve with regard to the diaries; same behaviour was reported during calibration.According to the monthly efficiency indexes just 11 cases are not satisfactory (E monthly less than 0.5), 10 are acceptable (E monthly between 0.5 and 0.7) and rest of the 41 cases are excellent (E monthly over 0.7).On the daily scale 20 cases are not good, 17 can be considered as good and 25 are excellent.
The northern basins show different characteristics among them, at in the case of the Ibaizabal River basin, where the Karrantza and Herradura Rivers show a very different behav-ior with regard to the nearest Nerbioi and Ibaizabal basins, although the Karrantza River presented poor quality data.
Once the tables were analyzed with the results of the temporal, spatial and spatial-temporal validation it is possible to conclude that validation has been performed satisfactorily and it is feasible to proceed with a worthy simulation.

Long-term simulation
The Basque Country Region includes a very extensive area, which has been divided into hydrographic basins.In this way, the number of points used during calibration was different to the number of points used in validation and different to the points used during the simulation.For example, the total area of the Deba River basin is 531.75 km 2 at the Bay of Biscay.The basin area at the most downstream flowgauge "A3Z1 Altzola" is 469 km 2 , where TETIS model was calibrated using flow data at this site.The final Hydrol.Earth Syst.Sci., 13,[229][230][231][232][233][234][235][236][237][238][239][240][241][242][243][244][245][246]2009 www.hydrol-earth-syst-sci.net/13/229/2009/ correction factors set was validated in the upstream subbasins "A2Z1Aixola", "A1Z1 San Prudentzio", "A1Z2 Zubillaga" and "A1Z3 Urkulu", but the simulation process was also carried out in the total basin area at the Bay of Biscay.The same procedure was carried out at all the Country Basque basins and then simulations were performed at 123 points, including calibration and validation points.
The simulation process used the available information from 1951 or 1961 until 2000 and was calculated from all northern basins to the Bay of Biscay, and all southern basins to Ebro River.At each basin the water balance was performed.In the case of the Deba River basin (531.75 km 2 ), the mean annual rainfall is 1614 mm, of which the evapotranspiration is approximately 761 mm, according to TETIS (with PET of 851 mm).The simulated flow value is 851 mm, without losses.The model distributed the flow in 29.1% as overland flow, 46.8% as interflow and 24.6% as base flow.In Fig. 11 one sees the results of the simulation for the 50 studied years.These results indicate that the Deba River basin never exhausts the flow; the recessions are short in time, indicating that it is a basin with a permanent flow.A similar  analysis was carried out to all basins in the Basque Country Region.In Table 6 are shown some results after the simulation process at the outlet of the main basins in the Basque Country, all of them are ungauged points and the simulation results performed at gauged points are not included in this table.During the simulation process groundwater outflows were considered only at some southern basins where water balance could not be adjusted properly and the northern basin shown more water resources than southern, as was expected.

Discussion and conclusions
In order to solve the equifinality problem it has been proposed the use of the "hydrological sense" in: the search range of each calibration variable (correction factors and β in the automatic calibration process, see Francés et al., 2007) and a final manual correction of these variables, in most cases with the cost of reducing the calibration efficiency.The meaning of "hydrological sense" refers to the participation in the calibration process of local experts, the increase of our basins knowledge by several field trips and analyzing carefully the neighbor basins behavior, the basin water balance and the internal flow distribution into the different flow components.
The main advantages in the use of distributed models are the capacity to represent the space variability of hydrological processes and the capacity to understand the main processes at the hillslope and river basin scale.The TETIS model estimates all hydrological flows that take place at hillslope and river basin scale.In this way, it is possible to explore the origin of the water resources and their implications in water management.The distributed models are spatially robust because they can give results at any point in the basin, most of them ungauged; they do not aggregate information and they can be useful at different time scales, from flash flood to water resources management.Another benefit of distributed modelling is the avoidance of the parameter regionalization (an implicit process in traditional models).This aspect does not result exclusively in an economic saving but avoids sources of additional errors and uncertainties.
The split-parameter structure proposed is fundamental in order to calibrate the distributed conceptual models, basically because the number of parameters needed calibrate is drastically reduced from cell number times parameter maps to just nine correction factors.But previously it required an initial estimation of parameter maps and then the correction factor calibration could be carried out.Therefore, this parameter structure could be adapted to other distributed models.
Distributed models required a relatively expensive initial calibration and validation process.However, the simulation was processed at many points with a comparatively low differential cost.For this reason increasing the output points made an economically and competitive approach.In fact, later studies at the Basque Country Region involved 567 points with a reduced cost allowing one to overcome different studies, such as: risk analysis on small hydraulic projects, environmental flows in regulated rivers, maps of specific Hydrol.Earth Syst.Sci., 13,[229][230][231][232][233][234][235][236][237][238][239][240][241][242][243][244][245][246]2009 www.hydrol-earth-syst-sci.net/13/229/2009/ environmental necessities, underground water resources and environmental necessities, etc.The Basque Country Region River basins were a typical example where it was observed that water quantity estimation is feasible and reliable using the TETIS model.The number of flow gauge stations and the hydrological differences between the northern and southern basins make the case study much more general than usual experimental basins case studies.Also, it was possible to obtain different runoff components and to observe the water balance inside each basin and subbasins.This type of study served to establish the water conditions at natural basins and areas of great concern, where using the water balance of the observed 40 or 50 years could be feasible to forecast the future behavior of the hydrological system, which serves to make better decision about planning, management, conservation and remediation.
The selected model is not a groundwater model and it is not able to simulate karstic effects and explain the bad performance of the model in small areas.The information quality was an important factor to be considered because the basins located in the north showed better and consistent records than the southern basins.The hydrologic behaviour of the north basins is similar among them and different compared to the southern basins, as initially assumed.All differences in flow components among the basins could be finally explained.
The uncertainty in the data makes necessary the use of some groundwater outflow, because in some basins it was not possible to adjust the water balance, which could be due to short-time series, inconsistency and lack of data.The scale effect observed in the small basins during calibration and validation could be lower if the cell size is reduced, but more studies are required.
The rainfall interpolation factor was important in some Basque Country Region basins because it eliminated the groundwater outflow and allowed for closure of the water balance in most basins.Therefore, it was necessary to include it during the simulation process.
In general, the calibration and validation process can be considered as successful because efficiency indexes and volume error were found to be in admissible ranges in most cases; in those cases where it was not satisfactory the reason was identified.The calibration and validation results indicated that the worst calibration set was obtained in small, karstic or poor data basins.The detailed results of this study can be observed in INTECSA-INARSA-UPV (2004).
In general, the use of distributed modeling for the estimation of the water resources in natural basins is feasible and efficient, as shown in the Basque Country Region application.These results suggest that the extrapolating parameter used during the split-parameter structure was successfully carried out and could be used in other basin projects in order to extend and generalize this methodology.

Fig. 1 .
Fig. 1.Division of the study area in 21 systems (different colours) and 41 basins.Thick red line divides the study area into northern and southern basins.Circles and triangles correspond to calibration and validation points, respectively.
Fig. 2.Runoff production in TETIS model and conceptual scheme of vertical movement at cell scale(Francés et al., 2007).

Fig. 6 .
Fig. 6. Results of the initial estimation of parameter map, Hu, in mm.The frequency histogram is in the lower left corner.

Fig. 7 .
Fig. 7. Relationship between elevation and annual precipitation atOria River basin.(The circles belong to the lower zone; and triangles belong to the upper zone, the continuous line represents the adjustment including all points, the dashed line represents the adjustment of the lower zone and dashed-dot line represents the adjustment including the upper zone filled triangle points).

Fig. 8 .Fig. 9 .
Fig. 8. Calibration results, monthly efficiency index versus (a) Daily efficiency index, (b) Area.(North zone is identified by circles and South zone is represented by small squares).

Fig. 10 .
Fig. 10.Validation results, monthly efficiency index versus (a) Daily efficiency index and (b) Area.(North zone is identified by circles and South zone is represented by squares points), (c) Coefficient of variation of monthly efficiency index versus area.

Table 1 .
Parameter structure and information used for the initial estimation of them at cell scale.

Table 2 .
Range and initial value of the correction factors used by the SCE-UA calibration method at Basque Country Region.

Table 3 .
Summary of calibration results at selected gauge stations at Basque Country Region (R: dimensionless correction factors, β is the rainfall interpolation factor in mm/m).

Table 4 .
Summary of water balance after calibration process (PPT: Precipitation, PET: Potential evapotranspiration, RET: Real evapotranspiration, Q obs: Observed discharge, Q sim: Simulated discharge, UL: Groundwater outflow, DR: Overland flow, IF: Interflow, BF: Base flow, BE: percentage of error in volume, E daily: daily efficiency index and E monthly: monthly efficiency index).

Table 5 .
Spatial, temporal and spatial-temporal validation results at the Basque Country Region.(No. years: validation period length, Obs Vol: observed volume, Sim Vol: Simulated volume, BE: Error in volume (%), E daily: daily efficiency index and E monthly: monthly efficiency index).

Table 6 .
Summary of mean results after 40 or 50 years of simulation at ungauged points.PPT: Precipitation in mm/year, PET: Potential evapotranspiration in mm/year, RET: Real evapotranspiration in mm/year, Q sim: Observed discharge in mm/year, V sim: Simulated volume, UL: Groundwater outflow in mm/year, DR: Overland flow, IF: Interflow, BF: Base flow).