Predicting streamflows in snowmelt-driven watersheds using the flow duration curve method

Predicting streamflows in snow-fed watersheds in the Western United States is important for water allocation. Since many of these watersheds are heavily regulated through canal networks and reservoirs, predicting expected natural flows and therefore water availability under limited data is always a challenge. This study investigates the applicability of the flow duration curve (FDC) method for predicting natural flows in gauged and regulated snow-fed watersheds. Point snow observations, air temperature, precipitation, and snow water equivalent were used to simulate the snowmelt process with the SNOW-17 model, and extended to streamflow simulation using the FDC method with a modified current precipitation index. For regulated watersheds, a parametric regional FDC method was applied to reconstruct natural flow. For comparison, a simplified tank model was used considering both lumped and semi-distributed approaches. The proximity regionalization method was used to simulate streamflows in the regulated watersheds with the tank model. The results showed that the FDC method is capable of producing satisfactory natural flow estimates in gauged watersheds when high correlation exists between current precipitation index and streamflow. For regulated watersheds, the regional FDC method produced acceptable river diversion estimates, but it seemed to have more uncertainty due to less robustness of the FDC method. In spite of its simplicity, the FDC method is a practical approach with less computational burden for studies with minimal data availability.


Abstract.
Predicting streamflows in snow-fed watersheds in the Western United States is important for water allocation. Since many of these watersheds are heavily regulated through canal networks and reservoirs, predicting expected natural flows and therefore water availability under limited data is always a challenge. This study investigates the applicability of the flow duration curve (FDC) method for predicting natural flows in gauged and regulated snow-fed watersheds. Point snow observations, air temperature, precipitation, and snow water equivalent were used to simulate the snowmelt process with the SNOW-17 model, and extended to streamflow simulation using the FDC method with a modified current precipitation index. For regulated watersheds, a parametric regional FDC method was applied to reconstruct natural flow. For comparison, a simplified tank model was used considering both lumped and semi-distributed approaches. The proximity regionalization method was used to simulate streamflows in the regulated watersheds with the tank model. The results showed that the FDC method is capable of producing satisfactory natural flow estimates in gauged watersheds when high correlation exists between current precipitation index and streamflow. For regulated watersheds, the regional FDC method produced acceptable river diversion estimates, but it seemed to have more uncertainty due to less robustness of the FDC method. In spite of its simplicity, the FDC method is a practical approach with less computational burden for studies with minimal data availability.

Introduction
Snow accounts for a significant portion of precipitation in the mountainous Western United States and snowmelt plays an important role in forecasting streamflow (Serreze et al., 1999). Extreme amounts of snowfall can result in a flood in the melting season, and sometimes snow accumulation alleviates drought by natural redistribution of precipitation in a high water-demand period. In such regions, snowmelt controls the hydrologic processes and water relevant activities such as irrigation. Therefore, the reliable prediction of snowmelt is crucial for water resources planning and management (He et al., 2011;Mizukami et al., 2011;. Conventionally, conceptual snowmelt models developed by combining rainfall-runoff models with temperature index models using a parameterized melting factor (e.g., Anderson, 2006;Albert and Krajeski, 1998;Neitsch et al., 2001), have been used to predict daily streamflows in snow-fed watersheds. Conceptual modeling is an attractive solution to daily streamflow simulation not only for rainfall-fed but also for snow-fed watersheds due to its flexibility and applicability (Uhlenbrook et al., 1999;Smakhtin, 1999). Examples include models such as SSARR (Cundy and Brooks, 1981), PRMS (Leavesley et al., 1983), NWSRFS (Larson, 2002), UBC (Quick and Pipes, 1976), CEQUEAU (Morin, 2002), HBV (Bergström, 1976), SRM (Martinec, 1975), and TANK (Sugawara, 1995), among others.
However, a significant simplification is necessary when complex hydrological behavior of a watershed is implicitly parameterized into a conceptual model . Such simplifications make it difficult to relate model parameters directly to measured watersheds properties (Beven, 2006). Hence, the parameters of conceptual models are usually identified by streamflow observations with calibration techniques such as the shuffled complex evolution or genetic algorithm. In truth, calibration is the major part of conceptual modeling, and it is still typically labor-consuming; however, computational efficiency has improved with advances in computer technology. In spite of the effort involved, uncertainty in conceptual models is always an important issue (Kuczera and Parent, 1998;Uhlenbrook et al., 1999;Panday et al., 2013). Furthermore, the parameter set calibrated by streamflow observations is usually not unique because there can be other sets of parameters providing similar model performance (Beven, 1993;Seibert, 1997;Oudin et al., 2006;Perrin et al., 2007). Particularly in snowmelt runoff modeling, calibration can produce less uniqueness, less robustness, and more uncertainty than rainfall-runoff modeling because additional inputs (e.g., air temperature) and parameters (e.g., melting factor) are required to define the snowmelt process.
As an alternate approach, linking point snow observations to streamflow can be a pragmatic option. A common statistical approach for simple generation of daily streamflow is the flow duration curve (FDC) method. A FDC gives a summary of streamflow variation and represents the relationship between streamflow and its exceedance probability (Vogel and Fennessey, 1994). For streamflow generation, one or multiple sets of donor variables are transferred to a target station by corresponding exceedance probability of the donor sets with that of the target. A number of variations of the FDC method have been used for the generation of daily streamflow data. Hughes and Smakhtin (1996), for instance, suggested a FDC method with a nonlinear spatial interpolation method to extend observed flow data. Smakhtin and Masse (2000) developed a variation of the FDC method to generate streamflow using rainfall observations as the donor variable instead of streamflow data. Recently, the FDC was used not only for generating streamflow directly, but also for calibrating conceptual models (Westerberg et al., 2011). Westerberg et al. (2011) used the FDC as a performance measure to circumvent uncertainty in discharge data and other drawbacks in model calibration with traditional methods. Despite the numerous applications with the FDC, there is still no good approach using the FDC method to generate daily streamflow from point snow observations. Given the simplicity of the FDC method, a suitable approach using the FDC method to predict snowmelt-driven runoff using point snow observations could be practical and cost efficient due to the reduced computational effort.
If the target station is ungauged, a regional FDC can estimate the FDC of the target station. The regional FDC is generally developed using the relationships between selected percentile flows in gauged FDCs and climatic or physical properties of the watersheds. Thus, the regional FDC estimates the unknown FDC of an ungauged watershed only with its physical properties. Many regional FDC methods have been proposed for generating streamflows in ungauged wa-tersheds. Shu and Ouarda (2012) categorized the regional FDC methods as a statistical approach (e.g., Claps et al., 2005), a parametric approach (e.g., Yu et al., 2002;Mohamoud, 2008), and a graphical approach (e.g., Smakhtin et al., 1997).
The regional FDC can be used not only for generating streamflows in ungauged watersheds, but also for reconstructing natural flows of watersheds regulated by reservoir operations, river diversions and other human activities. Smakhtin (1999), for example, evaluated the impact of reservoir operations by comparing between regulated outflows from a reservoir and natural flow estimated by a regional FDC. In the Western United States, the prior appropriation doctrine, the water right of "first in time, first in right," has produced many river basins with impaired streamflows. These impairments are particularly significant in watersheds with high aridity, low precipitation, and relatively large water demands. The regional FDC method can represent flow impairments by reconstructing natural flows using minimal data. The reconstruction of natural flow provides additional information to water managers for efficient water allocation during the high-demand periods. The volume difference between reconstructed natural flows and impaired streamflow observations can simply indicate the combined effects of reservoir operations, river diversions, and other humandriven activities. Thus, the effect of regulation in a watershed can be approximately evaluated from this comparison.
As discussed earlier, prior studies using the FDC method with precipitation data focused on predicting streamflows in natural and managed watersheds under typical rainfallrunoff conditions and not with snowmelt-driven streamflow. Therefore the goals of this work are twofold: first to assess the applicability of the FDC method in predicting streamflows in semi-arid snowmelt-driven watersheds through the comparison with conceptual rainfall-runoff models incorporating a temperature index-based snowmelt model; and second to assess the possibility of extending the work through regionalization to predict natural streamflows in regulated watersheds to determine water availability. In this work, a modified approach to the FDC method for streamflow generation from rainfall observations (Smakhtin and Masse, 2000) is proposed. The simplified SNOW-17 model was used here with point snow observations to estimate snowmelt discharge required by the FDC method and the conceptual model. Also, a parametric regional FDC method was applied for the reconstruction of natural flows and a proximity-based regionalization approach was used in the conceptual rainfall-runoff models for comparison with the regional FDC. By comparing with impaired streamflows and observed managed flows, water use in a watershed was estimated.

Description of the study area and data
The study area is the Sevier River basin, located in South Central Utah, and the details are given Fig. 1. The Sevier River basin is a semi-arid basin with relatively high ET (evapotranspiration). The watersheds in or adjacent to the Sevier River basin are dominantly fed by snowmelt from the high-elevation region. Particularly, the Sevier River is significantly regulated by diversions and reservoir operations along the major channel for agricultural water use. Hence, a realtime streamflow monitoring system along the main channel is operated by the Sevier River Water Users Association, but it is difficult to estimate the natural discharge from the regulated watersheds using this monitoring system. This study used the US Geological Survey (USGS) streamflow stations for the FDC method and conceptual modeling. Because only five watersheds in the Sevier River basin have natural streamflow observations, eight adjacent watersheds were included as well for generating streamflows in gauged watersheds. In addition, two USGS stations in the main Sevier River with significant impairments were selected for reconstructing natural flows using the regionalization methods. These two stations were assumed as ungauged watersheds although these have continuous daily observations. Hence, "gauged" watersheds in this study refer to watersheds with natural flow observations only, while "regulated" watersheds indicate watersheds with impaired flows and therefore these watersheds are treated as ungauged watersheds.
Precipitation, maximum and minimum air temperature, and snow water equivalent (SWE) data from the SNOTEL stations operated by US Department of Agriculture (USDA) were used as inputs to the FDC method and conceptual modeling. The details of the USGS stations and corresponding SNOTEL stations are given in Table 1 with corresponding data periods and watershed areas. Additionally, the records of canal diversions from the Utah Division of Water Rights were used to compare streamflows simulated by regionalization with actual river diversions. For the conceptual modeling, point SNOTEL data were adjusted to spatially averaged inputs using data from the PRISM database (PRISM Climate Group, 2012). The procedure included a comparison between a pixel located in a SNOTEL station and the areal average of pixels in a watershed or an elevation zone using 30 arcsec annual normals from 1981 to 2010. The ratio of the average of pixels to the pixel at a SNOTEL station was multiplied by the point precipitation at the SNOTEL station, while the difference between these was added to the point temperature. For the regional FDC, the SNOTEL data adjusted by PRISM data were also used for calculating climatic variables. The USGS National Elevation Dataset (2012) and US General Soil Map served by USDA (2013) were used to obtain geomorphologic and soil properties of the watersheds.

SNOW-17 snowmelt model
This study uses SNOW-17 as the snowmelt model which has been used for river forecasting by the National Weather Service (NWS). SNOW-17 is a single-layered, conceptual snowmelt model. This model estimates SWE and snowmelt depth as outputs. Input data required are precipitation and air temperature only. Although the original SNOW-17 model has 10 parameters for point-scale simulation, this study used the simplified model similar to Raleigh and Lundquist (2012). For simplification, temperature for dividing rainfall and snowfall (PXTEMP), base temperature for non-rain melt (MBASE), and the liquid water holding capacity (PLWHC) were assumed at typical values of 1.5 • C, 0 • C, and 5 %, respectively. Rain on snowmelt and daily melt at the snow-soil interface were deactivated since these contribute minimally to the energy budget of the snowmelt process (Raleigh and Lundquist, 2012;Walter et al., 2005). The simplified version has only five parameters, which are SCF, MFMAX, MFMIN, NMF, and TIPM. SCF is a multiplying factor to adjust new snow amounts. MFMAX and MFMIN are the maximum and minimum melting factors to calculate melting depths, respectively. NMF and TIPM are parameters for simulating energy exchange when there is no snowmelt. A detailed description of the model was given by Anderson (2006). This study used Nash-Sutcliffe Efficiency (NSE) for performance evaluation of SNOW-17 and model where Q SWE (t) andQ SWE (t) are observed and simulated SWEs (mm) at time t, respectively, Q SWE is the mean observed SWE (mm), and T is the number of observations.

Modified FDC method with precipitation index
The FDC method is a non-parametric probability density function representing the relationship between magnitude of streamflow and its exceedance probability. The FDC method is typically used to generate daily streamflow at a station from highly correlating donor streamflow data sets with a target station. A drawback of this approach is that streamflow generation is dependent on the availability of donor data sets. Hence, in a region with a low density of stream gauging stations, the FDC method may face the difficulty of not having adequate donor streamflow data. Smakhtin and Masse (2000) developed a modified FDC method with a precipitation index to overcome the limited availability of donor variable sets. Their method included transforming the time series of precipitation into an index having similar properties to streamflow data. The transformation was to avoid zero values in precipitation data caused by the intermittency of precipitation events, which therefore produce a different shape of duration curve from a typical FDC. The duration curve of transformed precipitation could indicate the exceedance probability at the outlet, which determines the magnitude of streamflow.
This study modified the original concept as follows. First, the outflow depth simulated by SNOW-17 was used for constructing the FDC instead of precipitation data to represent the snowmelt process. Second, a constant recession coefficient was applied for the calculation of precipitation index of Smakhtin and Masse (2000), but different coefficients were used to represent the different hydrologic responses of rainfall and snowmelt to streamflow. The modified approach is given below.
The current precipitation index at time t, I CP (t) in mm d −1 was defined in the original work as where k is the recession coefficient (d −1 ), P (t) is daily precipitation at time t (mm d −1 ), and t is the time interval (d).
Recession coefficient, k, represents the similar concept to the baseflow recession coefficient and needs to be determined by observed streamflow. According to previous studies, k varies from 0.85 to 0.98 d −1 (Linsley et al., 1982;Fedora and Beschta, 1989). In addition, the initial value of I CP can be assumed as the long-term mean daily precipitation because of the fast convergence of calculations (Smakhtin and Masse, 2000).
To consider the snowmelt process, outflow calculated by SNOW-17 was divided into two time series, since it was important to stipulate different recession coefficients for snowmelt and rainfall processes given the different timesscales of these processes for generating streamflow (De-Walle and Rango, 2008). Time series of snowmelt depth and rainfall depth were separated based on the existence of snow cover (when SWE > 0). Finally, the two indices were summed for simulating I CP . Hence, the I CP (t) is redefined as where I CS (t) is the current snowmelt index (mm) at time t, S(t) is the snowmelt depth (mm) at time t, I CR (t) is the current rainfall index (mm) at time t, R(t) is the rainfall depth (mm) at time t, k S and k R are recession coefficients (d −1 ) for snowmelt and rainfall, respectively. Generally, k S is greater than k R because snowmelt runoff varies more smoothly with time than quick flow caused by rain storms. In this study, k S and k R were selected by values showing maximum correlation between I CP and observed streamflow data. Figure 2 shows the proposed FDC method used in this work. The selection of a snow observation station when multiple stations are present in a watershed was based on high correlation between calculated I CP and observed streamflow. Although Smaktin and Masse (2000) commented that the effect of weights in the case of multiple stations was not a significant factor in their original FDC method with the precipitation index, a high correlation between I CP and streamflow supports better performance in the generation of streamflow because of the significant climatic variation of snow-fed watersheds located in high-elevation regions.

Simplified tank model
This study used the simplified tank model proposed by Cooper et al. (2007) to compare the performance under the conditions of similar and limited data availability. The simplified tank model reduced the number of parameters of the original tank model (Sugawara, 1995) to help minimize overparameterization when the tank model was combined with the snowmelt model. This simplified tank model shown in Fig. 3a has two vertical layers with the primary soil moisture layer in the upper tank. This study did not consider the secondary soil moisture layer in the simplified tank model because it was not sensitive to runoff simulations (Cooper et al., 2007). Evapotranspiration (ET) in the tank model was independently estimated using the modified complementary method proposed by Anayah (2012). The combined model has 12 parameters (5 for snowmelt, 7 for runoff). The structure of the tank model is adequately flexible to be calibrated by streamflow observations. It has more parameters than the Snowmelt Runoff Model with eight parameters (Martinec et al., 2008).
The model produces several modes of response representing the different conditions that may prevail in a watershed. The upper tank has a non-linear response in the rainfall-runoff process because of its multiple horizontal outlets, whereas the lower tank has a linear response. There are three thresholds to determine the four modes of hydrologic response, which are H S , H 1 , and H 2 . H S represents the soil moisture-holding capacity (mm). H 1 and H 2 represent the lower and upper thresholds for generating direct runoff (mm). The detailed procedure for calculating streamflow is available from Cooper et al. (2007).
This study used two approaches with the proposed tank model (as depicted in Fig. 3) for evaluating the performance with and without the consideration of climatic variation in a watershed. The first approach was a completely lumped model with a single set of climatic inputs that disregards the climatic variation of a watershed (Fig. 3a). The second approach was a semi-distributed tank model with five different tanks for the upper layer to accommodate climatic variation due to elevation (Fig. 3b). All of the upper tanks in both approaches were assumed to have same parameters for both snowmelt and runoff modeling. For the semi-distributed tank model, a watershed was divided into five zones with the aid of the area-elevation relationship. Inputs for each zone were individually computed from the corresponding SNOTEL station and PRISM data as explained earlier.
The parameters were optimized using the genetic algorithm in Matlab for both the lumped and the semi-distributed tank models with the objective function of minimizing the sum of weighted squared residuals shown as below.
where w(t) is weight (unitless) varying with magnitude of runoff data, Q(t) andQ(t) are observed and simulated streamflows (m 3 s −1 ), respectively, and T is the number of observations. The weights can be determined empirically with observed data for equalizing residuals in low flows with those in high flows. The weights used in previous studies (e.g., Kaluarachchi, 2008, 2009) ranged from 4 to 10. The average streamflows of gauged watersheds in the high flow season (April to June) were about 2 to 10 times (with median of 5.17) those in the low flow season (March to June). Hence, this study used a weight of 5 for the low runoff season and 1 for the high runoff season. Although Cooper et al. (2007) proposed two constraints to calibrate the tank model parameters with wide ranges, incorporating SNOW-17 into the tank model made it difficult to apply the constraints to the combined model. Hence, in the optimization with genetic algorithm, the ranges of parameters were identified using Monte Carlo simulations with uniform distributions. One of the best 100 parameter sets obtained by sorting the values of the objective function was selected to set the parameter ranges for genetic algorithm.

Regionalization
This study applied regionalization to simulate natural streamflows in regulated watersheds with impaired observations. A parametric approach was selected for constructing the regional FDC. The model proposed by Shu and Ouarda (2012) was used and given as where Q P is percentile flows, V 1 , V 2 , V 3 ,. . . are selected physical or climatic descriptors, b, c, d,. . . are model parameters, and a is the error term. Logarithmic transformation of Eq. (5) can help solve the model through linear regression. By step-wise regression, independent variables can be selected.
Meanwhile, a proximity-based regionalization method was used for the tank model. In the case of conceptual modeling, regionalization of parameters for ungauged watersheds were categorized by three approaches (Peel and Blöschl, 2011): (a) regression analysis between individual parameters and watershed properties (e.g., Kim and Kaluarachchi, 2008;Gibbs et al., 2012); (b) parameter transfer based on spatial proximity (e.g., Vandewiele et al., 1991;Oudin et al., 2008); and (c) physical similarity (e.g., McIntyre et al., 2005;Oudin et al., 2008Oudin et al., , 2010. Even if the performance of these three approaches was dependent on climatic conditions, performance and complexity of the model, and other factors, several studies concluded that the spatial proximity method was attractive due to its better performance and simplicity (Oudin et al., 2008;Parajka et al., 2013). Hence, this study used the proximity-based regionalization for regulated watersheds. Parameter sets were transferred from multiple gauged watersheds for better precision, and the average of streamflows simulated by the parameter sets was taken as the natural flow estimates for the regulated watersheds.

SNOW-17 modeling
SNOW-17 was calibrated and verified by SWE observations at SNOTEL stations. Figure 4 shows the results of SNOW-17 modeling where the comparison between simulated and observed SWE is excellent. The average NSE values between simulated and observed SWE for calibration and validation were 0.942 (a range of 0.867 to 0.984) and 0.933 (a range of 0.793 to 0.967), respectively. The loss of NSE from calibration to validation was not significant and therefore the model was unlikely to be over-parameterized. Also, the simple objective function of maximizing NSE (equivalent to minimizing the sum of squared residuals) seems to provide adequate performance as long as accumulated precipitation shows a consistent trend with observed SWE in the snow accumulation period. Simultaneous monitoring of precipitation and SWE at the same location may provide quality inputs to SNOW-17 modeling.
However, a temperature index snowmelt model can have errors from strong winds and dew-point temperature (Anderson, 1976). In other words, good calibration by SWE observations does not necessarily guarantee accurate simulation of outflow depth. The loss of SWE by winds or sublimation, for instance, is not contributing to the melting depth while some SWE reduction is observed. Thus, in a region with high possibility of such errors, caution is required to link point snowmelt observations to streamflow.

Streamflow generation in gauged watersheds
The time series of outflow depth from SNOW-17 was used to calculate I CP . Since the rationale behind the FDC method is that exceedance probability of I CP is same as that of streamflow, the data periods of both point snow observations and streamflow data should be same. In fact, I CP calculation is mathematically equivalent to the computation of storage in a single linear reservoir such as the lower tank in the tank model. Hence, the hydrological meaning of I CP is liquid water availability in a watershed with the assumption of a single linear reservoir. Through the I CP computation, the intermittent time series of outflow depth was transformed to a smooth time series.
The computed recession coefficients of snowmelt varied from 0.97 to 0.98 d −1 , while the range for rainfall was 0.85 to 0.86 d −1 . These results demonstrate that snowmelt runoff was slowly changing during the year, unlike rainfall runoff that showed a relatively large fluctuation due to the intermittent storm events. In the study area, snowmelt runoff accounted for a large portion of streamflow and therefore the recession coefficient of snowmelt played a major role in the high correlation between I CP and streamflow. However, if there was noticeable contribution of rainfall runoff to streamflow observations, then the recession coefficient of rainfall would be more important and sensitive. Particularly, rainfall runoff can be crucial in the non-melting season, and therefore, the separation of recession coefficients is necessary for high correlation between I CP and streamflow.
When calibrating the lumped and semi-distributed tank models, Monte Carlo method was used to identify the parameter ranges of the tank model for optimization with genetic algorithm as commented earlier. The random simulations were to avoid local parameter sets providing unrealistic or poor streamflow simulation when using genetic algorithm with wide parameter ranges. To decide on the required number of simulations, the Clear Creek watershed was selected and tested among the given gauged watersheds. By increasing the number of simulations from 1000 to 20 000, it found that 20 000 simulations provided the efficient number of simulations with the initial parameter ranges. From the best 100 parameter sets of the 20 000 simulations, a parameter set with an acceptable NSE and a low reduction of NSE between calibration and validation was chosen. For optimization with genetic algorithm, the parameter ranges were rescaled with the ranges of approximately 50 to 200 % of each parameter of the chosen set. With the rescaled param- eter ranges, the genetic algorithm produced the optimal parameter set. It was later found that the optimal parameter set showed better performance than the best 100 parameter sets of the 20 000 simulations for all gauged watersheds. From this observation, the optimal parameter set was assumed as the calibrated parameter set.
As expected, the semi-distributed tank model performed better than the others with NSE, as shown in Table 2. Figure 5 depicts the simulated streamflow at several stations using the FDC method and the tank model. Due to the high climatic variation in mountainous watersheds, ignoring the elevation distribution could result in poor streamflow generation. These results confirmed the earlier studies (e.g., Martinec et al., 2008;Uhlenbrook et al., 1999) that discussed the importance of the elevation distribution on snowmelt runoff modeling. Theoretically, it is natural to expect poor performance from point snow observations of the FDC method and Creek, Beaver River, and Mammoth Creek, which had fairly high correlation between I CP and streamflow data, showed good performance in streamflow prediction. Even the semidistributed tank model did not show better results than the FDC method for Ferron Creek and Beaver River. Typically, watersheds showing good performance with the FDC method have good performance with the lumped and semi-distributed tank models too. Since both methods used linear reservoir coefficients for simulating streamflow, they performed well in watersheds with linear behavior and such watersheds were likely to have relatively homogenous climatic conditions. In addition, the FDC method showed the highest performance reduction from calibration to validation among the three methods. This may be due to the unstable correlation between I CP and streamflow and the uncertainty of the FDCs. Figure 6 shows a comparison between field discharge measurements and simulated streamflows in the calibration period. In order to avoid potential errors in streamflow observations converted from water stage, streamflow simulations by three methods were directly evaluated by field measurements. Table 3 summarizes the NSE and correlation coefficient values between field measurements and three simulations. Streamflow values for this evaluation were normalized by watershed area to remove the influence of watershed scale. On average, the performance trend from the poorest to the best watersheds was similar to the calibrations with the continuous streamflow data in terms of NSE. However, Vernon Creek and Salt Creek experienced a large re- duction of NSE when compared with field measurements. It means that these two watersheds had relatively large observational errors in the continuous streamflow data. In addition, Muddy Creek and Sevenmile Creek had better NSE for the lumped and semi-distributed tank models with field measurements. It also means the two watersheds possibly had considerable observational errors, but the conceptual models produced more precise streamflows than water stage data and rating curves. Also, Mammoth Creek, Sevier River at Hatch, and Coal Creek were likely to underestimate high flows with all three methods, but this was not experienced with continuous streamflow data. This indicates precipitation data for the three watersheds were also underestimated, or I CP and the model parameters were adapted by the underestimated high flows.

Regional FDC for regulated watersheds
The FDC method and the tank model were upscaled to watersheds affected by river diversions and reservoir operations to predict the natural flows at impaired streamflow stations. As mentioned earlier, regionalization was used for upscaling of regulated watersheds. The regulated station near the Piute Reservoir (Fig. 1) is Seveir River near Kingston, and the other near the Sevier Bridge Reservoir is Sevier River below San Pitch River near Gunnison (hereafter Sevier River near Gunnison).Water use in agricultural areas through river diversions significantly affect streamflow observations in the two stations. Streamflow observations at Sevier River near Kingston only include river diversions while the diversions and reservoir operations are included in streamflow observations at Sevier River near Gunnison. The two watersheds were divided into several sub-watersheds because these were too large to fall within the areas of gauged watersheds used for developing the regional FDCs. Hence, the sum of streamflows of each sub-watershed simulated by regionalization was the volume of natural flow at each target station. Climatic, geomorphologic, land cover and soil properties of the gauged watersheds were used to identify independent variables in determining the percentile flows of the parametric regional FDC. The candidate properties are listed in Table 4. The step-wise regression was implemented for each percentile flow in the Matlab environment. The variable with the largest significance among the candidates was taken as an independent variable for the first step. Then, other variables were added step by step based on the p value of F statistics. The selected variables for each percentile flow and the statistics of the regression analysis are given in Table 5. Overall, the regional FDC reproduced minimum, average, and stan-dard deviation well, but underestimated the maximum of percentile flows. This means the regional FDC may underestimate percentile flows of large watersheds; therefore it is not recommended to use the regional FDC for an ungauged watershed with an area larger than the largest watershed of the regression model.
As expected, watershed area was included in every percentile flow as an independent variable. Watershed area was positively related to percentile flows, and its multipliers ranged from 0.5 to 1.0. The multiplier had an increasing tendency as percentile increases. The routing effect on high flow (low percentile) may cause less proportionality to watershed area than low flow (high percentile).
Also, mean elevation was selected as another crucial independent variable. The multiplier of elevation varied from 2.2 to 3.7. Elevation was considered to be a geomorphologic property, but it represented the climatic variation of the watersheds because every climatic candidate had high correlation with elevation. It is a natural observation because more precipitation and lower air temperature are expected in the higher elevations.
Proportion of clay, dry bulk density, and saturated hydraulic conductivity were chosen to explain the variance of the regression errors remained from watershed area and mean elevation. The higher proportion of clay means lower permeability of soil, and saturated hydraulic conductivity controlling infiltration. Hence, the proportion of clay seems to affect high flows while saturated hydraulic conductivity was selected for low flows. The higher dry bulk density produces less porosity and less water-holding capacity in soils, thus a positive relationship was obtained between dry bulk density and 30 and 40 percentile flows. Drainage density was included as an additional significant variable for low flows with negative relationships. The negative relationship is probably because the higher drainage density means more distribution of streamflow in a watershed. When using the regional FDC approach, I CP was not necessarily used as the only donor variable to transfer exceedance probability to the target stations. In fact, the best donor variable is a data set that can show the best correlation with gauged streamflow at the target station. However, it is impossible to check the correlation between donor variables and ungauged streamflow. Thus, one or multiple donor variables close to the target station have been typically used in the regional FDC approaches. Shu and Ouarda (2012) suggested using multiple donor variables to minimize the uncertainty of using a single donor variable. This study used two sets of neighboring streamflow observations as well as I CP to generate streamflows in sub-watersheds. The recession coefficients of I CP were assumed to be 0.98 and 0.85 d −1 for snowmelt and rainfall, respectively. As commented earlier, parameters of both lumped and semi-distributed tank models were transferred from nearby gauged watersheds for streamflow simulation at the target stations. The parameter sets of Mammoth Creek, Sevier River at Hatch, Coal Creek, and Beaver River were used for Sevier River near Kingston while Salina Creek, Manti Creek, Ferron Creek, and Sevenmile Creek were selected for Sevier River near Gunnison. Figure 7 shows the simulated streamflows by the regional FDC and the tank models with regionalized parameters at both target stations. In the case of Sevier River near Gunnison, the outflow from the Rocky Ford Reservoir was subtracted from the observed streamflow to calculate the discharge produced by the watershed only. It could be easily recognized that these two watersheds were significantly regulated based on the irregular shapes of hydrographs. At Sevier River near Kingston, the regional FDC method estimated more volume of natural flow than the lumped and the distributed tank models. On the other hand, water volume estimated by the regional FDC was between the estimates of the lumped and semi-distributed models at Sevier River near Gunnison. Volume errors between the regional FDC method and the tank models varied from −17.1 to +21.8 %. The differences among the three methods were mainly in middle to high flows rather than low flows. The correlation coefficients between the simulations with the regional FDC and the lumped tank model were 0.94 and 0.70 at both stations, respectively, while those between the regional FDC and the semi-distributed tank model were 0.92 and 0.90, respectively. The larger difference between the lumped and semi-distributed models at Sevier River near Gunnison may be due to the higher climatic variation of this watershed, making the lumped assumption inappropriate. This is evident from the greater difference of NSE between the lumped and semi-distributed models of gauged watersheds transferred to Sevier River near Gunnison.

FDC method for gauged watersheds
The basis of the FDC method is point snowmelt modeling with SNOW-17. SNOW-17 performed well for the study area, but its parameter uncertainty could be a concern similar to conceptual runoff modeling. However, the five parameters used in SNOW-17 were small when compared to most classical hydrologic models. Indeed, a simpler snowmelt model (e.g., DeWalle and Rango, 2008) or observed snowmelt depth (equivalent to a reduction in observed SWE) could be an alternative for SNOW-17, while not necessarily reducing the uncertainty.
The performance of the FDC method was affected by the correlation between I CP and streamflow. Particularly, the correlation between I CP and middle to high flow determined the performance. Figure 8 shows the relationship between the performance and the correlation coefficient between I CP and streamflow with exceedance probability less than 0.2. Based on this knowledge, good performance (NSE > 0.8) could be expected when the correlation coefficient is greater than 0.8. The greater NSE in the validation period of Mammoth Creek and Sevier River at Hatch (Table 2) than in the calibration period could be explained by the correlation coefficient. These two watersheds had greater correlation coefficients (about 0.04 differences for both watersheds) in the validation period. The stable FDCs found for both watersheds also supported the better performance during validation.
It is also noted that the FDC method is not any more robust than the other methods. As shown in Table 2, the NSE of the FDC method has a much wider range from the poorest to the best performing watersheds than the others. Indeed, more watersheds showed better NSE, as the inputs were more distributed. This means that considering only point inputs with the FDC method could result in highly variable performance. Also, more distributed inputs would be better for more robust performance, even in the case of a simple model. With the FDC method, its low input requirement and computational burden has to be traded with some loss of robustness of performance.
In general, the FDC method had a poorer performance than the lumped and the semi-distributed tank models. One reason may be that the tank model was directly calibrated to streamflow observations, while the FDC method matched the magnitudes of I CP and streamflow based on an empirical probability density function. However, the main reason was that correlation between I CP and streamflow could be lower significantly from one period to another. Fish Creek, for instance, experienced a reduced correlation coefficient (about 0.35) from calibration to validation. On the other hand, the lumped and semi-distributed models that considered spatial variations did not have such large reductions in NSE. It means a point snow observation might not represent the behavior of an entire watershed. Hence, the first task is to assess the applicability of the FDC method by evaluating the correlation between I CP and streamflow.
There could be many reasons for the low correlation between I CP and streamflow. For example, Vernon Creek and Muddy Creek showed poor performances with the FDC method, but the reasons were different. Vernon Creek is close to the Sevier Desert, which has extremely low excess precipitation, unlike Muddy Creek. Thus, the consideration of other hydrological processes was necessary for Vernon Creek (ET in the lumped tank model) while the spatial variation of inputs is required for Muddy Creek. If ET is considered in the FDC method when computing I CP , the FDC method may perform better than the proposed approach.

Regional FDC method for regulated watersheds
It is impossible to evaluate the correlation between I CP and streamflow observation for regulated watersheds. With the low robustness of performance, using I CP as the only donor variable could result in a large bias in streamflow generation. Even in the case of transferring multiple I CP values, the bias would not be small due to the performance variability of the FDC method. Thus, the use of I CP was limited as one of the multiple donor variables. Neighboring streamflow observations were also transferred in order to make up the drawback of I CP . Hence, the role of I CP for regulated (or ungauged) watersheds was to capture the hydrologic responses not included in the neighboring streamflow observations. The simulated streamflows were higher than observed from April to October due to river diversions for agriculture at both regulated watersheds, except for year 2011 at Sevier River near Gunnison. Sevier River near Gunnison is located below the intersection between the Sevier River and the San Pitch River, but it was difficult to know the streamflow from the San Pitch River on a regular basis. Streamflow in the San Pitch River was negligible in dry and normal years due to the high agricultural water demand in the San Pitch River basin, but it could not be neglected in a wet year such as 2011. Thus the observed streamflows at Sevier River near Gunnison were greater than the simulated natural flows in a wet year as shown in Fig. 7b.
Conceptually, when the simulated streamflow is greater than the observed flow, the difference indicates the volume of diversions. However, a similar difference could be assumed to represent the volume of return flow from the agricultural areas when the observation is greater than the simulated value. As depicted in Fig. 7a, streamflow not decaying from November to March (the period of no diversions) demonstrated that the return flows through infiltration affected streamflow continuously. Return flows may affect streamflow during the period of diversions, but it was difficult to estimate the impact due to the complexity of combined flow. Simply, a positive difference between the simulated and observed flows in Fig. 7a indicated diversions including return flows, whereas a negative difference indicated return flow.
This study used observed diversions in the watersheds to validate the simulated natural streamflow. Most river diversions above Sevier River near Kingston were recorded for management purposes. Due to the high efficiency of water use in the agricultural area above this station, the effect of surface return flows may be small or negligible during the period of diversions. Even though the return flows through infiltration may affect streamflow, it was relatively small when compared to the total diversions and streamflow during the period of diversions. If one assumes that there is no significant return flows during the diversion season, the difference between simulated and observed flows could be considered to be the volume of diversions. Table 6 shows the sum of observed diversions in the main channel of the Sevier River above Sevier River near Kingston and the estimated volumes from the three methods. The actual volume of diversions would be a little greater than the observed because some diversions might not be observed in spite of the large coverage of the diversion monitoring in the watershed. Hence, although Table 6 shows that the regional FDC method provided a larger natural flow than the others, the estimated volume of diversions by the regional FDC method could be considered a possible prediction.
However, the volume difference between the regional FDC and the semi-distributed model in Table 6 ranged from 13 to 40 %. This relatively high variation may come from the low robustness of the FDC method, errors in the regional FDC, and uncertainty in the regionalized parameters of the conceptual models. With these error sources, the use of only one method may be inappropriate. It is apparent that the semi-distributed model provides the most trustworthy results due to its better performance. Shu and Ouarda (2012) recommended at least four streamflow observations as donor variables for good precision with the FDC methods. Thus, the regional FDC with two streamflows and I CP in this study could add more uncertainty than a case with more donor variables. An important goal of this work in using the regional approaches was to estimate the amount of water from streamflow without actual diversion data. In most of these situations data are limited, yet water managers require such information to better manage water demands. The results of this analysis, especially from Table 6, shows the regional FDC method could produce acceptable estimates with less time and effort than conceptual modeling. There are several limitations in the regional FDC method. For every regionalization approach, including the regional FDC method, adequate streamflow observations are necessary to have good estimates. Parajka et al. (2013) commented that studies with more than 20 gauging stations produced better and stable performance with deterministic models. The regional FDC method is also sensitive to the number of gauging stations. Although the density of gauging stations was low in this study, gauged watersheds in the regional analysis should be adequate in terms of the watershed scale and climatic characteristics to minimize bias. As mentioned earlier, multiple donor variables can also minimize errors caused by bias of a single donor set.

Conclusions
In this study, a conceptual snowmelt model, SNOW-17, using point snow observations, was extended using a modified FDC method to simulate streamflows in the semi-arid and mountainous Sevier River basin of Utah. The FDC method was later extended to simulate natural streamflows in regulated watersheds by incorporating a parametric regional FDC method. The FDC method could be a simple practical approach for streamflow generation for watersheds with limited data. The FDC method was compared with the lumped and semi-distributed tank models under similar data availability to simulate streamflows and later extended via regionalization to estimate natural flows in regulated watersheds.
The results show that the FDC method could be a practical option for snow-fed watersheds with high correlation between I CP and streamflow. Of course, the performance of the snowmelt model was a prerequisite for good performance. With streamflow observations, I CP could be correlated and can be a good donor variable without other neighboring streamflow observations. In spite of the simplicity of the FDC method, it could provide approximate estimates of natural flow in terms of water volume. The spatial variation of climatic variables in a watershed could determine the performance of the FDC method. High ET could result in low correlation between I CP and streamflow. Thus, the consideration of ET in the calculation of I CP can enhance the accuracy of the FDC method. As seen here, when I CP and streamflow are highly correlated, the FDC method is able to outperform the lumped and semi-distributed models. Without the burden of parameter optimization and related computations of hydrologic processes, the FDC method could generate approximate streamflows with comparable precision to conceptual modeling. Importantly, checking the correlation between I CP and streamflow would be a key step for good performance. In the case of regulated or ungauged watersheds, a regional FDC should replace the gauged FDC. In snow-fed watersheds of the study area, drainage area and elevation were important to characterize percentile flows. Soil properties such as proportion of clay, saturated hydraulic conductivity, and dry bulk density, were also significant variables for estimating percentile flows of the regional FDC. Streamflows simulated by the regional FDC produced acceptable streamflow estimates when compared to the other conceptual models. In this work, the simulated natural flow by regionalization was used to estimate the volume of river diversions in regulated watersheds with impaired streamflow observations. Both the regional FDC and regionalization of conceptual modeling estimated the approximate volumes of river diversions. Even though the regional FDC method produced more uncertain diversion volume, both estimation approaches could provide practical and acceptable values under data-limited conditions for water resources planning and management. In short, the FDC method can be a practical method for the simulation of natural flows in both gauged and ungauged or regulated watersheds, especially under limited data. However, the parameters of snowmelt modeling should be estimated using SWE observations as shown here. Other studies are necessary to determine the parameters of the snowmelt model for watersheds without SWE observations. Also, the difficulty of determining the recession coefficients for I CP calculation in ungauged watersheds is another remaining issue, since the typical values for gauged watersheds are assumed. In summary, the FDC approach used here could produce practical values of expected streamflows from point observations for watersheds with limited data.