Unit Hydrograph Modeling using Geomorphological Instantaneous Unit Hydrograph ( GIUH ) Method

Flood forecasting at Wonogiri Reservoir is restricted on the availability of hydrologic data due to limited monitoring gauges. This issue triggers study of unit hydrograph modeling using Geomorphological Instantaneous Unit Hydrograph (GIUH) which is based on Geographic Information System (GIS). Analysis of physical watershed parameters was conducted on Digital Elevation Model (DEM) data using software Watershed Modeling System (WMS) 10.1 and ArcGIS. Nash model and S-curve method were used to process triangular GIUH into hourly Instantaneous Unit Hydrograph (IUH) and Unit Hydrograph (UH) and then was compared with the observed UH of Collins method. A sensitivity analysis was conducted on parameter of RL and Nashmodel k. Evaluation of accuracy of the simulated GIUH runoff hydrograph was also conducted. The GIUH model generated UH with smaller peak discharge Qp, also slower and longer of tp and tb values than the observed UH. Accuracy test of the simulated GIUH runoff hydrograph using Nash-Sutcliffe Efficiency (NSE) shows that Keduang watershed gives a satisfying result, while Wiroko watershed gives less satisfactory result. The inaccuracies occur due to limited flood events used to derive the observed UH and stream tributaries that were not properly modeled based on Strahler method.


INTRODUCTION
One of the operational purposes of Wonogiri Reservoir is for flood control.However, flood forecasting based on unit hydrograph is still restricted on the availability of hydrologic data due to the limitation of monitoring gauges in watersheds.This issue triggers study of unit hydrograph modeling using Geomorphological Instantaneous Unit Hydrograph (GIUH) method based on the characteristics of physical watershed parameters.Recently, Geographic Information System (GIS) has been widely applied to estimate spatial parameters of watershed for hydrologic modeling purpose.
This study aimed to determine physical parameter and to derive unit hydrograph from the selected watersheds in the catchment area of Wonogiri Reservoir using GIUH method by applying GIS approach in determining the geomorphological characteristics of watershed parameters.This study examined the accuracy of the obtained unit hydrograph and flood hydrograph derived from the model.This study was conducted at the catchment area of Wonogiri Reservoir located in Wonogiri Regency, Central Java Province, Indonesia.The derivation of watershed geomorphological parameters was performed on Digital Elevation Model (DEM) data using ArcGIS and Watershed Modeling System (WMS) software.

Catchment Area of Wonogiri Reservoir
Wonogiri Reservoir has a total catchment area of about 1,343 km 2 and consists of 10 watersheds.Keduang watershed with an approximate area of 397.36 km 2 is the largest watershed, while Wiroko is the second one with an approximate area of 216.95 km 2 (Oktavia, 2013).Besides, there are four Automatic Water Level Recorder (AWLR) stations in the catchment area of Wonogiri Reservoir which are managed by Perum Jasa Tirta I installed in Keduang, Kali Tirtomoyo, Kali Temon, and Bengawan Solo Hulu Rivers.

Unit Hydrograph (UH)
Unit hydrograph is defined as direct runoff hydrograph at the outlet of watershed generated by 1 mm of effective rainfall occurring uniformly over the catchment area with constant intensity for a specific duration (Chow et al., 1988).

Instantaneous Unit Hydrograph
When the duration of effective rainfall is infinitesimal, the resulting hydrograph is an impulse response function namely Instantaneous Unit Hydrograph (IUH) (Rodriguez-Iturbe andValdes, 1979 in Chow et al., 1988).Response from the complete input of () is a direct runoff () which is stated in this convolution integral (Chow dkk., 1988) as shown in Equation (1). (1)

GIUH
GIUH is defined as a probability density function of a drop's travel time in a basin.This theory is introduced by Rodriguez-Iturbe and Valdes (1979) and then is enhanced by Gupta (1980) (Quan, 2006).Figure 1 illustrates the relation between hydrograph and topographic factors (Derbyshire, et al., 1981in Quan, 2006).
Rodriguez-Iturbe and Valdes (1979), in Quan (2006) assume an Instantaneous Unit Hydrograph (IUH) as a triangular that consists of peak discharge and time to peak which is formulated in Equation ( 2), (3), and (4) as follows.
= (3) with:   is peak discharge (hour -1 );   ,   are time to peak (hour) and base time (hour),   is length of the highest order stream (km),  is dynamic parameter velocity (m/s), and   ,  ,   are stream-area ratio, bifurcation ratio, stream-length ratio of Horton.
3 THEORETICAL FRAMEWORK

Stream Order
Based on Strahler classification method, the smallest recognizable channels with no tributaries are designated as first stream order (Chow et al., 1988).Stream order classification using Strahler method is shown in Figure 2 (Bras, 1990).

Horton's Ratio
Horton's ratios that consist of bifurcation ratio (RB), stream-length ratio (RL) and stream-area ratio (RA) are representative parameters of a given watershed and are fixed values for a given watershed system (Rai et al., 2009).Horton's ratios are obtained using Equation ( 5), (6), and (7).

Dynamic Parameter Velocity
For GIUH modeling, velocity value is required to represent the entire watershed.Dynamic parameter velocity (V) for a watershed can be estimated using combination of Kirpich formula and velocity relationship (Jotish et al., 2010) as shown in Equation ( 8), (9), and (10).

Nash Model
The Nash model (Nash, 1957in Rai et al., 2009) is one of the distributed rainfall-runoff model based on the concept of instantaneous inflow routing through a cascade of linear reservoir with equal coefficient storage.Karamouz et al. (2013) stated that relation between storage and discharge of each reservoir is assumed to be linear  = , where value of  is average delay time for each reservoir.
If there are  reservoirs for a given watershed, and then unit pulse of rainfall is inputted in a very short time ∆ → 0, resulted outflow is ordinate () of an IUH.Outflow resulted from the first reservoir is calculated with Equation (11): Outflow  1 () of the first reservoir flows into a second reservoir and results Equation ( 12): By continuing process in Equation ( 12), outflow for th reservoir is derived in the function of Gamma distribution as shown in Equation ( 13) and is known as Nash model.
while  and  are parameters of Nash model, where  is the number of linier reservoir, and  is the storage coefficient (hour).

Geomorphological Parameter Estimation of Nash Model based on GIUH
The complete shape of GIUH is obtained by linking qp and tp of GIUH with scale () and shape parameter () of Nash model.In Rai et al. ( 2009) Equation ( 14) and ( 15) are obtained by substituting and simplifying Equation ( 13),: Parameter of n is obtained by solving Equation (15) using Newton Raphson method.Parameter of k for a certain value of V is calculated using Equation ( 16).
3.6 Derivation of UH from IUH Derivation of UH from IUH is conducted using two methods.First is lagging method which sums two identical IUHs with a lagging time, tr in certain duration and identic IUH.UH is obtained by averaging the resulted ordinates.Second is S-curve method which sums some IUHs in sequence until fix discharge is obtained.The difference of similar S-curve of each time interval is the total sum of unit hydrographs during the time interval.Final UH is obtained by dividing the ordinates with time interval (US Army Corps of Engineers, 1994).

Statistic Method to Evaluate Model's Accuracy
a) Nash-Sutcliffe efficiency (NSE) NSE is calculated using Equation ( 17): with:   is the ith observation discharge value,   is the ith simulated discharge value,  ̅ is the mean of observed discharge data, and  is the number of observed data.
b) Relative Mean Error (RME) RME between peak discharge value of simulated hydrograph and observed hydrograph is calculated using Equation ( 18) and ( 19).
with:   is the percentage of relative error of each event,   is the peak discharge of observed runoff hydrograph, and   is the peak discharge of simulated runoff hydrograph.
c) Root of Mean Square Error (RMSE) RMSE of the peak discharge is obtained using Equation ( 20) and ( 21) as follows. (20) with:   is the relative error of each event,   is the peak discharge of observed runoff hydrograph, and   is the peak discharge of simulated runoff hydrograph.

Research Methodology
This research was done according to general flowchart as shown in Figure 3.

Accuracy Analysis on Resulted GIUH
For the verification of unit hydrograph and simulated direct runoff hydrograph, the accuracy indicators are peak discharge (Qp), base time (tb), and time to peak (tp) (Figure 4).Verification is performed using statistic methods of Nash-Sutcliffe Efficiency (NSE), Relative Mean Error (RME), and Root of Mean Square Error (RMSE).Considering flow accumulation process, minimum threshold value of 0.5 km 2 , 1 km 2 , 1.5 km 2 , and 2 km 2 were used to compare the identified stream area.In ideal condition, threshold values are optimized so that the first stream order tributaries are properly modeled based on criteria conforming the actual watershed condition according to Strahler classification method.
The next process was determination of AWLR stations as an outlet point to delineate the watershed.Results of the delineation process for Keduang is given in Figure 5.The area for Keduang and Wiroko watershed are 360.73km 2 and 183.92 km 2 respectively.While in Octavia (2013), the total area for Keduang is 364.043km 2 and Wiroko is 183.131km 2 .Total area of watershed was estimated according to the selected outlet point.Coordinate of AWLR stations as outlet points, in this case were different.
Figure 5 shows how minimum threshold value affects the number of identified stream tributaries.At the watershed, minimum threshold value of 0.5 km 2 gives greater number of stream tributaries than the threshold value of 2 km 2 .Less threshold value will yield a greater number of stream network.In order to determine appropriate threshold value, actual stream condition needs to be investigated through field observation.
Result of delineation should be verified first so that it fulfills Strahler's criteria.Meanwhile, in this study, delineation result with minimum threshold value of 0.5 km 2 was chosen because it gives better tributaries and best at representing the actual condition of the watershed.
Stream network classification was conducted following Strahler scheme order.The analysis results conclude that the highest stream-order for Keduang is 5, while for Wiroko is 4 as shown in Figure 6 and 7.

Analysis of Geomorphological Parameter of Watershed
Geomorphological parameters for Keduang and Wiroko watersheds are given in Table 1 and Table 2 The physical parameters are used to calculate Horton ratios (RA, RL, dan RB) estimated by semi-logaritmic regression curve.The values of physical parameter characteristics are listed in the following Table 3. Ordinates of Nash's IUH of mm/hour unit were converted into Nash's IUH of m 3 /s unit using watershed area as conversion factor.Then, IUH was derived into UH by lagging and S-curve method.

Comparison between Modeled Unit Hydrograph and Observed Unit Hydrograph
The observed unit hydrograph is an average of several selected events in each watershed.Several selected events and averaging method refers to previous study of Pradipta (2014) according to the most updated watershed area from the present study.The selected flood events for Keduang watershed are 8 cases, while in Wiroko are only 4 cases.Results of GIUH unit hydrograph modeling and average observed unit hydrograph of each watershed are shown in Figure 8 and Figure 9.The results show that the shape of the GIUH unit hydrograph has lower peak discharge and longer recession limb, while the observed unit hydrograph has steep rising limb, greater Qp, and relatively shorter tb.This occurs because the stream networks were not modeled properly according to the real condition in watershed, so the discharge of hydrograph becomes slower.Therefore, determination of minimum threshold value of watershed during flow accumulation modeling becomes important in order to gain stream networks that represents more accurately the actual watershed condition according to Strahler method.
Table 4 shows summary of the GIUH, the observed UH, and also the performance of calculation accuracy.
The result shows that the error percentage of each modeled UH by S-curve method is less than the error percentage of those by lagging method.The difference is affected by the variability of flood events used for the calculation of the observed unit hydrograph by Collins method.In Table 4, the accuracy calculation of Keduang watershed with 8 flood events gives a better result than of Wiroko watershed that only used 4 food events.

Sensitivity Analysis on RL Parameter
Sensitivity analysis was performed by setting RL parameter of Wiroko watershed at maximum normal value of 3.50.Wiroko watershed was used in the sensitivity analysis since the difference between observed discharge and the simulation results were significant with smaller area compared to the Keduang watershed.RL and K are length and reservoir factors, respectively which both define shape of the flow discharge hydrograph.Result of GIUH calculation is given in Figure 10 below.It shows value of RL = 3.50 creates greater peak discharge compared to the result of RL = 2.66.This is an accordance with Equation (2), it states RL value is proportional with value of qp.Besides, value of tp and tb become shorter, which is also an accordance with Equation (3) that states RL value is inversely proportional with value of tp.Moreover, value of RL used to determine k parameter of Nash model is proportional with the ordinate of resulted unit hydrograph.Meanwhile, value of RL is inversely proportional with resulted k value.It proves that value of RL affects the shape of the calculated unit hydrograph and it also proves that the accuracy of stream network modeling takes an important role in determining value of RL.
5.9 Sensitivity Analysis on Nash Parameter Sensitivity analysis was performed on k parameter of Nash model.In this analysis, k parameter was given in the value of 1, 2, 3, 4 and 5, while the other parameters were assumed constant.The analysis was applied on Wiroko watershed and the resulted unit hydrographs are presented as follows (Figure 11).In Figure 11, the smallest value of k = 1 results in unit hydrograph with the greatest peak discharge and the shortest tp and tb.The greatest value of k = 5 results unit hydrograph with the smallest peak discharge and the longest tp and tb.This is an accordance with Nash model that states value of k is inversely proportional with the ordinate of unit hydrograph.Moreover, as k also represents storage coefficient of each reservoir, the greater value of k lengthens time of the flow retained in reservoir and it makes the discharge is released more slowly.It is proved by the characteristic of hydrograph with a short peak and long slope, and vice versa.The sensitivity analysis shows that the GIUH unit hydrograph is sensitive to the value of k parameter.

Evaluation of Direct Runoff Hydrograph from GIUH
The unit hydrograph from GIUH modeling was chosen to simulate direct runoff hydrograph.The rainfall events used for simulation of direct runoff hydrograph were the same rainfall events used to calculate observed unit hydrograph.

Recommendations
Recommendations of this research are that verification to the stream network during the delineation process is important in order to obtain real condition of watershed, so that unit hydrograph will give a better accuracy.Besides, the GIUH model needs to be applied to all watershed in Wonogiri Reservoir catchment area, and it can then be used to forecast inflow flood hydrograph more accurately.Moreover, further study is necessary to observe the effect of DEM resolution and minimum threshold value of watershed area during flow accumulation process towards value of RA, RB, RL and parameter qp, tp, tb of GIUH.In addition, more flood events data to derive the observed unit hydrograph is crucial to represent a better watershed hydrologic condition.

Figure 4 .
Figure 4. Illustration of verification between unit hydrograph of GIUH and observed unit hydrograph.

Figure 8 .
Figure 8.Comparison between GIUH unit hydrograph and observed unit hydrograph of Keduang watershed.

Figure 9 .
Figure 9.Comparison between GIUH unit hydrograph and observed unit hydrograph of Wiroko watershed.

Figure 10 .
Figure 10.The effect of change on RL value on the modeled unit hydrograph in Wiroko watershed.

Figure 11 .
Figure 11.Results of sensitivity analysis on k parameter of Nash model of resulted GIUH unit hydrograph at Wiroko watershed.

Table 3
IUH was derived by calculating qp and tp of GIUH using scale parameter of (k) and shape (n) from Nash model.At Keduang watershed, value of n is 3.034 and k is 2.586.Meanwhile, at Wiroko, the value of n is 3.085 and k is 3.835.Then, ordinates of Nash's IUH at time t can be calculated by using values of n and k.

Table 4 .
Comparison between GIUH and observed UH Evaluation of the direct runoff hydrograph was conducted using statistic criteria of NSE, RME and RMSE.The application of GIS in determining watershed's physical parameter characteristic is able to derive unit hydrograph of GIUH method with limited hydrologic data or unavailability of rainfall-runoff data.The unit hydrograph modeling by GIUH and Nash model approach conducted on Keduang and Wiroko watersheds found smaller Qp, later tp and longer tb than by the observed unit hydrograph.The accuracy analysis shows the rainfall-runoff simulation in Keduang watershed gives good and satisfying results, while in Wiroko watershed the results are less good and satisfying.Results of the sensitivity analysis show that RL parameter and k parameter of Nash model affect the shape of GIUH unit hydrograph.The inaccuracies because of the limited flood events used to derive observed unit hydrograph and of the stream tributaries that are not properly modeled because there is no verification process to calibrate the model with the actual stream based on criteria conforming the actual watershed condition according to Strahler method.