Streamflow simulation methods for ungauged and poorly gauged watersheds

Introduction Conclusions References


Introduction
The planning, design, and management of water resources projects require good estimates of streamflow and peak discharge at certain points within a basin.Observed meteorological and streamflow data are initially used for the understanding of the hydrological processes and thus for modelling these processes in order to estimate the streamflow of a watershed.It is likely that most watersheds or basins of the world are ungauged or poorly gauged.There is a whole spectrum of cases which can be collectively embraced under the term "ungauged basins".Some basins are genuinely ungauged, whereas others are poorly gauged or were previously gauged, where measurements discontinued due to instrument failure and/or termination of a measurement programme.Also, the term "ungauged basin" refers to a basin where meteorological data or river flow, or both, are not measured.The international community has recognized this challenging problem and as a result the International Association of Hydrological Sciences (IAHS) had declared the previous decade (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) the "decade of the ungauged basin" (Sivapalan et al., 2003).The IAHS Decade on Prediction in Ungauged Basins (PUB) was a major new effort and an international research initiative to promote the development of science and technology to provide hydrological data where the ground-based observations are needed but missing.This initiative included theoretical hydrology, remote sensing techniques, in situ observations and measurements, and water quantity and quality modelling (Hrachowitz et al., 2013).
In ungauged watersheds, where there are no data, the hydrologist has to develop and use models and techniques which do not require the availability of long time series of meteorological and hydrological measurements.One option is to develop models for gauged watersheds and link the Published by Copernicus Publications on behalf of the European Geosciences Union.model parameters to physiographic characteristics and apply them to ungauged watersheds, whose physiographic characteristics can be determined.Another option is to establish regionally valid relationships in hydrologically similar gauged watersheds and apply them to ungauged watersheds in the region.This approach holds both for hydrograph and flood frequency analysis.The various methods proposed for hydrological prediction in ungauged watersheds can be categorized into statistical methods and hydrological and stochastic modelling methods (Blöschl et al., 2013;Hrachowitz et al., 2013;Parajka et al., 2013;Salinas et al., 2013b).Regionalization techniques are usually applied for statistical methods.These techniques include the regression analyses of flood statistics (statistical moments of flood series) or flood quantiles of gauged watersheds within a homogenous region against geographical and geomorphologic characteristics of the watersheds (Kjeldsen and Rosbjerg, 2002), the combination of single site and regional data, the spatial interpolation of estimated flood statistics at gauged basins using geostatistics (Blöschl et al., 2013), and the region of influence (ROI) approach (Burn, 1990).Then, the established relationships are applied to ungauged watersheds of the region.
In hydrological modelling methods, hydrological models of varying degrees of complexity are used to generate synthetic flows for known precipitation (Singh and Woolhiser, 2002;Singh and Frevert, 2005;Singh, 2012).The complexity of the models can vary from simple event-based models to continuous simulation models, lumped to distributed models, and models that simulate the discharge in sub-daily, daily, or larger time steps.In this approach, a hydrological model is firstly calibrated to gauged watersheds within a region and the model parameters are linked through multiple regression to physiographic and/or climatic characteristics of the watersheds or are spatially interpolated using geostatistics or even using the average model parameter values (e.g.Micovic and Quick, 1999;Post and Jakeman, 1999;Merz and Blöschl, 2004).At the ungauged watersheds of the region, the model with the estimated model parameters is used for hydrological simulation (Wagener et al., 2004;Zhang and Chiew, 2009;He et al., 2011;Wagener and Montanari, 2011;Bao et al., 2012;Razavi and Coulibaly, 2013;Viglione et al., 2013).
The stochastic modelling methods employ a hydrological model which is used to derive the cumulative distribution function of the peak flows.These methods use a stochastic rainfall generation model, which is linked to the hydrological model.The cumulative distribution function of peak flows could be estimated analytically (Iacobellis and Fiorentino, 2000;De Michele and Salvadori, 2002) in the case of a simple hydrological model being used.However, the simplifications and the assumptions made in the analytical derivation of the cumulative distribution function of peak flows may result in poor performance.To overcome this problem the peak flow frequency could be estimated numerically using either an event-based model (Loukas, 2002;Svensson et al., 2013) or a continuous model (Cameron et al., 2000;Engeland and Gottschalk, 2002).
There are difficulties in universally applying the above methods for hydrograph simulation and peak flow estimation of ungauged watersheds.These difficulties arise from the definition of the homogenous regions, the number and the areas of the gauged watersheds, and the different runoff generation processes.The definition, or delineation, of homogeneous hydrologic regions has been a subject of research for many years, and it is necessary for the application of regionalization techniques.The definition of homogeneous regions enables uncorrelated data to be pooled from similar watersheds.A hydrological homogeneous region can be defined by geography, by stream flow characteristics, and by the physical and climatic characteristics of the watersheds.However, problems may arise when an ungauged watershed is to be assigned to a region.The assignment of the watershed to a region is unambiguous when the geographical classification is used and the regions are delineated clearly.On the other hand, the hydrological response of the ungauged watershed may be similar to the response of watersheds belonging in more than one region.This is particularly true for watersheds that are close to region boundaries.In the case of a classification based on stream flow and watershed characteristics, the regions commonly overlap each other.For a classification of regions based on the physical and climatic characteristics of the watersheds, the ungauged watershed could be erroneously assigned to a region.Furthermore, even if a homogenous region is correctly defined and an ungauged watershed is assigned in that region, there should be enough watersheds with extended length of meteorological and streamflow records in order to develop statistically significant regional relationships.However, this is not the case in many parts of the world, where data are very limited, both spatially and temporally.Additionally, the physiographic characteristics, such as slopes, vegetation coverage, soils, etc., and the runoff generation processes (rainfall runoff, snowmelt runoff, glacier runoff, etc.) change as the size of the watershed increases, even in the same region.
The streamflow of a watershed is often measured for a limited period and these streamflow data are inefficient for hydrological model calibration and statistical analysis.In this paper, a technique that couples a hydrological model with artificial neural networks (ANNs) is proposed to improve the streamflow simulation and estimation of peak flows for watersheds with limited streamflow data.In recent years, ANNs have become extremely popular for prediction and forecasting of climatic, hydrologic, and water resource variables (Govindaraju and Rao, 2000;Abrahart et al., 2004) and are extensively reviewed for their effectiveness in the estimation of water quantitative and qualitative variables (Maier and Dandy, 2000;Maier et al., 2010) and in hydrological modelling and forecasting applications (ASCE, 2000;Dawson and Wilby, 2001;Abrahart et al., 2010Abrahart et al., , 2012)).In the context of hydrological modelling, ANNs have mainly been used as rainfall-runoff models for the prediction and forecasting of streamflow in various time steps (Coulibaly et al., 1999;ASCE, 2000;Dawson and Wilby, 2001;Jain et al., 2009;Abrahart et al., 2010).Abrahart et al. (2012) present recent ANN applications and procedures in streamflow modelling and forecasting, which include modular design concepts, ensemble experiments, and hybridization of ANNs with typical hydrological models.Furthermore, ANNs have been used for combining the outputs of different rainfall-runoff models in order to improve the prediction and modelling of streamflow (Shamseldin et al., 1997;Chen and Adams, 2006;Kim et al., 2006;Nilsson et al., 2006;Cerda-Villafana et al., 2008;Liu et al., 2013) and the river flow forecasting (Brath et al., 2002;Shamseldin et al., 2002;Anctil et al., 2004a;Srinivasulu and Jain, 2009;Elshorbagy et al., 2010;Mount et al., 2013).
The objectives of the study are therefore to develop rainfall-runoff modelling procedures for ungauged and poorly gauged watersheds located on different climatic regions.A well-established rainfall-runoff model (Singh, 2012), the University of British Columbia (UBC) watershed model, is selected and applied in five different river basins located in Canada, Cyprus, and Pakistan.Catchments from cold, temperate, continental, and semiarid climate zones are included to demonstrate the procedures developed.In the present study, the term "ungauged" watershed refers to a watershed where river flow is not measured, and the term "poorly gauged" watershed indicates a watershed where continuous streamflow measurements are available for three hydrological years.Two streamflow modelling methods are presented.The first method is proposed for application at ungauged watersheds using a conceptual hydrological model, the UBC watershed model.In this method, most of the parameters of the UBC watershed model take constant values and the precipitation gradients are estimated by analysis of available meteorological data and/or results of earlier regional studies.A second modelling procedure that couples the UBC watershed model with ANNs is employed for the estimation of streamflow of poorly gauged watersheds with limited meteorological data.The coupling procedure of UBC ungauged application with ANNs is an effort to combine the flexibility and capability of ANNs in nonlinear modelling with the physical modelling of the rainfall-runoff process acquired by a hydrological model.

Study basins and database
For the assessment of the developed methodologies, preferably a large number of undisturbed data-intensive catchments located in different climate zones should be studied.However, data for these catchments are very difficult to obtain, which is why the study is limited to five river basins located in different continents.The main selection criteria were accessible hydrometeorological data of good quality and that the studied watersheds represent various climatic types with diverse runoff generation mechanisms.Hence, the developed methodologies are applied to five watersheds located in various geographical regions of the world and with varying physiographic, climatic, and hydrological characteristics, as well as quality and volume of meteorological data.The runoff of all study watersheds contributes to the inflow of local reservoirs.
Two watersheds are forested watersheds located in British Columbia, Canada.The first watershed, the Upper Campbell watershed, is located on the east side of the Vancouver Island Mountains and drains to the north and east into the Strait of Georgia.The 1194 km 2 basin is very rugged, with peaks rising to 2235 m and with mean basin elevation of 950 m (Table 1).The climate of the area is characterized as a maritime climate with wet and mild winters and dry and warm summers.Most of precipitation is generated by cyclonic frontal systems that develop over the North Pacific Ocean and move eastwards.Average annual precipitation is about 2000 mm and 60 % of this amount falls in the form of rainfall.Significant but transient snowpacks are accumulated, especially in the higher elevations.Runoff and the majority of peak flows are generated mainly by rainfall, snowmelt, and winter rainon-snow events (Loukas et al., 2000).The runoff from the Upper Campbell watershed is the inflow to the Upper Campbell Lake and Buttle Lake reservoirs.Daily maximum and minimum temperatures were available at two meteorological stations, one at 370 m and the other at 1470 m, and daily precipitation at the lower-elevation station.In total, seven years of daily meteorological and streamflow data (October 1983-September 1990) were available from the Upper Campbell watershed.
The second study watershed is the Illecillewaet watershed, which is located on the west slopes of the Selkirk Mountains in southeastern British Columbia, 500 km inland from the Coast Mountains.The size of the watershed is 1150 km 2 and its elevation ranges from 400 to 2480 m (Table 1).The Illecillewaet River is a tributary of the Columbia River and contributes to the Arrow Lakes reservoir.The climate of the area is continental, with cold winters and warm summers with frequent hot days, and is influenced by the maritime Pacific Ocean air masses and by weather systems moving eastwards.Average annual precipitation ranges from 950 mm at the mouth of the watershed to 2160 mm at the higher elevations.Substantial snowpacks develop during winter at all elevations of the watershed.The snowpack at the valley bottom is usually depleted by the end of April, but permanent snowpacks and a glacier with an area of 76 km 2 exist at the highest elevations.Streamflow is mainly generated during spring, by rain and snowmelt, and summers, by snowmelt and the contribution of glacier melt (Loukas et al., 2000).Goodquality daily precipitation and maximum and minimum temperature data are measured at three meteorological stations at 443, 1323, and 1875 m elevation, respectively.Twenty years of meteorological and streamflow data (October 1970- September 1990) were used to assess the simulated runoff from the watershed.
The third study basin is the Yermasoyia watershed, which is located on the southern side of mountain Troodos of Cyprus, roughly 5 km north of the city of Limassol.The watershed area is 157 km 2 and its elevation ranges from 70 m up to 1400 m (Table 1).Most of the area is covered by typical Mediterranean-type forest and sparse vegetation.A reservoir with storage capacity of 13.6 million m 3 was constructed downstream of the mouth of the watershed in 1969 for irrigation and municipal water supply purposes (Hrissanthou, 2006).The climate of the area is of Mediterranean maritime climate, with mild winters and hot and dry summers.Precipitation is usually generated by frontal weather systems moving eastwards.Average basin-wide annual precipitation is 640 mm, ranging from 450 mm at the low elevations up to 850 mm at the upper parts of the watershed.Mean annual runoff of the Yermasoyia River is about 150 mm, and 65 % of it is generated by rainfall during winter months.The river is usually dry during summer months.The peak flows are observed in winter months and produced by rainfall events.Good-quality daily precipitation from three meteorological stations located at 70, 100, and 995 m elevation were used.Data of maximum and minimum temperature measured at the low-elevation station (70 m) were used in this study.In total, 11 years of meteorological and streamflow data (October 1986-September 1997) were available for the Yermasoyia watershed.
The fourth and fifth study watersheds, the Astor and the Hunza watersheds, are located within the upper Indus River basin in northern Pakistan.The Astor watershed spans elevations from 2130 to 7250 m and covers an area of 3955 km 2 , only 5 % of which is covered with forest and 10 % covered with glaciers (Table 1).Precipitation is usually generated by westerly depressions, but occasionally monsoon storms produce heavy precipitation.Average basin annual precipitation is about 700 mm and more than 90 % of this amount is snow (Ahmad et al., 2012).Runoff and the peak streamflows are mainly generated by snowmelt and glacier melt (Loukas et al., 2002;Archer, 2003).Mean annual streamflow is about 120 m 3 s −1 , which amounts to 5 % of the inflow to the downstream Tarbela reservoir.Daily precipitation and maximum and minimum temperature data are measured at one meteorological station located at an elevation of 2630 m.In total, nine years of meteorological and streamflow data (October 1979-September 1988) were available from the Astor watershed.The Hunza watershed lies within the Karakoram Mountain Range.The Hunza River flows southwest from its headwaters near the China-Pakistan border and through the Karakoram to join the Gilgit River near the town of Gilgit.The Hunza watershed has a total drainage area of 13 100 km 2 (Table 1) and the entire area is a maze of towering peaks, massive glaciers, and steepsided gorges.The highest mountain peaks within the Hunza Basin area are Batura (7785 m), Rakaposhi (7788 m) and Disteghil Sar (7885 m).The elevation of the Hunza Basin ranges from 1460 to 7885 m.Twenty-three percent of the watershed area is covered by glaciers, including the large Baltoro and Hispar glaciers (Bocchiola et al., 2011;Ahmad et al., 2012).The Hunza Basin is arid and annually receives less than 150 mm of precipitation, mainly in the form of snow, from westerly weather systems.More than 90 % of the annual runoff and peak streamflows are generated by glacier melt (Loukas et al., 2002;Archer, 2003).Mean annual streamflow is about 360 m 3 s −1 , which amounts to more than 13 % of the inflow to the downstream Tarbela reservoir.Daily precipitation data measured at two meteorological stations located at 1460 and 2405 m elevation were used.Data of maximum and minimum temperature measured at the low-elevation station (1460 m) were used in this study.Four years of meteorological and streamflow data (October 1981-September 1985) were available from the Hunza Basin.

Method of analysis
Two methodologies are proposed in this paper for the simulation of daily streamflow of the five study watersheds.The first methodology uses the UBC watershed model with estimated universal model parameters and estimates of precipitation distribution, and it is proposed for use in ungauged watersheds.The second methodology proposes the coupling of UBC watershed model with ANNs, and is intended for use in watersheds where limited streamflow data are available.The UBC watershed model and the two methodologies are presented in the next paragraphs.

The UBC watershed model
The UBC watershed model was first presented 35 years ago (Quick and Pipes, 1977), and has been updated continuously to its present form.The UBC is a continuous conceptual hydrologic model which calculates daily or hourly streamflow using precipitation and maximum and minimum temperature data as input data.The model was primarily designed for the simulation of streamflow from mountainous watersheds, where the runoff from snowmelt and glacier melt may be important, apart from the rainfall runoff.However, the UBC watershed model has been applied to variety climatic regions, ranging from coastal to inland mountain regions of British Columbia, including the Rocky Mountains, and the subarctic region of Canada (Hudson and Quick, 1997;Quick et al., 1998;Micovic and Quick, 1999;Loukas et al., 2000;Druce, 2001;Morrison et al., 2002;Whitfield et al., 2002;Merritt et al., 2006;Assaf, 2007).The model has also been applied to the Himalayas and Karakoram Mountain Ranges in India and Pakistan, the Southern Alps in New Zealand, and the Snowy Mountains in Australia (Singh and Kumar, 1997;Singh and Singh, 2001;Quick, 2012;Naeem et al., 2013).This ensures that the model is capable of simulating runoff under a large variety of conditions.
The model conceptualizes the watersheds as a number of elevation zones, since the meteorological and hydrological processes are functions of elevation in mountainous watersheds.In this sense, the orographic gradients of precipitation and temperature are major determinants of the hydrologic behaviour in mountainous watersheds.These gradients are assumed to behave similarly for each storm event.Furthermore, the physiographic parameters of a watershed, such as impermeable area, forested areas, vegetation density, open areas, aspect, and glaciated areas, are described for each elevation zone and can be estimated from analogue and digital maps and/or remotely sensed data.Hence, it is assumed that the elevation zones are homogeneous with respect to the above physiographic parameters.In a recent study, the UBC watershed model was integrated into a geographical information system that automatically identifies and estimates the physiographic parameters of each elevation zone of a watershed from digital maps and remotely sensed data (Fotakis et al., 2014).A certain watershed can be divided in up to 12 homogeneous elevation zones.The UBC watershed model provides information on snow-covered area, snowpack water equivalent, potential and actual evapotranspiration, soil moisture interception losses, groundwater storage, and surface and subsurface runoff for each elevation zone separately and for the whole watershed.Figure 1 presents the flow diagram of the UBC watershed model.
The model is made up of several sub-routines: the subroutine for the distribution of the meteorological data, the soil moisture accounting sub-routine, and the flow-routing sub-routine.The meteorological distribution sub-routine distinguishes between total precipitation in the form of snow and rain using the temperature data.If the mean temperature is below 0 or above 2 • C, then all precipitation is in the form of snow or rain, respectively.When the mean temperature is between 0 and 2 • C, then the percentage of total precipitation which is rain is estimated by and the remaining percentage of precipitation is snow.Snow is stored until it melts, whereas rain is immediately processed by the soil moisture routine accounting to a sub-routine.Each meteorological station has two representation factors, one for A. Loukas and L. Vasiliades: Streamflow simulation methods snow, P0SREP, and one for rain, P0RREP.These factors are introduced because precipitation data from a meteorological station are point data and they may not be representative of a larger area or zone.If the data are representative, then these parameters are equal to zero.The point station data of precipitation are distributed over the watershed using the equation where PR i,j,l is the precipitation from meteorological station i for day j and elevation zone l, P0GRAD is the percentage precipitation gradient, and elev is the elevation difference between the meteorological station and the elevation zone.
The UBC model then adjusts the precipitation gradient according to the temperature, where ST(T ) is a parameter, which is affected by the stability of the air mass.It can be shown (Quick et al., 1995) that the ST(T ) parameter is related to the square of the ratio of the saturated and dry adiabatic lapse rates, L S and L D , re- versus temperature reveals an almost linear variation between −30 and +20 • C. The gradient of this linear approximation is 0.01; thus ST(T ) can be estimated as where T mean is the mean daily temperature.
The UBC watershed model has the capability of using three different precipitation gradients in a single watershed, namely P0GRADL, P0GRADM, and P0GRADU.The lowelevation gradient, P0GRADL, applies to elevations lower than the elevation E0LMID, whereas the upper-elevation gradient, P0GRADU, applies above the elevation E0LHI and the middle-elevation gradient, P0GRADM, applies to elevations between E0LMID and E0LHI.
The temperature in the UBC watershed model is distributed over the elevation range of a watershed according to the temperature lapse rates.Two temperature lapse rates are specified in the UBC watershed model, one for the maximum temperature and one for the minimum temperature.Furthermore, the model recognizes two conditions, namely the rainy condition and the clear-sky and dry-weather condition.Under the rainy condition, the lapse rate tends to be the saturated adiabatic rate.Under dry-weather conditions and during the warm part of the day, the lapse rate tends to be the dry adiabatic rate, whereas the lapse rate tends to be quite low, and occasionally zero lapse rates may occur during dry weather and night.The lapse rate is calculated for each day using the daily temperature range (temperature diurnal range) as an index.A simplified energy budget approach, which is based on limited data of maximum and minimum temperature and can account for forested and open areas, as well as aspect and latitude, is used for the estimation of the snowmelt and glacier melt (Quick et al., 1995).
The soil moisture accounting sub-routine represents the nonlinear behaviour of a watershed.All the nonlinearity of the watershed behaviour is concentrated into the soil moisture accounting sub-routine, which allocates the water from rainfall, snowmelt, and glacier melt into four runoff components, namely the fast or surface runoff, the medium or interflow runoff, the slow or upper zone groundwater runoff, and the very slow or deep zone groundwater runoff.The impermeable area, which represents the rock outcrops, the water surfaces, and the variable source saturated areas adjacent to stream channels, divides the water that reaches the soil surface after interception and sublimation into fast surface runoff and infiltrated water.The total impermeable area at each time step varies with soil moisture, mainly due to the expansion or shrinkage of the variable source riparian areas.The percentage of the impermeable areas of each elevation zone varies according the Eq. ( 5): where C0IMPA is the maximum percentage of impermeable areas when the soil is fully saturated, S0SOIL is the soil moisture deficit in the elevation zone, and P0AGEN is a parameter which shows the sensitivity of the impermeable areas to changes in soil moisture.
The water infiltrated into the soil must first satisfy the soil moisture deficit and the evapotranspiration and then continues to infiltrate into the groundwater or runs off as interflow.This process is controlled by the "groundwater percolation" parameter (P0PERC).The groundwater is further divided into an upper and deep groundwater zones by the "deep zone share" parameter (P0DZSH).This water allocation by the soil moisture accounting sub-routine is applied to all watershed elevation zones.Each runoff component is then routed to the watershed outlet, which is achieved in the flow-routing sub-routine.However, a different mechanism is employed in the case of high-intensity rainfall events, which can produce flash flood runoff.The runoff from these events is controlled by the soil infiltration rate.For these high-intensity rainfall events, some of the rainfall infiltrates into the soil and is subject to the normal soil moisture budgeting procedure previously presented.The remaining amount of rainfall which is not infiltrated into the soil is considered to contribute to the fast runoff component, which is called FLASHSHARE and is estimated with where FMR is the percentage of the flash share with range from 0 to 1 and is estimated with The UBC watershed model has more than 90 parameters.However, application of the model to various climatic regions and experience have shown that only the values of 17 general parameters and 2 precipitation representation factors (e.g.P0SREP and P0RREP) for each meteorological station have to be optimized and adjusted during calibration, and the majority of the parameters take standard constant values.These varying model parameters can be separated into three groups: the precipitation distribution parameters (namely P0SREP(i), P0RREP(i), P0GRADL, P0GRADM, P0GRADU, E0LMID, and E0LHI), the water allocation parameters (namely P0AGEN, P0PERC, P0DZSH, V0FLAX, and V0FLAS), and the flow-routing parameters (namely P0FSTK, P0FRTK, P0ISTK, P0IRTK, P0UGTK, P0DZTK, and P0GLTK).These parameters are optimized through a two-stage procedure.However, in this paper, the water allocation parameters and the flow-routing parameters are given constant universal values, whereas the precipitation distribution parameters are estimated from the meteorological data and/or using the results of earlier regional studies on precipitation distribution with elevation, as will be presented below.The total number of model parameters for the Upper Campbell and Astor watersheds is 19, for Illecillewaet and Yermasoyia 23, and for Hunza 21, as will be shown below.

Methodology for ungauged watersheds
The five study watersheds were initially treated as ungauged watersheds, assuming that streamflow measurements were not available.However, meteorological data were used from the available stations at each study watershed.The UBC watershed model was used to simulate the streamflow from the five study watersheds.Twelve out of the 17 general varying model parameters were assigned constant universal values, which were either estimated or taken as default (Tables 2 and   3).This work uses the results of a recent paper (Micovic and Quick, 1999) that applied the UBC watershed model in 12 heterogeneous watersheds in British Columbia, Canada, with different sizes of drainage area, climate, topography, soil types, vegetation coverage, geology, and hydrologic regime.Micovic and Quick (1999) found that averaged constant values could be assigned to most of the model parameters.Table 2 shows the averaged values of the model parameters that mainly affect the time distribution of the runoff.
Additionally, the UBC watershed model water allocation parameters P0AGEN, V0FLAX, and V0FLAS were assigned the default values suggested in the manual of the model (Quick et al., 1995).The flow-routing parameter of glacier runoff, P0GLTK, was assigned the value of the rainfall fast flow-routing parameter, P0FRTK, assuming that the response of the glacier runoff is similar to the response of the fast component of the runoff generated by rainfall.The values of these parameters are presented in Table 3. Apart from these parameters, the precipitation distribution parameters were estimated separately from the available meteorological data for each watershed.This estimation procedure is described in the next paragraphs for each one of the five study watersheds.

Estimation of model precipitation distribution parameters for the Upper Campbell watershed
Only one precipitation station was available in the Upper Campbell watershed.For this station the precipitation representation parameters for rainfall and snowfall, P0RREP and P0SREP, respectively, were set to zero.The results of earlier studies on the precipitation distribution with elevation in the coastal region of British Columbia (Loukas and Quick, 1994;Loukas and Quick, 1995) were used for assigning values of precipitation distribution model parameters.In these earlier studies, it was found that the precipitation increases 1.5 times from the coast up to an elevation equal to about two-thirds of the elevation of the mountain peak, and then levels off at the higher elevations.Using this information, the low precipitation gradient, P0GRADL, was estimated from Eq. ( 2), substituting the mean annual precipitation of the lower meteorological station located at 370 m for PR i,j,l , PR i,j,l+1 the increased 1.5 times the mean annual precipitation of the lower meteorological station, and elev the elevation difference between the elevation of the maximum precipitation (two-thirds of the maximum mountain peak, 1490 m) and the elevation of the lower meteorological station (370 m) which equals 1120 m.Hence, the estimated value of P0GRADL was equal to 3.7 %.The elevation where the maximum precipitation occurs (1490 m) defines the value of model parameter E0LMID.The middle and upper precipitation gradients (i.e.P0GRADM and P0GRADU) were set to zero.In this case, it was not necessary to define the model parameter E0LHI, because the precipitation was assumed constant above the E0LMID elevation (1490 m).

Estimation of model precipitation distribution parameters for the Illecillewaet watershed
Three precipitation stations were available at the Illecillewaet watershed located at elevations of 443, 1323, and 1875 m, respectively.The model precipitation representation parameters for rainfall and snowfall and for all three stations were set to zero (i.e.P0RREP(1 . The low precipitation gradient, P0GRADL, was estimated from Eq. ( 2) using the mean annual precipitation at the low-and middleelevation stations and the elevation difference between the two stations ( elev=1323-443 = 880 m).P0GRADL was found to equal 6 %.Similarly, the middle precipitation gradient, P0GRADM, is estimated to equal 5.5 % considering the mean annual precipitation of the middle-and upper-elevation station.The upper precipitation gradient, P0GRADU, was set to zero.The parameter E0LMID was set equal to the elevation of the middle-elevation station, which is 1323 m.The parameter E0LHI was set equal to the highest elevation of the watershed, 2480 m.

Estimation of model precipitation distribution parameters for the Yermasoyia watershed
Precipitation data from three meteorological stations located at 70, 100, and 995 m elevation were available at the Yermasoyia watershed.The precipitation representation parameters for snowfall and for all three stations were set equal to zero, because snowfall is rarely observed (i.e.P0SREP(1) The annual precipitation data of the three stations were compared with the annual precipitation of other stations in the greater area of the watershed.This comparison showed that the three meteorological stations record 30 % more annual rainfall than other stations located at similar elevations.For this reason the rainfall representation parameters for all three stations were set equal to −30 % (i.e.P0RREP(1) = P0RREP(2) = P0RREP(3) = −30 %).The low precipitation gradient, P0GRADL, was estimated using Eq. ( 2) as well as the mean annual precipitation of the lower-elevation station and the mean annual precipitation at the upper-elevation station.The precipitation gradient between the two lowerelevation stations is essentially zero because of the small elevation difference.The lower precipitation gradient parameter, P0GRADL, was estimated to equal 4.9 %.The parameter E0LMID was set equal to the elevation of the upper-elevation station, which is 995 m.The middle and the upper precipitation gradients, P0GRADM and P0GRADU, respectively, were set equal to zero.This means that the simulation was performed with one precipitation gradient.In this case, it was not necessary to define the model parameter E0LHI.

Estimation of model precipitation distribution parameters for the Astor watershed
In the Astor watershed, only the precipitation data of one meteorological station located at 2630 m were available.For this reason and because it was not any information on the distribution of precipitation with elevation, all the model precipitation representation and distribution parameters, i.e.P0RREP, P0SREP, P0GRADL, P0GRADM, and P0GRADU, were set equal to zero.In this case, it was not necessary to define the model parameters E0LMID and E0LHI, which were set equal to zero.

Estimation of model precipitation distribution parameters for the Hunza watershed
Daily precipitation data from two meteorological stations located at 1460 and 2405 m elevation were available at the Hunza Basin.The mean annual precipitation at the two stations was estimated, and it indicated that the precipitation gradient between the two stations was essentially zero.For this reason, and because there was no information on the distribution of precipitation with elevation, all the model precipitation representation and distribution parameters were set equal to zero (i.e.P0RREP(1) = P0SREP(1) = P0RREP(2)

Methodology for poorly gauged watersheds
The streamflow is frequently measured for a limited period of time.These streamflow data are inadequate for peak flow analysis and validation of the simulated streamflow.Unfortunately, there are no specific guidelines about the precise Nat.Hazards Earth Syst.Sci., 14, 1641-1661, 2014 www.nat-hazards-earth-syst-sci.net/14/1641/2014/calibration length of streamflow data needed for optimal hydrological model performance in poorly gauged watersheds (Seibert and Beven, 2009).Several studies in gauged watersheds have shown that, for an acceptable rainfall-runoff model calibration, a large calibration record including wet and dry years (at least eight years) is needed for complex hydrologic models, and the minimum requirements are one hydrological year (Sorooshian et al., 1983;Yapo et al., 1996;Duan et al., 2003).For example, Yapo et al. (1996) stated that for a reliable and acceptable model performance, a calibration period with at least eight years of data should be used for NWSRFS-SMA hydrologic model with 13 free parameters.Harlin (1991) suggested that from two to six years of streamflow data are needed for optimal calibration of the HBV model with 12 free parameters.Xia et al. (2004) suggest that at least three years of streamflow data are required for successful application of their model (with seven parameters) for a case study in Russia.In this regard, few studies investigate the use of limited number of observations for calibration periods shorter than one year.Brath et al. (2004) suggest for flood peak modelling using a continuous distributed rainfallrunoff model that three months are the minimum requirement for flood peak estimation.However, their best results are acquired with the use of one year continuous runoff data.Perrin et al. (2007) found that calibration of a simple runoff model (the GR4J model with four free parameters) is possible using about 100-350 observation days spread randomly over a longer time period including dry and wet conditions.These results were also verified by Seibert and Beven (2009), who showed that a few runoff measurements (larger that 64 values) can contain much of the information content of continuous streamflow time series.The problem of limited streamflow data might be tackled if the data are selected in an intelligent way (e.g.Duan et al., 2003;Wagener et al., 2003;Juston et al., 2009) or by using information from other variables such as data from groundwater and snow measurements in a multiobjective context (e.g.Efstratiadis and Koutsoyiannis, 2010;Konz and Seibert, 2010;Schaefli and Huss, 2011).
The above studies give an indication of the potential value of limited observation data for constraining model prediction uncertainties even for ungauged basins.However, these studies indicated that the results diverge significantly between the watersheds, depending on the days chosen for taking the measurements, and misleading results could be obtained with the use of few streamflow data (Seibert and Beven, 2009).Furthermore, the conceptual hydrological models employed are simple and have a small number of free parameters, and more research is needed for complicated hydrological structures with more than 10 parameters such as the UBC watershed model.In a recent study, the impact of calibration length in streamflow forecasting using an ANN and a conceptual hydrologic model, GR4J, was assessed (Anctil et al., 2004b).The results showed that the hydrological model is more capable than ANNs for 1-day-ahead flow forecasting using calibration periods less than one hydrological year due to its internal structure, and similar results are obtained for calibration periods from one to five years.However, the ANN model outperformed the GR4J model for calibration periods larger than five years as a result of its flexibility (Anctil et al., 2004b).Based on the above studies and discussion, it is difficult to define the minimum requirements for model (conceptual or black-box) calibration for poorly gauged watersheds.Furthermore, model accuracy may also depend on the climatic zone, an aspect that is rarely explicitly analysed.Therefore, we developed a methodology that can make use of limited streamflow information with the internal memory of a non-calibrated semi-distributed rainfall-runoff model and the predictive capabilities of ANNs for poorly gauged watersheds as defined in this study.

UBC coupling with ANNs
The coupling of the UBC watershed model with ANNs is described in this section.ANNs distribute computations to processing units called neurons or nodes, which are grouped in layers and densely interconnected.Three different layer types can be distinguished: an input layer, connecting the input information to the network and not carrying any computation; one or more hidden layer, acting as intermediate computational layers; and an output layer, producing the final output.In each computational node or neuron, each one of the entering values (x i ) is multiplied by a connection weight, (w j i ).Such products are then all summed with a neuronspecific parameter, called bias (b j 0 ), used to scale the sum of products (s j ) into a useful range: A nonlinear activation function (sometimes also called a transfer function) to the above sum is applied to each computational node producing the node output.Weights and biases are determined by means of a nonlinear optimization procedure known as training that aims at minimizing an error function expressing the agreement between observations and ANN outputs.The mean squared error is usually employed as the learning function.A set of observed input and output (target) data pairs, the training data set, is processed repeatedly, changing the parameters of ANN until they converge to values such that each input vector produces outputs as close as possible to the observed output data vector.In this study, the following neural network characteristics were chosen for all ANN applications: 1. Structure of ANNs: feedforward ANNs were used, which means that information passes only in one direction, from the input layer through the hidden layers up to the output layer, allowing only feedforward connections to adjacent layers.

2.
Training algorithm: back-propagation algorithm (Rumelhart et al., 1986) was employed for ANNs training.In this training algorithm, each input pattern of the training data set is passed through the network from the input layer to the output layer.The network output is compared with the desired target output and the error according to the error function, E, is computed.This error is propagated backward through the network to each node, and correspondingly the connection weights are adjusted based on the following equation: where w j i (n) and w j i (n − 1) are the weight increments between the node j and i during the nth and (n − 1)th pass or epoch.A similar equation is employed for correction of bias values.In Eq. ( 9) the parameters ε and α are referred to as learning rate and momentum, respectively.The learning rate is used to increase the chance of avoiding the training process being trapped in a local minimum instead of global minima, and the momentum factor can speed up the training in very flat regions of the error surface and help prevent oscillations in the weights.
3. Activation function.Here, the sigmoid function is used: The sigmoid function is bounded between 0 and 1, and is a monotonic and nondecreasing function that provides a graded, nonlinear response.
The UBC watershed model, as has been previously discussed, distributes the rainfall and snowmelt runoff into four components, i.e. rainfall fastflow, snowmelt fastflow, rainfall interflow, snowmelt interflow, upper zone groundwater, deep zone groundwater, and glacial melt runoff.These runoff components due to errors in measurements and inefficiently defined model parameters may not be accurately distributed, affecting the overall performance of the hydrologic simulation.The UBC watershed model used the parameters with values described in the previous subsection of the paper.In order to take advantage of the limited streamflow data and achieve a better simulation of the observed discharge, the runoff components of the UBC watershed model are introduced as input neurons into ANNs.During the training period of ANNs, the simulated total discharge of the watershed is compared with the observed discharge to identify the simulation error.
The geometry or architecture of ANNs, which determines the number of connection weights and how these are arranged, depends on the number of hidden layers and the number of hidden nodes in these layers.In the developed ANNs, one hidden layer was used to keep the ANNs architecture simple (three-layer ANNs), and the number of the hidden nodes was optimized by trial and error.In this sense, the input layer of ANNs consists of four to seven input neurons, depending on the runoff generation mechanisms of the basin; one hidden layer with varying number of neurons; and one output layer with one neuron, which is the total discharge of the watershed (Fig. 2).Since the various input data sets span different ranges, and to ensure that all data sets or variables receive equal attention during training, the input data sets were scaled or standardized in the range of 0-1.In addition, the output variables were standardized in such a way as to be commensurate with the limits of the activation function used in the output layer.In this study, the sigmoid function (Eq.10) was used as the activation or transfer function, and the output data sets (watershed streamflow) were scaled in the range 0.1-0.9.The advantage of using this scaling range is that extremely high and low flow events occurring outside the range of the training data may be accommodated (Dawson and Wilby, 2001).
However, the final network architecture and geometry were tested to avoid overfitting and ensure generalization as suggested by Maier and Dandy (1998).For example, the total number of weights was always kept less than the number of the training samples, and only the connections that had statistically significant weights were kept in the ANNs. is the number of inputs for the node.The learning rate (ε in Eq. 9) was set fixed to a value of 0.005, whereas the momentum (α in Eq. 9) was set equal to 0.8 as suggested by Dai and Macbeth (1997).

Evaluation of the method
For the four study watersheds, namely the Upper Campbell, Illecillewaet, Yermasoyia, and Astor watersheds, the first three years of streamflow record were assumed to be available for training of ANNs.In this sense, the observed streamflow used as the target output of ANNs was the daily measured streamflow for the hydrological years 1983-1984 and 1985-1986 for the Upper Campbell watershed, the streamflow data for the hydrological years 1970-1971 and 1972-1973 were considered for the Illecillewaet watershed, the data for the hydrological years 1986-1987 and 1988-1989 were used for the Yermasoyia watershed and the streamflow data for the hydrological years 1979-1980 and 1981-1982 were used for the Astor watershed.For the fifth catchment, the Hunza watershed, streamflow data for two hydrological years (1981-1982 and 1982-1983) were used for ANN training.The daily streamflow measurements for the remaining years of record were used for the validation of the methodology in each study watershed.The modelling procedure with this configuration is termed UBCANN, or method with limited data.It should be noted that the early stopping technique was applied to UBCANN to prevent overfitting and to improve the generalization ability of the developed UBCANNs.The last year in each watershed of the training period was used as an indication of the error when ANN training should stop (test set).
For comparison purposes, the UBCANN method was compared with the ungauged application of the UBC model, termed UBCREG, and with the classical calibration of the UBC model in poorly gauged watersheds using the same calibration period for each watershed as defined previously.The latter method is termed UBCCLA and is used for evaluation of the proposed coupling method, UBCANN, for poorly gauged watersheds.The UBC free parameters are optimized through a two stage procedure.In the first stage, a sensitivity analysis of each parameter is performed to estimate the range of parameter values for which the simulation results are the most sensitive.In the second stage, a Monte Carlo simulation is performed for each parameter of each group by keeping all other parameters constant.The parameter values are sampled from the respective parameter range defined in the first stage of the procedure (sensitivity analysis).The parameter value that maximizes the objective function is put in the parameter file, and the procedure is repeated for the next parameter of the group and then for the parameters of the next group.The procedure starts with the optimization of the precipitation distribution parameters and ends with the optimization the flow-routing parameters.The objective function of the above calibration procedure is defined as where V sim and V obs are the simulated and the observed flow volumes, respectively, and NSE is the Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970), defined as where Q obs i is the observed flow on day i, Q sim i is the simulated flow on day i, Q obs is the average observed flow, and n is the number of days for the simulation period.The evaluation of all the applied methods is based on the combination of graphical results, statistical evaluation metrics, and normalized goodness-of-fit statistics.Furthermore, a comprehensive procedure proposed by Ritter and Muñoz-Carpena (2013) for evaluating model performance is tested for all applied methods.Approximated probability distributions for NSE and root-mean-square error (RMSE) are derived with bootstrapping followed by bias correction and enhanced calculation of confidence intervals.Statistical hypothesis testing of the indicators is done using threshold values to compare model performance.More details on the evaluation protocol can be found in Ritter and Muñoz-Carpena (2013).Finally, the streamflow simulation results of the applied methods for ungauged and poorly gauged watersheds were used for frequency analysis of the annual maximum peak flows.This analysis was performed only for the watersheds which have streamflow data for at least six consecutive years.Based on these criteria, the Hunza watershed is excluded for this comparison.The estimated peak flows were compared with the observed peak flows of the four study watersheds (Upper Campbell, Illecillewaet, Yermasoyia, and Astor).Furthermore, the results of frequency analysis of the estimated peak flow from the two methodologies were compared to the results of frequency analysis of the observed peak flows.The frequency analysis was performed using the extreme value type I theoretical distribution (EVI) due to the small sample of the streamflow observations, and due to its simple two-parameter estimation procedure.This distribution is a special case of the generalized extreme value (GEV) distribution, and the GEV distribution is considered in a recent study as a potential pan-European flood frequency distribution (Salinas et al., 2013a).Furthermore, the EVI has proven to give satisfactory and acceptable results for return periods less than 50 and 100 years, respectively, in estimating hydrometeorological extremes (Koutsoyiannis, 2004).

Application and results
The daily streamflow of the five study watersheds was simulated using the two proposed methodologies for ungauged watersheds and poorly gauged watersheds.The simulated and observed hydrographs compared graphically and statistically.Five statistical indices were used to assess the accuracy and performance of the two simulation methods, namely the NSE; the percent runoff volume error %DV = V sim −V obs V obs ×100; the correlation coefficient (CORR) between the simulated and the observed flows; RMSE (in m 3 s −1 ) between the simulated and the observed flows, and the average percent error of the maximum annual flows, where MaxQ sim j is the simulated maximum annual flow of year j , MaxQ obs j is the observed maximum annual flow of year j , and k is the number of hydrological years of the simulation period.
The model efficiency (NSE) is widely used in hydrological simulation studies.It compares the scale and the shape of the simulated and the observed hydrographs, and its optimal value is 1.The percent runoff volume (%DV) is a scale parameter which measures the percent error in volume under the observed and the simulated hydrographs for the period of simulation.Positive values of %DV indicate overestimation of the observed runoff volume, negative values of %DV indicate underestimation of the observed runoff volume, and %DV equal to zero indicate perfect agreement between simulated and observed runoff volumes.The correlation coefficient (CORR) is a shape statistical parameter that measures the linear correlation between the observed and simulated flows with optimal value of 1.The RMSE measures the residual or error variance between the simulated and the observed flows, and its optimal value is 0. The average percent error of the maximum annual flows (%AMAFE) estimates the average percent error in the simulation of the maximum annual peak flows for the simulation period.Positive values of %AMAFE show the average overestimation of the maximum annual flow, whereas negative values indicate the average underestimation of the maximum annual flow; its optimal value is 0.
Firstly, the five study watersheds were treated as ungauged and the UBCREG methodology for ungauged watersheds was applied.The daily streamflows of the study watersheds were simulated using the uncalibrated UBC watershed model with the estimated values of model parameters presented previously.The results of these simulations are shown in Fig. 3 and Table 4.The simulation was performed for the whole period of available data in each study watershed since the UBC watershed model was uncalibrated, and thus the whole simulation period is a validation period for the performance of the method.However, the training and validation periods indicated in Fig. 3 and Table 4 are indicated for comparison with the results of the second methodology intended for use in poorly gauged watersheds with limited streamflow measurements.
The graphical and the statistical comparison of the simulated hydrographs with the observed hydrographs (Fig. 3 and Table 4) show that, in general, the ungauged UBCREG method estimates the observed hydrograph with reasonable accuracy.For the Upper Campbell watershed, the value of CORR (CORR = 0.84) indicates that the method reproduced the shape of the observed hydrograph reasonably well but the annual peak streamflows were severely underestimated (%AMAFE = −32.06% in Table 4).The method performed better in the Illecillewaet watershed, for which there was a significant improvement in the simulation of hydrograph (NSE = 0.84 and CORR = 0.96 in Table 4).However, in the Illecillewaet watershed, the method overestimated the total runoff volume and the maximum annual peak flows (%DV = 14.6 3% and %AMAFE = 11.26 % in Table 4).The simulation results for the Yermasoyia watershed indicate that the method reproduced the shape and scale of the hydrograph reasonably well(NSE = 0.73 and CORR = 0.87 in Table 4) but overestimated the runoff volume and the annual peak discharge (%DV = 11.45 % and %AMAFE = 9.85 % in Table 4).The overall worst simulation results were acquired in the Astor watershed; however, the annual peak flows were generally overestimated (%AMAFE = 6.3 %), and the runoff volume was underestimated (%DV = −7.68%), leading to a low but acceptable value of model efficiency (NSE = 0.68) (Table 4).On the other hand, the best simulation results were found for the Hunza watershed.The statistical indices (Table 4) and the graphical comparison of the simulated and the observed hydrographs (Fig. 3) indicate that the shape and scale of the observed hydrograph were reproduced reasonably well.
The above results indicate that the simulation accuracy heavily depends on the quality and availability of meteorological data.This is obvious from the simulation results for the Illecillewaet watershed (Fig. 3b and Table 4).This watershed has three high-quality meteorological stations, and the hydrograph shape was simulated with improved accuracy, although the runoff volume and the annual peak flows were overestimated (Table 4).The performance of the method also seems to be dependant on the runoff generation mechanisms.As a comparison, better simulation results have discovered for watersheds that the runoff is mainly generated by snowmelt and glacier melt and not by watersheds where rainfall runoff is the dominant runoff generation mechanism.For example, the runoff simulation statistics for the Yermasoyia watershed is similar to the simulation statistics for the Upper Campbell watershed, although data from three precipitation stations were used for streamflow simulation of the small Yermasoyia watershed (157 km 2 ) and only one precipitation station was used in the Upper Campbell watershed, which is larger in area (1194 km 2 ).Furthermore, the best simulation results have been achieved for the Hunza and Illecillewaet watersheds (13 100 and 1150 km 2 in area, respectively).The runoff in the Yermasoyia watershed is generated by rainfall, whereas snowmelt is a significant percentage of total runoff in the Upper Campbell watershed.On the other hand, more than 90 % of the runoff in the Hunza Basin is generated by glacier melting, whereas snowmelt and glacier melt produces most of the runoff in the Illecillewaet watershed.The spatial variability of rainfall is much larger than the variability of snowfall.Also, the precipitation gradients are steeper and more consistent for snowfall than rainfall (Loukas andQuick, 1994, 1995).Hence, a larger number of precipitation stations is necessary in watersheds where rainfall-runoff is the dominant runoff generation mechanism in order to capture the spatial variability of rainfall and better simulate the streamflow (Brath et al., 2004).However, keeping in mind the very limited number of meteorological stations and data used, the overall results of the UBCREG methodology are judged satisfactory and show that the UBC watershed model can simulate reasonably well the watershed streamflow in various cli- matic and hydrological regions with a universal set of model parameters and making assumptions of precipitation stations representativeness and precipitation distribution.
The second proposed UBCANN methodology for poorly gauged watersheds was applied to the five study watersheds, assuming that only two or three years of daily streamflow data were available.The UBC watershed model was firstly run as in the first methodology for the years that streamflow data were assumed to be available, and the calculated runoff components were used as input to ANNs.The ANNs were optimized and trained for this initial period and then the UBC watershed model coupled with the trained ANNs was run and validated for the remaining period for validation.The final geometry or architecture of the optimized ANNs for the five  study watersheds is presented in Table 5. Figure 3 and Table 6 present the simulation results for the training and validation periods of the UBCANN methodology at the five study watersheds.Comparison of the graphical (Fig. 3) and statistical results (Tables 4 and 6) indicates that the coupling of UBC watershed model with ANNs greatly improves the simulation of hydrographs and maximum annual streamflow in all five watersheds compared to the first methodology.The discussion will be focused on comparison of the validation periods of UBCANN application since the ANNs of this methodology were optimized during the training period and an improvement in the simulation results is expected.Furthermore, to investigate the suitability of the UBCANN method for poorly gauged watersheds, the classical calibration method of the hydrological model is applied and compared.Table 7 presents the results of the UBCCLA method as a benchmark model for watersheds with limited information.
The simulation results of the UBCANN method for Upper Campbell watershed indicate that although there is significant improvement in the prediction of runoff volume and maximum annual peak flows (Table 6), the model efficiency (NSE = 0.68) has the same value with the first method (Table 4).On the other hand, the runoff simulation is greatly improved in the other four study watersheds.All statistical indices of the hydrological simulation have been improved in the Illecillewaet, Yermasoyia, and Astor wa-  tersheds (Table 6).Only the percent runoff volume error (%DV = −11.26% in Table 6) is not improved compared to the results of the UBCREG method (%DV = 0.25 % in Table 4) for the Hunza watershed.The improvement of the hydrograph simulation leads to better estimation of runoff volume and peak streamflow.The improvement of runoff simulation with the second methodology depends upon the volume and the range of the available streamflow data, since ANNs are a data intensive technique.When the available data cover a large range of streamflows, then the trained ANNs can accurately and efficiently simulate the unknown streamflows.
Application of the UBCCLA method shows that UBC is a reliable hydrological model in streamflow modelling for diverse climatic environments because the statistics are improved using streamflow data for calibration (Table 7).However, from Tables 6 and 7 it is difficult to assess the superiority of the UBCANN method using the UBCCLA method.For example, the validation NSE values show that the UBCANN method in the Yermasoyia and Astor watersheds greatly outperforms the UBCCLA method, in the Upper Campbell and Illecillewaet watersheds is marginally inferior to the UBCCLA method, and in the Hunza watershed both methods perform similarly (Tables 6 and 7).These contradictory results are also in agreement with the study of Anctil et al. (2004b), which showed that similar results are with the bias-corrected and accelerated method, which adjusts for both bias and skewness in the bootstrap distribution.The calculation procedure of these figures is described analytically in Ritter and Muñoz-Carpena (2013).Careful examination of scatterplots, NSE classes, and 95 % CI of the selected evaluation metrics NSE and RMSE showed that the UBCANN method is less effective in streamflow modelling than the UBCCLA in two watersheds (Figs. 4 and 5), whereas in the other three watersheds is superior to the UBC-CLA method .For these watersheds no prior information was used for the distribution of precipitation distribution and ANNs, with input the UBC flow components showing great skill in reproducing the daily streamflow patterns.However, in cases where prior hydrological knowledge was incorporated in the UBC model, such as in the two Canadian watersheds, ANNs showed similar capabilities with UBC-CLA approach due to expert knowledge "optimization" of the ungauged UBC flow components.The second step of the analysis was to compare the simulated and observed maximum annual peak flows and to perform a simple frequency analysis using the EVI theoretical distribution.It should be noted that the EVI distribution was selected to demonstrate that the methods employed for ungauged and poorly gauged watersheds and other candidate distributions could be used.This analysis was performed only for the four study watersheds (Upper Campbell, Illecillewaet, Yermasoyia, and Astor) which have streamflow data for at least six consecutive years.Application of the non-parametric Kolmogorov-Smirnov test for checking the adequacy of the selected distribution with the observed and simulated values showed that the EVI distribution is acceptable at the 5 % significance level for all observed and simulated streamflow values at the study watersheds.Figure 9 shows the comparison of the fitted EVI distributions using the three methodologies (UBCREG, UBCANN, and UBC-CLA) with the observed data and the fitted observed EVI for the four study watersheds.For Upper Campbell watershed these results indicate that the methodology for ungauged watersheds underestimates the observed maximum annual peak flows.Comparison of the UBCANN and UBCCLA methods for flood frequency estimation in poorly gauged basins showed that high peak flows are more accurately represented by the UBCANN method (Table 8 and Fig. 9a).Peak flow frequency analysis for Illecillewaet watershed indicates that the UBCREG methodology overestimate the observed peak flows.The best flood frequency curves for this watershed are acquired with the use UBCANN method, whereas the UBC-CLA method underestimates the peak flows for all examined return periods (1-100 years) (Table 8 and Fig. 9b).Peak flow frequency analysis for the poorly gauged Yermasoyia watershed again shows the superiority of the UBCANN method compared to the UBCCLA method.Flood frequency analysis of the UBCREG method suggests that caution is required for flood modelling since the method significantly underestimates the observed peak flows (Table 8 and Fig. 9c).Finally, in the Astor watershed, all applied methods perform similarly and the flood frequency estimation using simulated values underestimates the observed flows at larger return periods (Table 8 and Fig. 9c).However, except for the maximum annual peak of the last hydrological year of record 1996-1997 (Fig. 3), the simulated peak flows using the methodology for ungauged watershed underestimate the observed peak flows.For this particular year, the method severely overestimates the maximum annual peak flow.The result is that the estimated peak flows with return periods of 25, 50, and 100 years are quite similar with the applied methods for poorly gauged watersheds (Table 8).Overall the coupling of ANNs with the ungauged UBC flow model components is considered an improvement and an alternative method over the conventional calibration of a hydrological model with limited streamflow information based on the evaluation criteria employed for streamflow modelling and flood frequency estimation.

Conclusions
Rainfall-runoff modelling procedures for ungauged and poorly gauged watersheds are developed in this study.A well-established hydrological model (Singh, 2012), the UBC watershed model, is selected and applied in five different river basins located in Canada, Cyprus, and Pakistan.Catchments from cold, temperate, continental, and semiarid climate zones are included to demonstrate the procedures developed.Two methodologies for the modelling of streamflow are proposed and analysed.The first methodology, proposed for ungauged watersheds, uses the UBC watershed model with a set of universal constant values of model parameters and makes assumptions and estimates regarding the representativeness of precipitation stations and precipitation distribution.This methodology requires a good description of the watershed (area, elevation bands, vegetation coverage, soils, etc.) and limited meteorological stations as well as data to estimate the distribution of precipitation over the elevation range of the watershed, or even regional information about the orographic precipitation gradients of a watershed.The second methodology is an extension of the first method, and couples the UBC watershed model with ANNs.This method is proposed for poorly gauged watersheds.The limited streamflow data are intended for training of ANNs.For comparison purposes, this method is compared with the classical calibration of the UBC model in poorly gauged watersheds.The evaluation of all the applied methods is based on the combination of graphical results, statistical evaluation metrics, and normalized goodness-of-fit statistics.
Application of the methods employed to five watersheds with areas that are in the range of 157 to 13 100 km 2 , have different runoff generation mechanisms, and are located in various climatic regions of the world resulted in reasonable results for the estimation of streamflow hydrograph and peak flows.The first methodology for ungauged watersheds performed quite well despite the very limited available meteorological data.The second hybrid method is a significant improvement on the first methodology because it takes advantage of the limited streamflow information.The coupling of the UBC regional model with ANNs provides a good alternative to the classical application (UBC calibration and validation) without the need for optimizing UBC model parameters.The ANNs coupled to the UBC watershed model improve the streamflow modelling at poorly gauged basins.Furthermore, using the non-calibrated UBC flow components as input to ANNs in a coupling or hybrid procedure combines the flexibility and capability of ANNs in nonlinear modelling with the conceptual representation of the rainfall-runoff process acquired by a hydrological model.Hence, the black-box constraints in ANN modelling of the rainfall-runoff are minimized.Overall the coupling of ANNs with the regional UBC flow model components is considered to be a successful alternative method over the conventional calibration of a hydrological model with limited streamflow information based on the evaluation criteria employed for streamflow modelling and flood frequency estimation.In the future, the two methodologies should be compared with other regional techniques or hydrologic models and could be applied in other regions to generalize the results.Another step further could be the more rigorous estimation of flood frequency by additionally incorporating the uncertainty of the state variables.

Figure 1 .
Figure 1.Flow diagram of the UBC Watershed model.

Figure 2 .
Figure 2. Typical ANN geometry for combining the outputs of the UBC watershed model in the methodology for poorly gauged watersheds.
The developed ANNs were operated in batch mode, which means that the training sample presented to the network between the weight updates was equal to the training set size.This operation forces the search to move in the direction of the true gradient at each weight update; however, it requires large storage.The mean squared error was used as the minimized error function during the training.The initial values of weights for each node were set to a value, a = 1

Table 1 .
Characteristics of the five study watersheds.
PMXIMP is the percentage of impermeable area of the elevation zone and is estimated by Eq. (5); RNSM is the summation of rainfall, snowmelt, and glacial melt of the time step; V0FLAS is a parameter showing the threshold value of precipitation for flash runoff; and V0FLAX is the parameter showing the maximum value of precipitation, which limits the FMR range.The last two parameters (i.e. ) Nat. Hazards Earth Syst.Sci., 14, 1641-1661, 2014 www.nat-hazards-earth-syst-sci.net/14/1641/2014/

Table 3 .
Default values for the water allocation and flow-routing parameters of UBC watershed model.

Table 4 .
Statistical indices of streamflow simulation with the proposed methodology for ungauged watersheds -UBCREG method.

Table 5 .
Geometry of optimized ANNs used in the methodology for poorly gauged watersheds.

Table 6 .
Statistical indices of streamflow simulation with the proposed methodology for poorly gauged watersheds -UBCANN method.

Table 7 .
Statistical indices of streamflow simulation with the classical methodology for poorly gauged watersheds -UBCCLA method.