160 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 A New Algorithm For The Grid Cell-Based Runoff Routing Model Based on Travel Time Concept Baina Afkril1*, M. Pramono Hadi2 and Slamet Suprayogi2 1Doctoral Program of Geographical Science, Faculty of Geography, Gadjah Mada University, Jl. Kaliurang, Sekip Utara, Bulaksumur, Yogyakarta, 55281, Indonesia 2Faculty of Geography, Gadjah Mada University, Jl. Kaliurang, Sekip Utara, Bulaksumur, Yogyakarta, 55281, Indonesia *Corresponding Author : bafkril@gmail.com Received 8 April 2020/ Revised 5 May 2020 / Accepted 20 May 2020/ Available Online 1 June 2020 Abstract The grid cell-based routing model has recently been used to simulate direct runoff hydrographs at catchment scales. This study develops a flexible event-based runoff routing algorithm to simulate a direct runoff hydrograph (DRH). The experiment was based on the spatiotemporal inputs of a hydrological data set. The flexibility is based on the time step and grid cell size applied in the original STORE-DHM. Rainfall distribution was obtained using radar data adjusted by the measured point ground, while the runoff yield was determined using the NRCS-CN method. The parameter distribution was captured in the GIS environment as raster data formats. Furthermore, it was converted into ASCII data formats for scripting the routing algorithm using Matlab programming codes. The model algorithm was tested for storm events within two small study river systems in Yogyakarta, Indonesia. One event in each catchment was selected and calibrated to the observed hydrograph, treating the Curve Number (CN) and Manning coefficient (n) values as parameter calibrations. In the end, two events were selected for validation. The proposed routing model algorithm simulates DRHs of all selected events in the study areas with excellent performance. The Nash- Sutcliffe coefficient was greater than 0.75 for all DRH during validation, and the volume bias and peak discharge error were less than 25%. Keywords: Algorithm; Cell-based runoff routing; Travel time; GIS; Direct runoff hydrograph. 1. Introduction Rainfall-runoff models at catchment scales are simplified hydrological processes and mechanisms. However, capturing runoff phenomena in catchments is still complex. These models are constructed based on mathematical descriptions of the hydrologic cycle. Their architectures are determined according to their purposes (Singh & Woolhiser, 2002). GEOSFERA INDONESIA p-ISSN 2598-9723, e-ISSN 2614-8528 Vol.5 No. 2 (2020), 160-185, August, 2020 https://jurnal.unej.ac.id/index.php/GEOSI DOI : 10.19184/geosi.v5i2.17351 Accredited by the Ministry of Research , Technology , and Higher Education of the Republic of Indonesia, No. 30/E/KPT/2019. https://jurnal.unej.ac.id/index.php/GEOSI https://drive.google.com/file/d/1rSnVAs6cuHOwHl5BJ87CL2L6K5DQz7s6/view 161 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Runoff generation and routing are the essential components that should be prepared in every modeling (Beven, 2012). The runoff generation counts how much rainwater turns into overspill and becomes part of a storm hydrograph. It is related to rainfall and catchment characteristics, which are significant issues in hydrological abstraction. Moreover, it shows how the catchment responds to rainfall by producing excess rainwater flowing to the surface and subsurface. The routing component represents the runoff distribution to shape the hydrograph at the outlet or any selected point of observations along the channel pathway in the catchment. The runoff yield can be counted using the curve number (CN) method, also called the Soil Conservation Service-Curve Number (SCS-CN) or Natural Resources Conservation Service-Curve Number (NRCS-CN) method (USDA, 2004a). To estimate the runoff, the approach relies on land covers and treatments, soil types, and antecedent hydrologic conditions. It is a conceptual rainfall loss calculation method supported by empirical exercises, which also represents the infiltration loss model (Ponce & Hawkins, 1996). Although the CN method was initially developed from experimental agricultural watersheds in the USA, it has been adapted by scientists worldwide in counting the runoff yield. The CN method does not involve spatial variability in counting runoff (Ponce & Hawkins, 1996). Soulis & Valiantzas (2012) introduced a simplified concept of a two-CN heterogeneous system to figure out its spatial variability in a watershed. The study established that the approach sufficiently describes the CN-rainfall variation in natural watersheds. Gonzalez et al. (2015) proposed a vegetation correction factor to adjust the vegetation- adjusted CNs and applied it for monthly runoff estimation. They recorded better results compared to standard approaches. Also, they argued that the adjustment is vital for flash flood monitoring and forecasting. It is quite challenging to represent the CN-rainfall variation in a real heterogeneous catchment. In hydrological modeling, calculating the composite CN using the traditional SCS-CN is tedious and time-consuming (Rawat & Singh, 2017). The use of GIS technology has made it easier to spatially present the values and calculate the runoff yields within the catchment. Bansode & Patil (2014) estimated the runoff in GP-3 watershed, India, using SCS curve number and ArcGIS on a daily, monthly, and yearly basis. The results showed linear correlations between rainfall and runoff, where the yearly runoff was the best-fitted correlation. Using the same method and approach, Rawat & Singh (2017) estimated the surface runoff from Jhagrabaria’s semi-arid ungauged agricultural watershed and established a strong linear correlation of annual rainfall-runoff. Additionally, Ahmad et al., (2015) 162 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 applied the curve number method with remote sensing and GIS to estimate runoff potential within the Sheonath river basin in India. They revealed that the remote sensing and GIS- based SCS-CN are essential in estimating runoff within catchments of similar geo- hydrological characteristics, as well as in land use planning and watershed management. Vojtek & Vojteková (2016) applied the curve number method, coupled with GIS, in estimating surface runoff to define potential flood risk areas in the Vycoma catchment. The study concluded that the approach is suitable for locating potential risk areas prone to flooding. Maina & Raude (2016) assessed suitability for harvesting surface runoff. They used the curve number method with geospatial techniques in Njoro Catchment, Kenya. The results showed that about 50% of the catchment had curve numbers between 82 to 89, indicating the potential for rainwater harvesting. Satheeshkumar et al., (2017) used the SCS-CN method and GIS approach to estimate rainfall-runoff in Pappiredipatti watershed, South India. The approach was proven efficient, consuming less time and facility to handle extensive data set. Site selection of artificial recharge structures could be identified by the watershed as a larger environmental area. Rohman et al., (2019) examined the impact of land-use changes on curve numbers to quantify the effectiveness of the Natural Flood Management (NFM) approach in the Ciliwung Basin, West Java, Indonesia. The results showed that flood risks are primaril y affected by the changes in the CN values. Al-Juaidi (2018) evaluated the impact of land-use alterations on runoff volumes for the Gaza Strip using a simplified GIS-based SCS-CN method. The results showed that land-use changes play a significant role in CN number and runoff volumes. In this paper, the CN method was used as predictive values based on hydrological catchments. The variability of the runoff yields in a catchment was implied from the constructed grid cells. The routing structures and procedures can be simple or complex depending on the spatiotemporal considerations, physical processes involved, and computation resources. The development of Geographic Information System (GIS) technology and computing program codes have eased all the routing procedures, leading to the fast production of the model outputs. GIS technology can capture all the catchment boundaries. The spatial and temporal variability of any hydrological characteristics within the catchment can be depicted in a more detailed manner. The area fraction of the catchment in grid cells makes its application advanced in the distributed hydrological modeling. Furthermore, the routing program codes can be incorporated within the GIS environment or built separately and then interconnected using command tools available in the GIS software. 163 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 The grid-based runoff routing based on the travel time concept has recently been proposed to simulate direct runoff hydrograph (DRHs) at catchment scales. For instance, Melesse & Graham (2004) proposed a grid-based runoff routing model based on time- invariant. They assumed that the travel time in each cell is not varying during the storm event. The proposed method constructed the runoff hydrograph directly through the cell to cell routing. Du et al., (2009) studied the variability of travel time within the grid cells. They modified the method proposed by Melesse & Graham (2004) by involving time variability in runoff generation and spatially distributed direct hydrograph (SDDH) model. The travel time from the origin of the runoff cell to the outlet cell was determined cumulatively along the flow path. The direct runoff hydrograph (DRH) at the outlet was obtained by summing up the volumetric rate from total contributing cells for all time intervals. Zhao & Wu (2015) applied the concept of grid cell travel time on the soil surface with different micro-topography scales to simulate runoff hydrographs. The contributions of three parameters, including rainfall intensity, mean flow velocity, and ponding time depression were used to determine the flow time. Additionally, the duration from the most upstream to the outlet cell was defined as a sum of all travel times along the path. The runoff rate was estimated by the summation of the rain rate from all contributing cells for all time intervals. Asfaw et al. (2018) developed a simple runoff generation and routing in the GIS environment to construct an event-based model prediction of metaldehyde concentrations at certain abstraction sites. Using the curve number method in a surface runoff generation module, the model predicted the arrival of peak metaldehyde concentrations using the curve number method in surface runoff. To determine the total travel time from the source cell of runoff to the outlet, it is acceptable to physically sum up all cell travel periods along the flow path. However, Kang & Merwade (2011) established that counting for a volumetric flow rate in a cell using the SDDH approach is inconsistent with the mass balance principle in the cell. This is because the volumetric flow rate in a cell is counted several times from upstream cells of different paths. Consequently, the outlet cell consistently receives high volumetric flow rates due to the repeated computation of flow accumulation. To overcome these issues in the travel-time concept, Kang & Merwade (2011) proposed the storage-release based distributed hydrologic model (STORE DHM). The approach treats all cells in the raster grid as storages for water from adjacent cells. The stored water is then released to downstream cells using a continuity equation combined with Manning's formula. The model has the capability to maintain the water balance in each cell since the processes of incoming-storing-outgoing in all cells are 164 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 counted once in each time step. However, the approach has flexibility issues since the minimum travel time must be determined first using the critical cell travel time (CCT) condition. This ensures that the travel time in a cell always exceeds the chosen fractional time of rainfall. Using this approach means the cell size should be selected to satisfy the cell’s travel time for the chosen rainfall time fraction. This is especially the case in small catchments where large grid cells are inappropriate. In this study, a new algorithm for runoff routing based on the concept of the STORE DHM without relying on the CCT is presented. The process of storage-release involves using different approaches when the travel time is less than the model time step before proceeding to the next step. This is conducted by applying a looping mechanism as an iterative process in the model. A “loop” statement in a computation program is a control statement that executes a block of codes repeatedly to meet a condition given for the block. 2. Methods 2.1 Area of Study The model program was applied to two small study river systems in Yogyakarta Province, Indonesia. The study river catchments were delineated using a 15 m DEM that was converted from a 12.5 m contours map. The first river system is located at the upstream of the Sumur Mbandung catchment, as shown in Figure 1. The second river system is located at the upstream of the Sempor catchment, as shown in Figure 2. The study area at Sumur Mbandung catchment covers an area of about 1.84 km2. It is located at elevation ranges of 160 m – 390 m above sea level. The outlet is at 439,538.77 meters East (mE) and 9,126,419.00 meters South (mS) of the Universal Transfer Mercator (UTM) coordinate system of zone 49S. The second study area covers an area of about 1.47 km2 in a river system at the upstream of the Sempor catchment. It is located at the elevation ranges of 600 m – 831 m above sea level while the outlet is at 432,023.00 mE and 9,155,540.00 mS of the UTM coordinate system of zone 49S. 165 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Figure 1. Study river system in Sumur Mbandung catchment 166 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Figure 2. The study river system in Sempor catchment 2.2 Data Resources The model needs rainfall data, land use and land cover (LULC) types, along with soil types to be processed later as input data. Also, flow data is used for model testing. The rainfall data were collected from two sources, including the portable automatic rainfall recorder rain gage (ARR) and the meteorological radar. One ARR of 0.2 mm for 1 tick measurement was installed within each of the study areas, as shown in Figure 1 and Figure 2. The rainfall radar data was from the X-Band Multiparameter Radar (XMPR) operated by the 167 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Hydraulic Laboratory of the Engineering Faculty, Gadjah Mada University. The radar is installed at the Merapi Volcano Museum, about 3 km from Sempor Catchment and about 30 km from Sumur Mbandung catchment. The rainfall radar data is available at http://data.hydraulic.lab.cee-ugm.ac.id/. After brief verifications in the field, the LULC types were delineated from the Google Earth Map (Image of August 2016). Figure 3 and Figure 4 show the LULC maps of the study areas in Sumur Mbandung and Sempor catchments, respectively. Figure 3. Map of land uses and land covers in the study area of Sumur Mbandung catchment http://data.hydraulic.lab.cee-ugm.ac.id/ 168 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Figure 4. Map of land uses and land covers in the study area of Sempor catchment The forest predominantly covers the study area in Sumur Mbandung catchment. In the middle part, there is a wide area of non-irrigated paddy fields with high runoff production. In Sempor catchment, the study area is predominantly covered by agricultural land. The Zalacca palm plantation covers the most significant area among all agricultural lands. 169 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 The soil types were obtained from the map provided by the Bureau of Planning and Development of Yogyakarta. Based on the soil map, the study area in Sumur Mbandungcatchment has latosol soil with a sandy, loamy clay texture. In contrast, Sempor catchment is mostly covered by regosol soil with a sandy texture. The flow data were obtained by installing two automatic water level readings (AWLRs) at the outlet of the river system. The readings were installed both at upstream and downstream sections. The flow was determined using a continuous slope area method by applying Manning’s formula of two cross-sectional areas. The DRHs were obtained using the constant slope method of baseflow separation. 2.3 Model Construction 2.3.1 Water Balance Concept In the grid-based water balance concept, a cell is perceived to be a bucket that stores, receives, and releases water. Representation of mass flow balance in the grid cells routing can be written as follows: (Kang & Merwade, 2011) 𝑆(𝑡)𝑗,𝑖 = 𝑄(𝑡)𝑗,𝑖 ∆𝑡 + 𝑆(𝑡)𝑗,𝑖−1 + ∑ 𝐼(𝑡)𝑘,𝑖∆𝑡 8 𝑘=1 − 𝑂(𝑡)𝑗,𝑖 ∆𝑡 (1) where S(t) is water stored in a cell (m3), Q(t) is the runoff in a cell resulted from excess rainfall within Δt, or excess rainfall intensity (m3/s). I(t) is the inflow from the adjacent cells (m3/s), O(t) is the outflow to the downstream cell (m3/s), Δt is the time step interval (s), i is the incremental time step (i = 1, 2,…, n), j the cell’s index (j = 1, 2,…, m), while k is the adjacent cell’s index (k = 1, 2,…, 8). The outflow, O(t), in eq. (1) depends on the residence time or flow travel time in each cell. The following conditions are applied to calculate the outflow 𝑂(𝑡)𝑗,𝑖 = { 𝑆𝑗,𝑖−1 ∆𝑡 , if 𝑇𝑗,𝑖−1 < ∆𝑡 𝑆𝑗,𝑖−1 𝑇𝑗,𝑖−1 , if 𝑇𝑗,𝑖−1 ≥ ∆𝑡 (2) The conditions given in eq. (2) indicate that all stored water is released to a downstream cell when the travel time, Tj,i-1, is less than the time step (Δt). Furthermore, only a portion of or exactly all of the stored water is released when Tj,i-1 is greater than or equal to Δt. 170 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 In applying the routing concept of eq. (1), the cell networking representing the flow directions and the flow paths in the networking system needs to be prepared. Also, the runoff and the travel time in each cell should be calculated. 2.3.2 Flow direction In the GIS environment, the flow direction refers to the 3x3 grid cell rule using direction codes or the D8 grid method (Tarboton et al., 2009), as shown in Figure 5(a). To satisfy the water balance, a cell may receive incoming flows from 7 neighboring cells. Additionally, the water stored may only flow out to one downstream cell. In this study, recoding the flow directions was performed for simplification in building the computation program, as shown in Figure 5(b). Figure 5. Flow direction codes: (a) ArcGIS basic codes, (b) codes used in model construction. 2.3.3 Runoff Calculation The runoff yield in each grid cell was computed after the excess rainfall was calculated. The depths’ distribution within the study area was constructed using Inverse Distance Weight (IDW) in ArcGIS software. The data points used for applying the IDW method were obtained from the radar adjusted to the ground rain gage (ARR). Due to lack of equipment, it was the only ARR installed in each study area. It involved adjusting the total rainfall depth of a point radar by measuring the same coordinate as the position of the ARR to the total depth of rainfall. The depth was determined by the gage counted within the duration of an event to obtain a multiplication factor. Afterward, the multiplication factor was used to obtain other points of radar values within the study area. By this approach, it is assumed that the rainwater falls vertically from the sky to the earth. Therefore, the rainfall measured by the gage is recorded by the radar. Moreover, it is hypothesized that the storm occurs uniformly in time throughout the entire study area. The runoff in each cell is the volume of the excess rainfall in each time fraction, Δt, calculated as follows: 171 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 𝑄(𝑡)𝑖 = (𝑃𝑒(𝑡))𝑖 ∆𝑡 𝐴 (3) where Q(t) is the runoff in a cell at time step i (m3/s), Pe(t)I is the excess rainfall depth at time step i (m), Δt is the time fraction (s), and A the cell size (m2). Based on eq. (3), the runoff volume, V(t), in each cell can be written as follows: 𝑉(𝑡)𝑖 = (𝑃𝑒 (𝑡))𝑖 𝐴 (4) The excess rainfall, Pe(t) in eq. (3) was calculated using the NRCS-CN method as follows: (USDA, 2004a) 𝑃𝑒 = 𝑃 − 0.2𝑆 − 𝑆(𝑃−0.2𝑆) 𝑃+0.8𝑆 (5) where P is the rainfall depth (mm) and S is the maximum soil water retention parameter (mm) calculated by the equation below: 𝑆 = 25400 𝐶𝑁 − 254 (6) where CN is the curve number determined based on land characteristics and soil properties (the values were referred to USDA (2004b) and the antecedent rainfall condition (ARC) class (Table 1). The standard CN values were determined based on the ARC class II, i.e., CN(2). The CN values for ARC class I and III are calculated using eq. (7) and eq. (8) (Chow et al, 1988). 𝐶𝑁(1) = 4.2𝐶𝑁 (2) 10−0,058𝐶𝑁 (2) (7) 𝐶𝑁(3) = 23𝐶𝑁 (2) 10+0,13𝐶𝑁(2) (8) Table 1. Classification of antecedent runoff condition (McCuen, 1998; USDA, 2004b) ARC Class Soil condition 5-day antecedent rainfall depth (mm) Dry season Wet Season I Dry < 12.7 < 35.56 II Average 12.7 – 27.94 35.56 – 53.34 III Wet > 27.94 > 53.34 The Pe(t) and Qe(t) in eq. (3) in all cells for each incremental time step were computed using Matlab software, after converting the rainfall raster data proceeded in the GIS environment into ASCII data. 172 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 2.3.4 Travel Time Calculation Two kinds of flow were considered in calculating the flow travel time, including overland and channel. (1) Overland Flow Travel Time The overland travel time in any given cell j for each incremental time step i is estimated using Eq. (9). (𝑇𝑙)𝑗,𝑖 = 𝐿𝑗 (𝑣𝑙)𝑗,𝑖 (9) where Tl is the travel time (s), L is the flow length (m) (L = cell width in case the flow is in horizontal or vertical directions, or L = 20.5 of the cell width where the flow is in a diagonal direction). vl is the flow velocity (m/s), and the l subscript indicates the overland flow. The overland flow velocity in a cell can be estimated by combining the steady-state uniform flow expression with Manning's formula (Singh & Aravamuthan, 1996; Kang & Merwade, 2011). In the steady-state uniform flow, the unit width discharge in a cell can be written as follows: (Kang & Merwade, 2011) (𝑞𝑙)𝑗,𝑖 = (𝜑𝑙)𝑗,𝑖 . 𝐿𝑗 = 𝑦𝑗,𝑖 . (𝑣𝑙)𝑗,𝑖 (10) where q is the unit width discharge (m2/s), 𝜑 is the flux (m/s) as given by Eq. (11), and y is the flow depth (m). The flux is defined by the following equation: (𝜑𝑙 )𝑗,𝑖 = 𝑆𝑗,𝑖 𝐴𝑗 ×∆𝑡 (11) where S is the storage in the cell (m3), A is the cell’s area (m2), and Δt is the time step interval or rainfall time fraction (s). The overland flow velocity can be estimated using Manning’s formula (Chow et al., 1988) as given in eq. (12): 𝑣𝑙 = 𝑠 1 2⁄ × 𝑦 2 3⁄ 𝑛 (12) where s is the slope (m/m), and n is the Manning roughness coefficient. Rearranging Eq. (12) to obtain y and substituting it into eq. (10), the overland flow velocity can be written as follows: (𝑣𝑙)𝑗,𝑖 = 𝑠𝑗 0.3 × (𝜑𝑙)𝑗,𝑖 0.4 × 𝐿𝑗 0.4 𝑛𝑗 0.6 (13) For the sake of simplicity in programming, the velocity in eq. (13) is expressed differently as given by eq. (14). 173 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 (𝑣𝑙)𝑗,𝑖 = (𝐾𝑠𝑛 ) 0.6 × ((𝜑𝑙)𝑗,𝑖 × 𝐿𝑗 ) 0.4 (14) where Ksn is the Slope-Manning Coefficient written as: 𝐾𝑠𝑛 = 𝑠𝑗 0.5 𝑛𝑗 (15) Therefore, the overland flow travel time in any given cell is calculated using eq. (9) after substituting the velocity, eq. (14). (2) Channel Flow Travel Time As in overland flow travel time derivation, the travel time in any given cell for channel flow is estimated by the following equation: (𝑇𝑟 )𝑗,𝑖 = 𝐿𝑗 (𝑣𝑟)𝑗,𝑖 (16) where the subscript r refers to channel flow. The channel flow velocity in any given cell for each incremental time step can be estimated by combining the continuity equation with Manning's formula based on the assumption of a wide channel (Melesse, 2002; Kang & Merwade, 2011). The assumption is reasonable because the cell size is larger than the actual channel. The continuity equation for each incremental time step for channel flow is written as follows: 𝑆𝑗,𝑖 ∆𝑡 = 𝐴𝑗 × (𝑣𝑟 )𝑗,𝑖 = 𝑦𝑗,𝑖 × 𝐵𝑗 × (𝑣𝑟 )𝑗,𝑖 (17) where B is the width of the cell (m). The Manning's formula for channel flow is written as follows: (Melesse, 2002; Kang & Merwade, 2011) 𝑣𝑟 = 𝑠 1 2⁄ ×𝑅 2 3⁄ 𝑛 (18) where R is the hydraulic radius (m) as the ratio of the cross-sectional area to its wetted perimeter. By using a wide channel assumption, the hydraulic radius can be replaced by the flow depth (R = y). Similarly, the channel slope can be replaced by the slope as in overland flow (s = sj). Replacing these two flow parameters in eq. (18) and substituting it into eq. (17), then rearranging the substituted equation, the velocity expression for channel flow is obtained, as given by eq. (19). (𝑣𝑟 )𝑗,𝑖 = 𝑠𝑗 0.3 𝑛𝑗 0.6 × ( 𝑆𝑗,𝑖 𝐵𝑗∆𝑡 ) 0.4 (19) To organize the scripting of the computation program in a good pattern, eq. (19) can be rewritten in the same form as in the overland flow, eq. (14). 174 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 (𝑣𝑟 )𝑗,𝑖 = (𝐾𝑠𝑛 ) 0.6 × ((𝜑𝑟 )𝑗,𝑖 × 𝐵𝑗 × 𝐿𝑗) 0.4 (20) where Ksn is defined as the same as in eq. (14) and 𝜑𝑟 as the flux in wide channel given by (𝜑𝑟 )𝑗,𝑖 = 𝑆𝑗,𝑖 𝐴𝑗 × ∆𝑡 (21) Therefore, by substituting the velocity in eq. (20) into eq. (16) the channel flow travel time is obtained. 2.3.5 Routing Model Concept The concept of the routing is schematically depicted in Figure 6. Each cell, j, represents a single value of excess rainfall as a runoff yield in each incremental time step, i. The value is calculated based on the NRCS-CN method. The runoff from the most upstream cell is then routed to the outlet cell of the study catchment through a flow path constructed by flow direction connectivity using GIS software. During the routing processes, accumulations co-occur at cells that receive flows from upstream. This means that the outlet receives all these simultaneous flow accumulations at each time step. The two critical parameters in this model include the curve number (CN) and surface Manning roughness (n). The CN values are related to the amount of excess rainfall in each cell at every incremental time step. In this regard, the length of the time step, Δt, corresponds to the time fraction of the recorded rainfall. Furthermore, the n values are related to the speed of the flow. Since cell travel time is the critical issue, the consequence of selecting the grid cell size and the time step to the cell travel time is crucial in the proposed routing model development. The time step is used as the basis of a hydrograph. Additionally, the flow in a cell should only be the runoff and water from its neighboring cells stored within a duration that is equal to or greater than the given time step. This means that the travel time of the flow in the cell should at least be equal to the travel time in the neighboring cells. In the case of an unfavorable deviation, flows from upstream cells also contribute to the volumetric rate. This may cause inconsistency in determining the basis time as well as the calculated discharge. Therefore, this study mainly focused on constructing an algorithm to deal with this inconsistency. This was achieved using looping or iteration in the routing model. Importantly, the loop statement satisfied a cell travel time following the time step. The step block justifying looping for travel time at outlet cell is performed before executing the next time step (i+1). 175 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 DEM Rainfall Depth Slope (j) n (j) CN (j) Excess Rainfall (j,i) Area (j) Flow Length (j) Flow Direction (j) Inflow Volume (k,i) Hydrograph (joutlet,i) Soil Type LU/LC Velocity (j,i) Travel Time (j,i) Justifying Looping for Travel Time at outlet (joutlet,i) Routed Runoff (j,i) Routing (j,i) Runoff Volume (j,i) Outflow Volume (j,i) Storage (j,i) Variable ParameterProcess Calibrated Parameter ARC Figure 6. Routing model concept for hydrograph simulation To contextually apply the routing model concept, the algorithm was constructed and scripted using Matlab programming codes. The model results constructed without looping as the original model were also presented for comparison. 2.4 Model Performance Evaluation The proposed model algorithm was tested to simulate DRHs of three continuous storm events. One event was treated as a calibration while the other two were selected for validation. The CN and n values were treated as parameter calibrations. Furthermore, the calibration was performed manually based on the literature values of CN and n that might be appropriate for the study areas. The model performance was evaluated using statistical analysis criteria of the Nash- Sutcliffe efficiency (NSE), the relative error of peak flow (QpE), and the volume bias (VB) given by eqs. (22) to (24), respectively. 𝑁𝑆𝐸 = 1 − ∑ [(𝑄𝑠 )𝑖−(𝑄𝑜 )𝑖] 2𝑚 𝑖=1 ∑ [(𝑄𝑠 )𝑖−(�̅�𝑜)𝑖] 2𝑚 𝑖=1 (22) 𝑄𝑝 𝐸 = |𝑄𝑝𝑠 −𝑄𝑝𝑜 | 𝑄𝑝𝑜 𝑥100 (23) VB= 𝑉𝑠−𝑉𝑜 𝑉𝑜 𝑥100 (24) where Qs and Qo are the simulated and observed flows, respectively, (m 3/s), Qpo and Qps are the observed and simulated flows, respectively, (m3/s), and Vo and Vs are the observed and 176 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 simulated runoff volume total at the outlet, respectively, (m3). The approximate values for the acceptance of the model performance are greater than 0.75 for NSE and less than 25 % for both QpE and VB. 3. Results and Discussion 3.1 Using Model Algorithm The model program was constructed using Matlab software. The construction was based on the routing model concept (Figure 6), cell water balance (eq. 1 and eq. 2), and all the involved water flow components (eq. 9 – eq. 21). The algorithm for scripting the program is depicted in Figure 7. A new form of this model involves the application of loop or iteration statements, as enclosed by the red line in Figure 7. At the first step, indexed by i = 1, the flow parameters are calculated to justify the travel time in a cell. When the travel time in the outlet cell is less than the time step Δt, a looping process is started. The travel time is recalculated for each repetition. Furthermore, the looping process is automatically terminated when the accumulative travel time in the outlet cell is equal to or greater than the time step. The next step, i = 2, to the last are then consecutively executed. Each time the travel time in the outlet cell is less than the time step in each increment, the looping process is performed. By applying loop statements, the travel time is confined to the selected time step and acts as the basis for a simulated hydrograph. The next task involves justifying the discharge obtained at the given time step. Since the iterations in a loop result in more than one calculated outflow in the outlet cell, a single discharge value is counted for each incremental time step. The condition for outflow in eq. (1) as shown in line 2 of eq. (2) prevents direct summation of the iterative discharges in each looping. Although the numerator is Δt, the denominator depends on each iteration's travel time in a looping. Therefore, the final discharge of each looping is approximated by multiplying the average iterative discharges by time factor as follows: 𝑄𝑙 = ∑ 𝑄𝑚 𝑀 𝑚=1 𝑀 𝐹𝑇 (25) where Ql is the discharge in each time step (m 3/s), Qm is the iterative discharges in each looping (m3/s), m is the iteration index, M is the total iteration in each looping, and FT is the time factor calculated as the rainfall time fraction (seconds) divided by 60 seconds, or FT= Δt/60. 177 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 r,c (no. of rows and colums of grid area) times (# of time steps) Pe (r*c*times -- matrix of excess hytograph) A (r*c -- matrix of grid area) Sg (r*c -- matrix of slope grid) ∆t (Time step) L (r*c -- matrix of flow length) n (r*c -- matrix of Manning Roughness coef.) flodir (r*c --matrix of flow direction) (xo,yo) (outlet coordinate) FT (rainfall time fraction) Ksn = Sg0.5/n; no = n(xo,yo) Q1 = P1 /∆t (0.001*A) O1 = 0 S1 = Q1 ∆t φ1 = S1 /(A*∆t) vl1= Ksn0.6 (L*φ1)0.4 vr1= Ksn0.6 (B*L*φ1)0.4 T1 = traveltime(vl,vr,L,n) to = T1(xo,yo); Oo = O1(xo,yo) i=1 iter = 1 Plot(times,discharges) i= i+1 Qi = Pi /∆t (0.001*A) Oi = outflow(Ti-1,Si-1,∆t) Ii = inflow(flodir,Oi) Si = Qi*∆t + Si-1 + Ii*∆t - Oi*∆t φi = Si /(A*∆t) vli= Ksn0.6 (L*φi)0.4 vri= Ksn0.6 (B*L*φi)0.4 Ti = Travel time(vl,vr,L,n) Oo = Oi(xo,yo) to = Ti(xo,yo) iter = loop i≤ times? NO START YES END to<∆t? loop = iter Sloop = Si tloop = to tcum = to Ocum = Oo tcum<∆t? loop = loop+1 Oloop = outflow(Tloop,Sloop,∆t) Iloop = inflow(flodir,Oloop) Sloop = Sloop + Iloop *∆t - Oloop∆t φloop = Sloop/(A*∆t) vlloop = Ksn0.6 (L*φloop)0.4 vrloop = Ksn0.6 (B*L*φloop)0.4 Tloop = traveltime(vlloop,vrloop,L,n) Ocum = Ocum+Oloop(xo,yo) tloop = Tloop(xo,yo) tcum = tcum+tloop Discharge(i)=Oo*FT /iter NO iter = loop Oi = Oloop Si = Sloop Ti = Tloop Oo = Ocum NO YES YES Figure 7. Algorithm of the proposed model. A block bounded by the red line is the looping mechanism. 178 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 3.2 Model Testing; Calibration and Validation Three isolated storm events in each study area were selected to test the performance of the proposed routing model. The summary of measured rainfall and flow data for the events is shown in Table 2. Table 2. Summary of the rainfall and flow data in the study areas Study Catchment Storm Event Duration (min.) Time step (Δt) (min.) Total rainfall depth gauge (mm) Tot. 5 days antecedent rainfall depth (mm) Peak flow (m3/s) Time to peak (min.) Code Date Start time Sumur Mbandung B1 1-Mar-17 13:45 135 15 42.8 53.6 4.09 45 B2 3-Feb-17 2:10 150 15 39.8 134.2 4.63 45 B3 15-Jan-17 16:05 45 15 22.8 61.8 2.97 45 Sempor S1 1-Mar-17 16:50 165 15 115.2 188.8 4.17 60 S2 15-Mar-17 15:25 45 15 24.4 78.8 0.76 30 S3 9-Nov-16 14:15 300 15 147.2 7.8 1.20 150 The model was first tested for each study catchment using the B1 event for Sumur Mbandung and S1 for Sempor. For this reason, the CN and n values have been selected as in columns Test Value of Table 3 and Table 4. The values are adjusted for LULC and HSG in Sumur Mbandung and Sempor catchments, respectively. The tables also include the calibrated CN and n values. Simulated hydrographs of the initial test are shown in Figure 8, while their statistical evaluations are presented in Table 5. Table 3. Curve Number (CN) and Manning coefficient (n) values for the study area in Sumur Mbandung Catchment (initial test and calibration) LULC HSG Curve Number, CN Manning Coefficient, n Range CN(2)* Test value Calibrated value Range** Test value Calibrated value CN(2) CN(1) CN(3) Scarcely forest C 79 - 86 86 79 61 90 0.035 - 0.160 0.035 0.195 Seasonal agriculture C 77 - 84 84 83 67 92 0.030 - 0.500 0.030 0.205 Residential lots C 82 82 83 67 92 0.011 - 0.035 0.110 0.095 Grass C 80 - 89 89 80 63 90 0.030 - 0.500 0.030 0.175 Non-irrigated paddy field C 90 - 95 95 91 81 96 0.030 - 0.500 0.030 0.095 Open area C 90 - 95 95 95 89 98 0.025 - 0.033 0.025 0.075 Water - 100 100 98 98 98 0.025 - 0.033 0.025 0.040 *USDA (2014b); ** Kang & Merwade (2011) 179 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Table 4. Curve Number (CN) and Manning coefficient (n) values for the study area in Sempor Catchment (initial test and calibration) LULC HSG Curve Number, CN Manning Coefficient, n Range CN(2)* Test value Calibrated value Range** Test value Calibrated value CN(2) CN(1) CN(3) Scarcely forest B 55 - 73 73 50 30 70 0.035 - 0.160 0.035 0.155 Densely forest C 55 - 73 73 50 30 70 0.035 - 0.160 0.035 0.155 Seasonal agriculture C 69 - 80 80 70 49 84 0.030 - 0.500 0.03 0.175 Zalacca palm plantation C 69 - 80 80 70 49 84 0.030 - 0.500 0.03 0.175 Residential lots C 74 74 74 54 87 0.011 - 0.035 0.011 0.075 Irrigated paddy field C 90 - 95 95 47 27 67 0.030 - 0.500 0.030 0.055 Water - 100 100 98 98 98 0.025 – 0.033 0.025 0.030 *USDA (2014b); ** Kang & Merwade (2011) Figure 8. Initial test of model results for event B1 and S1 Table 5. Results of initial test models for B1 and S1 events Study Area Event Vo Vs Qpo Qps NSE QpE VB Tpo Tps (m3) (m3) (m3/s) (m3/s) % % (min.) (min.) Sumur Mbandung B1 15,729 9,105 4.04 2.23 0.30 44.89 -42.11 45 45 Sempor S1 22,335 20,543 4.166 3.71 0.95 10.96 -8.02 60 45 Vo = Observed runoff volume, Vs = Simulated runoff volume, Qpo = Observed peak flow, Qps= Simulated peak flow, NSE = Nash-Sutcliffe coefficient, QpE = Relative error of peak flow, VB = Volume bias, Tpo = Observed time to peak and Tps =Simulated time to peak. As shown in Figure 8, initially, the simulated DRHs of the two events have similar trends. The statistical evaluations of the model performances were critical in this regard. In the Sumur Mbandung catchment, the NSE is less than 0.75, while the total runoff volume that reaches the outlet is underestimated by 42%. Similarly, the peak flow is underestimated by Sumur Mbandung Sempor 180 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 45%. In the Sempor catchment, the NSE is greater than 0.75, while the total runoff volume is underestimated by 8%. Additionally, the peak flow is underestimated by 11%. Based on the initial test results, the proposed algorithm performs well in simulating DRHs. The calibration on the two events is then performed manually with CN and n as the parameters. The study area at Sumur Mbandung catchment has a latosol classified into the hydrological soil group C. Furthermore, the regosol in Sempor catchment is categorized in hydrological soil group B. Both events have a total of 5 days of antecedent rainfall depths greater than 53.34 mm. Therefore, they are classified as ARC III, and the CN(3) values of all land characteristics in the study areas were used to calculate the runoff yield. The values of CN and n for any combination of LULC and soil types after calibration in Sumur Mbandung and Sempor are presented in Table 3 and Table 4, respectively. The calibrated and the validated model DRHs for all events in both study areas are shown in Figure 9. The statistical evaluations of the model performances are presented in Table 6. Calibrated model DRHs (B1 and S1 events) and validated models (B2, B3, S2, and S3 events) Sumur Mbandung Sempor Figure 9. 181 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Table 6. Statistical evaluation results after calibration and validation for all selected events in the study areas Study Area Event Vo Vs Qpo Qps NSE QpE VB Tpo Tps (m3) (m3) (m3/s) (m3/s) % % (min.) (min.) Sumur Mbandung B1 15,729 16,348 4.04 4.17 0.97 3.06 3.94 45 45 B2 20,812 18,384 4.63 4.40 0.89 4.91 -11.67 45 45 B3 5,490 4,912 3.14 2.97 0.96 5.62 -10.53 45 45 Sempor S1 22,335 22,075 4.16 4.20 0.98 0.86 -1.16 60 45 S2 935 1,103 0.68 0.74 0.99 4.59 -17.95 30 30 S3 11,435 11,394 1.26 1.20 0.88 4.61 -0.36 150 150 Vo = Observed runoff volume, Vs = Simulated runoff volume, Qpo = Observed peak flow, Qps= Simulated peak flow, NSE = Nash-Sutcliffe coefficient, QpE = Relative error of peak flow, VB = Volume bias, Tpo = Observed time to peak and Tps =Simulated time to peak. After calibrations using CN(3) and n values, the shape of the simulated DRHs for both B1 and S1 are closer to the observed ones than initial test results. The model performance measures give excellent results. Precisely, the NSE is 0.97 for the B1 event and 0.98 for S1. The QpE is overestimated by 3% for B1 and 1% for S1. The VB is underestimated by 4% for B1 and is overestimated by 1% for S1. Although the time to peak in the simulated DRH for S1 is one step before its observed value, the routing model might produce the best fit of calibrated hydrographs. Model validations for the other four events, two in each study area, also produce comparable results. For instance, B2 and B3 in the Sumur Mbandung study area, both events have a total 5-day rainfall depth of greater than 53.34 mm. Therefore, CN(3) values were used for running the model, keeping the calibrated n values as they are. The validated model results for both events had good NSE, specifically, 0.89 for B2 and 0.96 for B3. Other model performance measures are also within the range of acceptance. In the Sempor catchment, validation for S2 was also performed using CN(3). The validated model result shows the best fit DRH compared to the observed one with the NSE is 0.99, the QpE is overestimated by 5%, and the VB is underestimated by 18%. The time to peak is also the same as the time to peak of the observed DRH, which is 30 minutes. S3 lasts in the longest duration of 300 minutes with the highest total rainfall depth gauge (147.2 mm) of all selected events in the Sempor area. However, it produced lesser runoff volume at the outlet compared with S1(Vo in Table 6). This is attributed to the hydrological condition before the event. As shown in Table (2), the total 5-day rainfall depth 182 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 before the S3 event is 7.8 mm. At low antecedent rainfall conditions, vegetation covers potentially retain a large amount of rainwater for the next event, especially in the early stage. The soil is in quite dry conditions and therefore, infiltration reduces runoff production. In applying the model for validation for S3, the CN(1) values were used as ARC laid on class I. The validated result shows that the simulated DRH for S3 is relatively similar to the observed one and its model performance measures are within the range of acceptance (Table 6). 3.3 Comparison with Original Model The original model, which does not involve looping and iteration mechanisms in the routing algorithm, was also tested using the calibrated CN and n values for all the selected events in the study areas. The results of the simulated DRHs of the model are shown in Figure 10. Original model results for the selected storm events in the study areas Figure 10 shows that the first model results of all simulated DRHs end up with much longer estimated runoff time than the observed ones. Also, all the peaks are consistently Figure 10. 183 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 underestimated. The results are principally consistent with Kang & Merwade (2011), which established that when the cell size is small and the time step, Δt, is long enough to produce travel time shorter than the time step, the model hydrograph always has a long basis time with low peak. To overcome this issue, the critical cell travel time (CCT) was used to select a time step for the corresponding grid cell size to ensure the flow travel time is greater than the selected step. Additionally, the proposed routing model with looping adjusts the cell travel time to the time step. 4. Conclusion This paper discussed hydrograph simulations constructed from a runoff routing model based on the cell water balance. The study used travel time as a conditional preference. The original model's simulation results showed consistently long runoff times and low peaks compared to the observed hydrographs. To overcome this problem, a new algorithm was developed using a looping mechanism. The routine guarantees that cell travel time is always greater than the selected time step. The simulation results of the new algorithm showed resemblance in shapes to observed hydrographs. The model performed well with excellent measures. Moreover, all simulated direct runoff hydrographs had a Nash-Sutcliffe coefficient of greater than 0.75, the volume bias, and peak discharge error of less than 25%. Although the results of applying the proposed model algorithm are promising, this paper only presents an initial effort. Therefore, future studies should focus on promoting the model validity. These can be achieved by testing the model for other catchments with different hydrological conditions. Also, the model can be tested in terms of grid cell size and time step variations. Conflict of Interest The authors declare that there is no conflict of interest with any financial, personal, or other relationships with other people or organizations related to the material discussed in the article. Acknowledgments This research was financially supported by the Directorate General of Higher Education, the Ministry of Education and Cultural Republic of Indonesia. The authors express gratitude to the Faculty of Geography, Gadjah Mada University, for facilitating the research. The authors also thank Mrs. Nurhaida, a postgraduate student at the Faculty of 184 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Basic Sciences, Gadjah Mada University, for valuable assistance in constructing the computation program. References Ahmad, I., Verma, V., & Verma, M. K. (2015). Application of the Curve Number Method for Estimation of Runoff Potential in GIS Environment. 2015 2nd International Conference on Geological and Civil Engineering, IACSIT Press, Singapore, 80, 16–20. https://doi.org/10.7763/IPCBEE. Al-Juaidi, A. E. M. (2018). A simplified GIS-based SCS-CN method for the assessment of land-use change on runoff. Arabian Journal of Geosciences, 11(269). https://doi.org/10.1007/s12517-018-3621-4 Asfaw, A., Maher, K., & Shucksmith, J. D. (2018). Modelling of metaldehyde concentrations in surface waters : A travel time based approach. Journal of Hydrology, 562, 397–410. https://doi.org/10.1016/j.jhydrol.2018.04.074. Bansode, A. ., & Patil, K. A. (2014). Estimation of runoff by using SCS curve number method and arc GIS. International Journal of Scientific & Engineering Research, 5(7), 1283–1287. Beven, K. (2012). Rainfall-Runoff Modelling: The Primer (2nd ed.). Chichester, UK: John Wiley & Sons, Ltd. Chow, V. T., Maidment, D. R., & Mays, L. R. (1988). Applied Hydrology. Singapore: McGRaw-Hill Book Co. Du, J., Xie, H., Hu, Y., Xu, Y., & Xu, C. Y. (2009). Development and testing of a new storm runoff routing approach based on time variant spatially distributed travel time method. Journal of Hydrology, 369(1–2), 44–54. https://doi.org/10.1016/j.jhydrol.2009.02.033. Gonzalez, A., Temimi, M., & Khanbilvardi, R. (2015). Adjustment to the curve number (NRCS-CN) to account for the vegetation effect on hydrological processes. Hydrological Sciences Journal, 60(4), 591–605. https://doi.org/10.1080/02626667.2014.898119. Kang, K., & Merwade, V. (2011). Development and application of a storage-release based distributed hydrologic model using GIS. Journal of Hydrology, 403(1–2), 1–13. https://doi.org/10.1016/j.jhydrol.2011.03.048 Maina, C. W., & Raude, J. M. (2016). Assessing Land Suitability for Rainwater Harvesting Using Geospatial Techniques : A Case Study of Njoro Catchment, Kenya. Applied and Environmental Soil Science, 2016. https://doi.org/dx.doi.org/10.1155/2016/4676435. McCuen, R. H. (1998). Hydrologic Analysis and Design (2nd ed.). New Jersey: Prentice- Hall. Melesse, A. M. (2002). Spatially distributed storm runoff modeling using remote sensing and geographic information system. The University of Florida. Melesse, A. M., & Graham, W. D. (2004). Storm Runoff Prediction Based on a Spatially Distributed Travel Time Method Utilizing Remote Sensing and GIS. Journal of the American Water Resources Association, 40(4), 863–879. Ponce, V. M., & Hawkins, R. H. (1996). Runoff Curve Number: Has It Reached Maturity? 185 Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185 Journal of Hydrologic Engineering, 1(1), 11–19. Rawat, K. S., & Singh, S. K. (2017). Estimation of Surface Runoff from Semi-arid Ungauged Agricultural Watershed Using SCS-CN Method and Earth Observation Data Sets. Water Conserv Sci Eng, 1, 233–247. https://doi.org/10.1007/s41101-017-0016-4. Rohman, A., Comber, A., & Mitchell, G. (2019). Evaluation of Natural Flood Management using Curve Number in the Ciliwung Basin, West Java. In P. Kyriakidis, D. Hadjimitsis, D. Skarlatos, & A. Mansourian (Eds.), Geospatial Technologies for Local and Regional Development : Proceedings of the 22nd AGILE Conference on Geographic Information Science (Vol. 2018, pp. 2–5). Cham, Switzerland: Springer. Satheeshkumar, S. Venkateswaran, S., & Kannan, R. (2017). Rainfall-runoff estimation using SCS–CN and GIS approach in the Pappiredipatti watershed of the Vaniyar sub-basin, South India. Model. Earth Syst. Environ., 3(24). https://doi.org/doi.org/10.1007/s40808- 017-0301-4. Singh, V. P., & Aravamuthan, V. (1996). Errors of kinematic-wave and diffusion-wave approximations for steady-state overland flows. Catena, 27, 209–227. Singh, V. P., & Woolhiser, D. A. (2002). Mathematical Modeling of Watershed Hydrology. Journal of Hydrologic Engineering, 7, 270–292. Soulis, X., & Valiantzas, J. D. (2012). SCS-CN parameter determination using rainfall-runoff data in heterogeneous watersheds – the two-CN system approach. Hydrology and Earth System Sciences, 16, 1001–1015. https://doi.org/10.5194/hess-16-1001-2012. Tarboton, D. G., Schreuders, K. A. T., Watson, D. W., & Baker, M. E. (2009). Generalized terrain-based flow analysis of digital elevation models. 18th World IMACS/MODSIM Congress. Retrieved from http://mssanz.org.au/modsim09. USDA. (2004a). Estimation of Direct Runoff from Storm Rainfall. In the National Engineering Handbook, Part 630: Hydrology. USDA. (2004b). Hydrologic Soil-Cover Complexes. In the National Engineering Handbook, Part 630: Hydrology. Vojtek, M., & Vojteková, J. (2016). GIS-Based Approach To Estimate Surface Runoff In Small Catchments : A Case Study. Quaestiones Geographicae, 35(3), 97–116. https://doi.org/doi:10.1515/quageo–2016–0030. Zhao, L., & Wu, F. (2015). Simulation of Runoff Hydrograph on Soil Surfaces with Different Microtopography Using a Travel Time Method at the Plot Scale. PLoS ONE, 10(6), 1– 14. https://doi.org/10.1371/journal.pone.0130794. 1. Introduction 2. Methods 2.1 Area of Study 2.2 Data Resources 2.3 Model Construction 2.3.1 Water Balance Concept 2.3.2 Flow direction 2.3.3 Runoff Calculation 2.3.4 Travel Time Calculation 2.3.5 Routing Model Concept 2.4 Model Performance Evaluation 3. Results and Discussion 3.1 Using Model Algorithm 3.2 Model Testing; Calibration and Validation 3.3 Comparison with Original Model 4. Conclusion Acknowledgments References