Hydraulic and Hydrologic Modeling of Steep Channel of Putih River, Magelang District, Central Java Province, Indonesia

Hydrologic and hydraulic modeling are important to be conducted to examine the watershed response based on a rainfall input, especially over disaster-prone watershed such as Putih River watershed in Magelang, Central Java Province. A GIS-based gridbased distributed rainfall-runoff model was used to simulate the rainfall-runoff transformation. A two-dimensional hydrodynamic flow modeling was then carried out to simulate the flood processes on the stream and floodplain area. A sensitivity analysis was conducted on infiltration rate, Manning’s n value, and rainfall intensity. Infiltration rate, Manning’s n value, and rainfall intensity give considerable effects to the resulted flow hydrographs. The modeling results show that the results of hydrologic-hydraulic modeling is in good agreement with the observed results.


INTRODUCTION
The development of Geographic Information System (GIS) has a significant impact on hydrologic and hydraulic modeling.As GIS-based hydrologic and hydraulic data are now widely available, many recent hydrologic and hydraulic models have been developed based on GIS.Lahar flood occur regularly in Putih River, Central Java Province.Since the upstream is located near the peak of active volcano Merapi, Putih River has a huge amount of lahar producing sediment materials on its upstream.Located on the steep upstream, the sediment materials, when saturated by rainfall, have a high probability to produce destructive lahar flood threatening the populated downstream area.This research is aimed to assess the lahar flood in Putih River by applying a whole modeling process ranging from hydrologic modeling to hydraulic modeling.

HYDROLOGY AND HYDRAULIC MODEL
Zhang and Savenije (2005) simulated saturation overland flow by relating the variable source area to both the topography and groundwater level.Koren, et al. (2005) evaluated the performance of a grid-based distributed hydrology model when transforming rainfall to stream flow which resulted in runoff volumes that was in good agreement with the observed data.Liu, et al. (2009) modeled runoff volume and peak discharge based on several flood events using grid-based distributed kinematic hydrological model with rainfall data derived from climate radar.The research was able to reproduce peak discharge conforming to the field observed data, and also show that the model can be used to evaluate the effect of change in land cover and soil on the stream flow due to human activity.Miyata, et al. (2014) examined the effects of intensive local rainfall on flashflood events using a grid based distributed hydrological model.The hydrology model is able to both evaluate the basin response due to rainfall, in form of flow, and to identify contributions of each sub-basins to the quickly rising stream discharge.

Watershed
A watershed is an area of land draining into a stream at a given location.One of the watershed morphology that reflects on flow discharge is watershed hydrologic length and slope.Hydrologic length is the length of the main stream from the outlet to the watershed boundary.Hydrologic length reflects on time behavior of discharge.Watershed slope is the total elevation change of the main stream divided by hydrologic length.The watershed slope is an important factor in that affects runoff momentum.The longitudinal profile of the river, which is extracted from topographic maps, can be used to divide stream reaches into slope categorization that reflects profile morphology as listed in Table 1 (Rosgen, 1994).

Design Hyetograph
Hyetograph is a plot of incremental rainfall intensity against time interval.There are several methods to develop hypothetical rainfall distribution, one of which is Alternating Block Method (Chow, et al., 1988).The general shape of design hyetograph by ABM is shown in Figure 1.In Japan, the equation to calculate rainfall intensity for any duration based on daily rainfall depth is as follow (Mori, et al., 2003): where i t is rainfall intensity for any duration t (mm/hour), R 24 is daily rainfall depth (mm), and t is rainfall duration (hour).

Hydrograph
A streamflow or discharge hydrograph is a graph or table showing the flow rate as a function of time at a given location on the stream.Figure 2 shows components of a streamflow hydrograph during a storm (Chow, et al., 1988).The rising limb is the limb that has a positive gradient indicating the rising of discharge (point A to B).The recession limb is the limb that has a negative gradient indicating the falling of discharge (point B to C). Point C to D represents baseflow recession.Lag time between peak rainfall and peak discharge indicates the time required for a surface runoff flow to reach an outlet point where the hydrograph is observed.The shape of the hydrograph describes the response of a basin toward a rainfall input which is different to another depending on several climatic, topographic, and geological factors.Precipitation contributes to various storage and flow processes.The precipitation which becomes streamflow may reach the stream by overland flow, subsurface flow, or both (Chow, et al., 1988).
The passage of overland flow into a channel can be viewed as a lateral flow.A discharge of the overland flow, q 0 per unit width, passes into the channel, with the length of the channel, L c , so the discharge in the channel is Q = q 0 L c .To find the depth and velocity at various points along the channel, an iterative solution of Manning's equation is necessary.
where u is flow velocity (m/s), n is Manning's roughness coefficient, R is hydraulic radius (m), and S is energy line gradient.
Infiltration is the process of water penetrating from the ground surface into the soil.Many factors influence the infiltration rate, including the condition of the soil surface and its vegetative cover, the properties of the soil, such as porosity and hydraulic conductivity, and the current moisture content of the soil (Chow, et al., 1988).
Hydrologic-hydraulic model is used to model rainfallrunoff transformation.The output of the model is flow hydrograph at a particular basin outlet.The model is a grid-based distributed model and consists of two domains.The first domain is hillslope and the second one is stream.The calculation in hillslope domain includes water balance calculation and flow routing on hillslope, while the calculation in stream domain consists of channel flow routing.The scheme of hydrology model in hillslope domain is shown by Figure 2. In each grid, there are 3 layers on which runoff flow, sub-surface flow, and vertical flow between layers are calculated.The relation of a grid with its 8 surrounding grids is defined by flow direction.Flow direction of a grid to adjacent grid is determined by a grid whose elevation level is the lowest among those 8 surrounding grids.
Figure 3.The illustration of hydrologic model scheme on hillslope domain (Miyata, et al., 2014) The illustration (Figure 3) explains the process of rainfall dropping on a surface of a hillslope grid which infiltrates to the layer A. The saturated flow of layer A then percolates to the layer B. The deeper percolation of flow on layer B to the deeper layer could occur, yet percolation flow to the deeper layer is assumed not flowing back to the slope or stream.
The hydrologic model calculates both runoff flow and lateral flow based on Manning equation and Darcy's law.The runoff flow calculation takes into account the water balance from rainfall, flow from adjacent upstream grid, and flow to adjacent downstream grid, infiltration flow to the lower layer, and evaporation.The relation between depth of surface runoff flow, h OF , and runoff discharge per unit width, q OF , is represented in the following equations: where r is intensity of effective rainfall (rainfallinfiltration), i OF is the energy line gradient (approached by surface slope gradient), and n M is Manning's roughness coefficient (Miyata, et al., 2014).
Sub-surface lateral flows occur between layer A and layer B. The equation to calculate lateral discharge per unit width, q A , and the flow depth, h A , in layer A, are: where q in_A is infiltration flow from surface to the layer A, q in_B is infiltration flow from layer A to layer B, i A is hydraulic gradient of flow in layer A, and K A is permeability coefficient of layer A. The calculation of lateral flow discharge and flow depth in layer B uses the equation same as of the layer A (Miyata, et al., 2014).
Furthermore, the calculation of discharge on stream domain, q CH , uses the same kinematic wave model which is being used in flow model in hillslope domain.The composition of stream grid is carried out as a series of cross sections with lateral coming from adjacent right and left grids.The following equation is used to calculate the flow in stream domain: where r is rainfall on a stream grid, L is lateral flows from left and right adjacent grid, i CH is hydraulic gradient of stream flow, and n M is Manning's roughness coefficient.
3.5 Two-Dimensional Hydrodynamic Model SIMLAR v1.0 SIMLAR is developed by Yogyakarta Sabo Centre in cooperation with Universitas Gadjah Mada (Hardjosuwarno, et al., 2012).The hydrodynamic model is used to model flow along a particular stream or channel.The flow which is modeled using SIMLAR is of debris flow type.The debris flow is calculated as a fluid unity of debris flow and suspended sediments.Some equations used in the model are described as follows.
To calculate the permanent flow velocity at a uniform channel cross section, Manning equation is used.
Based on flow velocity from the Manning equation, flow discharge at a particular channel cross section is calculated using the following equation.
where Q is flow discharge (m 3 /s), and A is flow area of a channel cross section (m 2 ).
For a non-permanent and non-uniform flow, the debris equations are developed considering the mass conservation principles and force and momentum conservation principles.
Mass conservation equation for two-dimensional flow is: Force and momentum equation for two-dimensional flow is: where H is depth of flow (m), M is flow discharge per unit width at x direction (m 2 /s), N is flow discharge per unit width at y direction (m 2 /s), g is gravity acceleration (m 2 /s), and τ is shear stress.

Location of Study
The research was carried out on Putih River watershed in Magelang District, Central Java Province.It has upstream end on the peak of Merapi Volcano.It flows downward the Merapi slope to the southwest direction until it intercepts with Blongkeng River on the downstream end.The boundary of Putih River watershed is shown by Figure 4. c) Physical simulation parameters Several physical parameters for hydrologic-hydraulic simulation such as Manning coefficient and properties of bed material are collected from previous researches, namely Widowati (2015), and Musthofa (2014).

Research Tools
Several basic data processing applications such as Microsoft Word, Microsoft Excel, Notepad++ are used to process the data and conduct the whole simulation.The GIS data preparation is done using ArcMap 10.1, QGIS and GRASSGIS application.The hydrology simulation of transforming rainfall-runoff is carried out using the grid-based distributed hydrology model (Miyata, 2014).The flood hydrodynamic simulation is conducted using 2 dimensional hydrodynamic model SIMLAR (Hardjosuwarno, et al., 2012).

Delineation of Putih River watershed
The result of the watershed delineation gives the total area of Putih River watershed (Figure 5) as 23.011 km 2 with a total length of the main river of 21.639 km.The sub-watershed has an area of 8.432 km 2 .The sub watershed of Putih River has riverbed slope ranging from 3.0° to 7.6°.

Rainfall -Runoff Modeling
The rainfall-runoff modeling is conducted using kinematic distributed model.The final results are flow hydrograph on the outflow point.A sensitivity analysis was conducted on infiltration rate, Manning's n value, and rainfall intensity.The rainfall event used during the simulation of sensitivity analysis is the rainfall event on December 8 th 2010 recorded at Ngepos station (Figure 6).
The simulated watershed boundary in Putih River covers Putih River upstream sub-basin with an outlet point located at approximately upstream of PUD-02.
The area of Putih River inflow sub-basin is 8.43 km 2 .Figure 7 presents the scheme of slope and stream domain of the simulation.

Sensitivity analysis on infiltration rate
The sensistivity analysis used a constant slope Manning's n = 0.7 and a constant channel Manning's n = 0.03. Figure 8 shows the results of sensitivity analysis on infiltration rate in Putih River upstream watershed.Table 2 provides the values of the resulted peak discharge and peak time.High infiltration rate enables more surface f low to infiltrate into the soil and thus result in low surface runoff discharges.High infiltration rate also diminishes the discharge of flows from far-grids that enter the stream at a later time which is indicated by later peaks in the hydrographs of small infiltration rate.The later rising parts of the hydrographs of small infiltration rate (f = 0.6, 1.2, and 2.2 mm/hour) appear due to the flow from slope-grids located in farther distance from the stream that just enter the stream at a later time.In the hydrograph with high infiltration rate (f = 5.0 mm/hour), this additional rising part of hydrograph is not observed because the flow is already being infiltrated into the soil along its path to the stream.Small infiltration rate also gives flow a chance to accumulate faster, thus hydrograph of small infiltration rate rises earlier.

Sensitivity analysis on Manning's slope n values
Referring to the Manning's n value from literature study (Engman, 1986), the surface condition in Putih River is considered ranging from grass to woods.
Thus, it was taken the Manning's n values ranging from 0.3 to 0.8 for sensitivity analysis of Manning's n value.The result of sensitivity analysis of Manning's n values is shown by Figure 9. Table 3 provides the values of the resulted peak discharge and peak time.For a given rainfall intensity and infiltration coefficient, smaller Manning's slope n values results in higher flow velocity than higher n values.High flow velocity also gives effect in reducing the chance of the flow to be infiltrated into the soil.Thus, smaller n values produce greater and faster flow discharge entering the stream than higher n values.
In addition, small n value which results in high flow velocity is correlated with a faster draining time.Flow with high velocity is easier to be drained out of slopes than slow flow, thus makes it quicker for the flow both to form the peak discharge and to decline.So, small n value makes the hydrographs sharper: shorter duration to reach the peak discharge, higher peak discharge rate, and steeper recession limb.

Modeling with variation of rainfall intensity
The design hyetograph was derived using a modified Mononobe equation and distributed using Alternating Block Method (ABM) technique.The results of modeling with variation of rainfall intensity in Putih River is shown in Figure 10.Although the infiltration rate and Manning's slope n value are set constant of f = 5 mm/hour and n = 0.7 respectively, each rainfall produces different shape of hydrograph.
The hypothetical rainfall of P 2-hours = 10 mm produces a very low hydrograph with maximum peak discharge of 0.7 m 3 /s at t = 215 minute.The hypothetical rainfall P 2-hours = 15 mm produces a terraced hydrograph with a discharge rising at t = 180 min and a peak discharge rate of 3,541 m 3 /s at t = 285 minute.The hypothetical rainfall P 2-hours = 20 mm starts to rise at t = 155 minutes and reaches maximum discharge rate 9,277 m 3 /s at t = 245 minutes.The hypothetical rainfall P 2- hours = 25 mm starts to rise at t = 145 minutes and reaches maximum discharge rate 15.99 m 3 /s t = 225 minute.
The hydrographs have different shape one another.The hydrographs of higher rainfall depth rise earlier have higher discharge.The maximum discharge also occurs earlier in the hydrographs of higher rainfall.By the specified rainfall distribution, the maximum discharge occurs after the rain stops.The rainfall intensity affects the flow discharge greatly, and produces hydrographs of rainfall intensity distinctly.

Hydrodynamic Modeling
The hydrodynamic debris flow modeling was carried out on Putih River basin.The inflow point is assigned at the downstream of PUD-02.The downstream boundary of the modeling is the upstream of PU-C0 Sukowati.The inflow hydrograph (Figure 11) was obtained from rainfall-runoff modeling in Putih River upstream sub-basin of January 9 th 2011 rainfall event from Ngepos ARR station, using infiltration rate of f = 5 mm/hour and Manning's slope n value of n = 0.7.The hydrodynamic modeling of debris flow in Putih River was conducted for 210 minutes.The result of hydrodynamic lahar modeling is a spread of water and sediment inundation on the simulated river and riverbank.Figure 12  Four cross-sectional stream profiles (CS 1, CS 2, CS 3, and CS 4) at the Putih River reach which intersects Magelang public road are extracted to observe the effect of sediment erosion and inundation near the vulnerable public road.Figure 14 shows the riverbed changes due to sediment erosion and sediment inundation at the specified cross-sectional stream profiles.The cross-sectional stream profiles are extracted at t = 180 minutes and t = 210 minutes only based on the time the sediment flow has reached the Magelang road.
At CS 1, it is observed that an erosion occurs with the maximum erosion depth of 0.73 m at t = 210 minutes.The riverbed has been lowered due to the erosion.At CS 2, it is observed that a sediment inundation occurs with the maximum sediment inundation depth of 0.20 m at t = 210 minutes.At CS 3, it is observed that a sediment inundation occurs with the maximum sediment inundation depth of 0.49 m at t = 210 minutes.At CS 4, it is observed that a sediment inundation occurs with the maximum sediment inundation depth of 0.94 m at t = 210 minutes.
According to the field event record of lahar flow in Putih River at January 19 th 2011 (Delson, 2012;SID Jumoyo, 2015), the lahar flow sediment material inundated the Magelang road with the sediment inundation depth ranging from 2 to 3 meters.Compared to the simulation results, the highest simulated sediment inundation depth around Magelang road is 0.94 m at t = 210 minutes.During the event of January 19 th 2011 lahar flow, the lahar flow was also reported to overtop the road and flow into the old stream course.The simulated results also show an indication of the water and sediment flow overtopping the road and starting to flow into the old stream course.However, due to the simulation duration was limited to 210 minutes only, the farther propagation of water and sediment flow could not be observed yet.

Recommendations
Of all the research results, there are several recommendation for a better future research: a) The rainfall-runoff model should be applied on basins with different topography and stream structure, such as basins with finer slope and more stream tributaries.b) The radar-derived rainfall data should be applied as rainfall input in order to produce a spatially better result.c) A hydrodynamic modeling of lahar flow using DEM with smaller resolution should be conducted to have a more accurate result.d) A longer duration of hydrodynamic modeling of lahar flow should be conducted to observe the longer debris flow propagation over a stream reach.

Figure 2 .
Figure 2. The flow hydrograph at a basin outlet during a rainfall-runoff simulation by a particular observed rainfall event

Figure 4 .
Figure 4. Location of study area in Putih River watershed

Figure 5 .
Figure 5.The delineation of Putih River watershed and subwatershed.

Figure 6 .
Figure 6.Rainfall event of December 8 th 2010 in sensitivity analysis in Putih River.

Figure 7 .
Figure 7.The scheme of stream and watershed used in rainfall-runoff modeling in Putih River.

Figure 8 .
Figure 8. Simulated channel flow discharge hydrographs at outlet point using different values of infiltration rate in Putih River.

Figure 9 .Figure 10 .
Figure 9. Simulated channel flow discharge hydrographs at outlet point by different Manning's slope n values in Putih River.

Figure 12 .
Figure 12.Simulated sediment spread of debris flow in Putih River at t = 210 minutes
rate affects the magnitude of flow discharge.The higher infiltration rate, the lower the resulted runoff flow discharge.Infiltration rate also affects the shape and peaks hydrograph b) Manning's n value affects the flow velocity.The small Manning's n value results in fastaccumulated flow and high discharge.c) Rainfall intensity dominantly affects the typical hydrograph, yet, the effect of rainfall intensity toward the hydrograph parameter is different from one rainfall intensity and other intensities.

Table 2 .
Peak discharge (Q) and and peak time (t p ) of simulated hydrographs in sensitivity analysis of infiltration rate (f) in Putih River.