A New Algorithm For The Grid Cell-Based Runoff Routing Model Based on Travel Time Concept

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 NashSutcliffe coefficient was greater than 0.75 for all DRH during validation, and the volume bias and peak discharge error were less than 25%.

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 GISbased SCS-CN are essential in estimating runoff within catchments of similar geohydrological 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)  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.
Baina Afkril et al / GEOSI Vol 5 No 2 (2020)   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 timeinvariant. 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 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.

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 km 2 . 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 km 2 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.

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 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.  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.
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.

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) 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.
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.

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).

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: where Q(t) is the runoff in a cell at time step i (m 3 /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 (m 2 ). Based on eq. (3), the runoff volume, V(t), in each cell can be written as follows: The excess rainfall, Pe(t) in eq. (3) was calculated using the NRCS-CN method as follows: (USDA, 2004a) where P is the rainfall depth (mm) and S is the maximum soil water retention parameter (mm) calculated by the equation below: 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).

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).
where Tl is the travel time ( where q is the unit width discharge (m 2 /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: where S is the storage in the cell (m 3 ), A is the cell's area (m 2 ), 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): 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: For the sake of simplicity in programming, the velocity in eq. (13) is expressed differently as given by eq. (14).
where Ksn is the Slope-Manning Coefficient written as: 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: where the subscript r refers to channel flow.
where B is the width of the cell (m).
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 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).
where Ksn is defined as the same as in eq. (14) and as the flux in wide channel given by Therefore, by substituting the velocity in eq. (20) into eq. (16) the channel flow travel time is obtained.

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. 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). 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.

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.

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. 9eq. 21). The algorithm for scripting the program is depicted in Figure 7. (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: where Ql is the discharge in each time step (m 3 /s), Qm is the iterative discharges in each looping (m 3 /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. Baina Afkril et al / GEOSI Vol 5 No 2 (2020) 160-185

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. 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.   Figure 8. Initial test of model results for event B1 and S1 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, 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. 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 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).

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 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.

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.