An innovative approach of GSSHA model in flood analysis of large watersheds based on accuracy of DEM, size of grids, and stream density

Distributed modeling approach may have much better performance and accuracy compared with lumped-parameter hydrologic models. The main goals of this research are: investigating the possibility of combining distributed hydrological models with an one-dimensional hydraulic model and simulating waterways in large watersheds with limited hydrological and hydraulic data. Then performing sensitivity analysis on different parameters in order to identify the parameters containing the major influences on results. In the current research, an innovative approach in Gridded Surface Subsurface Hydrologic Analysis (GSSHA) model, the cross-sections of all 414 waterways in the 3450 km2 Karvandar watershed, used for flow routing calculations, are uniquely extracted. Then, the effect of three essential factors are evaluated. These factors are accuracy of the digital topographic model, cell size of grid network, and density of streams, on the results of GSSHA model simulations. This watershed is located in southeastern Iran, has a dry climate with limited available hydrological data. Results showed that peak discharges obtained from the GSSHA model, developed based on a DEM with a spatial resolution of 12.5 m, are slightly (< 4%) lower than the corresponding values in the GSSHA model with a 30 m DEM resolution. This fact confirms that the use of the topographic model with a lower spatial resolution has no substantial effects on the accuracy of simulation. Also, the peak discharges increased significantly (44% to 57%) by increasing the density of waterways in the GSSHA model. Furthermore, results showed that peak discharge obtained from three models with grid cell sizes of 100, 150, and 200 m (base model), are close together. Comparing with two models of coarser grids (250 and 300 m), significant differences observed, which indicated that the grids larger than 200 m could induce substantial errors in results.


Introduction
For many hydrologic studies, distributed modeling approaches may have many significant advantages in rainfall-runoff models compared with traditional lumpedparameter hydrologic models. Gridded Surface Subsurface Hydrologic Analysis (GSSHA) is a physically-based hydrologic model. This model simulates water surface flows in watersheds with miscellaneous runoff production processes, including runoff due to excess infiltration and saturated aquifers, as well as direct interaction between waterways and saturated groundwater; using mass-conserving equations. The hydrologic components are closely related to ensure overall mass balance (Downer et al. 2002a). GSSHA uses similar two-step explicit finite volume schemes to route the flow for one-dimensional (1-D) channels and two-dimensional (2-D) overland flow. Compared with the kinematic wave 33 Page 2 of 19 approach, this diffusive wave approach allows GSSHA to route runoff through pits or depressions and areas of adverse slope. Despite the widespread use of the Muskingum method (Farzin et al., 2018, Node Farahani et al., 2018, Farahani et al., 2019, the diffusive wave approach has high accuracy. The Manning formula relates flow depth to discharge and is employed in many studies (Farzin et al., 2021). Features of the GSSHA model include 1-D infiltration, 1-D streamflow, 2-D overland flow, and 2-D groundwater flow.
The GSSHA model has been successfully applied to diverse problems and can provide data necessary for the potential effects of projects, land-use change, best management practices (BMPs), climate change, and related issues (Downer 2008). Downer and Ogden (2004) calibrated and verified the GSSHA model versus outlet discharges. The model accurately reproduces soil moisture values during the growing season. With additional enhancements, soil moisture during the non-growing season is obtained within a root-mean-square error of 0.1. Downer and Ogden (2004) tested the GSSHA model on a watershed with different runoff mechanisms. Results showed that the model is helpful to investigate the stream flow-producing processes in the watershed scale. Jaiswal et al. (2020) compared two conceptual models including TANK and Australian water balance model (AWBM) and a physically distributed but lumped on HRUs scale SWAT model for Tandula basin of Chhattisgarh (India). They deduced that distributed methods could provide convenient results. Kalin and Hantush (2006) analyzed two distributed watershed-scale sediment models (Kinematic Runoff and Erosion (KINEROS-2) and GSSHA) by application on an experimental watershed. Even though GSSHA's flow component slightly performed better than KINEROS-2, the latter outperformed GSSHA in simulations for sediment transport. Sharif et al. (2010) presented a hydrologic analysis of a flood event in a small urbanized watershed using the GSSHA model in Austin, Texas. Observed discharge was compared to hydrographs obtained from the accuracy of model decreases by increasing the grid sizes. Surosoa et al. (2013) modeled flood hazard with GSSHA in the South Sumatra Province of Indonesia. Based on the risk map, several flood control measures are identified: canalization, polder, retention pond, and infiltration measurement for downstream areas; detention basin and dam construction for midland areas; reforestation for upstream areas. Zhang et al. (2014) detected the impacts of optimization on the Yong-Ding watershed in China by comparing the land use pattern characteristics using the GSSHA model. The results of the hydrological evaluation revealed that instead of the land-use location control, the land-use composition and configuration control might be a more powerful method for minimizing the negative hydrological impact of urbanization. Yang et al. (2016) analyzed flash flooding in small urban watersheds, with a major focus on rainfall variability, antecedent soil moisture, and urban stormwater management infrastructure over Harry's Brook watershed in the USA, during 2005USA, during -2006. Hydrologic modeling simulation with GSSHA indicated that the sensitivity of hydrologic response to spatial rainfall variability decreases by increasing storm intensity. Temporal variability of rainfall is relatively more important than spatial variability, especially for extreme storm events. Tauhidur and Adel (2016) determined the physical, social, and overall composite flash flood vulnerability indices (CFVI) for 153 neighborhoods of Riyadh, Saudi Arabia. The CFVI map showed that the central and southern half of the city are highly vulnerable; northern and northeastern districts are moderate to highly vulnerable, and the mountainous western neighborhoods are the least vulnerable to flash flooding. Furl et al. (2018a) assessed the use of satellite-based precipitation products versus a Next-Generation Radar (NEXRAD) Stage IV product during big storms across three spatial scales ranging from ~ 200 to 10,000 km 2 using a GSSHA model. As a final result, the products significantly underestimated strong storm events in all spatial scales. Furl et al. (2018b) suggested that storm location and its motion had an important role in the peak flows. When the intense headwater cells were started over the middle portion of the basin, peak flows were greatly attenuated, while the runoff volume remained constant. Pradhan and Loney (2018) analyzed the unit hydrograph peaking factor (UHPF) using the GSSHA model. They concluded that design events with return periods greater than 5-year are needed for the UHPF to fall within the guidance range and that UHPF becomes less sensitive to rainfall intensity with increasing accumulation time. Lugina et al. (2018) analyzed the effect of moving rainstorms in increasing peak discharge in the 2013 flood event in Ciliwung Rivers Basin (CRB) based on 13 scenarios of synthetic rain, including three scenarios for nonmoving rainfall and ten scenarios of the moving rainfall to-downstream or to-upstream direction, which all based on the equivalent volume. The GSSHA hydrological simulation was conducted, and the results showed that moving rainfall has important effects in increasing the peak discharge and decreasing time to peak discharge in CRB. Hwan et al. (2020) evaluated the applicability of GSSHA for long-term discharge and water quality simulation for the ungauged Peace Dam Watershed in South Korea. Results showed good capabilities to be used as a watershed model even with its overestimation of peak discharges for small storms. Li et al. (2020) tested a data-driven rainfall-runoff model using a long-short-term-memory (LSTM) network for a watershed in Houston, Texas, for severe flood events. The model results were also compared with the output of a process-driven GSSHA model. Results showed that the data-driven model was more efficient in terms of prediction and calibration.  and  utilized a new method to set the flood mitigation strategies in a coastal wastewater treatment plant in New York City, USA. The strategies are selected and the reliability of the Hunts Point plant is estimated with (and without) using best management practices (BMPs). The results present 46% increase on reliability by using BMPs. Anaraki et al. (2021) used data-driven methods for uncertainty analysis of flood frequency under climate change conditions. In flood frequency, analyzing the hydrological modeling is one of the most important sources of uncertainty. Zhou et al. (2021), using a flood frequency analysis framework with the GSSHA model, examined the role of rainfall spatial and temporal variability in flood frequency in the highly urbanized Dead Run watershed outside of Baltimore, USA. The impact of impervious areas on flood response decreased by increasing rainfall return period. For extreme storms, maximum discharge was closely linked to the spatial structure of rainfall. de Arruda Gomes et al. (2021) analyzed the Capibaribe River basin's hydrological response to extreme rainfalls and its impact on the city of Recife in Brazil using three lumped and distributed hydrologic and hydraulic models. A reservoir was also considered to investigate its effect on watershed hydrology. The study showed that the reservoir reduced the flood severity by 70.3% for the 2011 event, but this event still caused a flood covering an area of 6.01 km 2 . In order to investigate the effects of urbanization on increasing the flood risk in urban areas, Feng et al. (2021) applied six simulated land-use scenarios into coupled distributed hydrologic and hydraulic models. Results indicated that urbanization increases surface runoff and river discharge and decreases time to peak discharge. Karamouz and Mahani (2021) used Gridded surface subsurface hydrologic analysis (GSSHA) to model flood delineation in a coastal city of the USA. They showed that the development of an uncertainty-based flood inundation modeling could significantly improve the ability for better flood preparedness planning in coastal cities.
A review of the previous researches shows that they include evaluating and analyzing the sensitivity of GSSHA results in distributed hydrologic simulations of flood events with variable spatial and temporal parameters of rain, for different grid sizes and use of satellite rain models data. Still, less attention has been paid to one-dimensional hydraulic routing of flow by diffusive wave method in the GSSHA model and its vital role in the results of distributed hydrological simulations. Also, there was no idea concerning modeling the cross-sections of waterways, especially in large watersheds, which are necessary for hydraulic routing in GSSHA. In the current study, using an innovative approach, the unique cross-sections of all waterways in a large watershed of 3450 km 2 , including 414 waterways with a drainage area more than 5 km 2 , were extracted from the digital elevation model (DEM). This approach has been used in different sensitivity analyses in the models with the variable density of streams.
So far, the integrated hydrological and distributed hydraulic conceptual model has not been used in large watersheds due to the several of required parameters and the unavailability of most of these parameters. In this study using innovative approaches and in terms of uncertainty, and favorable results have been obtained. The most important innovation used in this research, as emphasized in the article, is the extraction of unique cross-sections for 414 different waterways with a drainage surface of more than 5 km 2 , which has not been done in other researches. As explained in the article, the comparison of the influence of modeling of waterways in one-dimensional hydraulic routing calculations using the diffusive wave method in three different cases (1) modeling the main river and its tributaries, (2) all Waterways with a drainage area of more than 5 km 2 , and (3) without modeling the waterways and performing hydraulic routing. The results have a very wide range of changes and indicate a very high sensitivity of the model to how performing calculations of hydraulic routing.

Infiltration model
The Green-Ampt (GA) equation has been used as an infiltration model in the current study. Four soil hydraulic parameters are needed for modeling the infiltration process using the GA equation: a. Soil saturated hydraulic conductivity, K s (cm h −1 ). b. Soil capillary suction head parameter, S f (cm). c. Effective porosity, θ e . d. Initial soil moisture content, θ.
The GA equation is as follows: in which f t , infiltration loss at time t, K, φ, θ i , (φ − θ i ), S f , F t are saturated hydraulic conductivity, saturated water content, initial water content, volumetric moisture deficit, wetting front suction, and cumulative loss at time t, respectively. This model is based on the following assumptions: • Homogeneous soil with a uniform θ i . • Pressure head at the wetting front, H f , is constant.
(1) The Green-Ampt with soil moisture Redistribution (GAR) method expands the capability of the GA method by redistributing soil moisture during non (or low intensity rainfall) periods. This allows infiltration capacity to recover for the next burst of storm intensity, and makes the GAR method suitable for simulating multiple rainfall events in series.
Considerable previous studies have been conducted to relate soil infiltration parameter values to soil textural classification. Hydraulic conductivities for all GA-based approaches are half of the saturated values (Rawls et al. 1982). These values have great uncertainty and low correlation with textural classification and soil texture. However, these values are useful and supply an initial estimate of infiltration parameters.
The standard approach in creating GSSHA model is the use of digital soil textural classification data to develop an index map of soil types. Soil textural maps may be combined with land use or vegetation maps. Table 1 is used to assign initial parameters to the soil types in the index map of current study.

1-D Routing scheme in GSSHA
The 1-D channel routing scheme is shown in Fig. 1. Inter-cell flows, Q i − 1/2 and Q i + 1/2 (m 3 /s) in longitudinal, x, direction are computed from depths, d, at the nth time level using the Manning equation for the head discharge relationship (Downer et al. 2002b): where n is roughness coefficient, A is channel cross sectional area the area (m 2 ), R is the hydraulic radius (m), and S f is the friction slope (m), calculated in x direction as (Downer et al. 2002b): where d is water depth (m), S ox is the land surface slope in the x direction. For negative slopes the energy head loss is used to calculate the discharge as (Downer et al., 2002b): The flow direction is determined around each node and the locally upstream cell properties are used to compute the flow discharge. This simple local determination of the upstream cells prevents flood routing in channels with adverse slopes. This method also allows better simulations of backwater effects.
Inter-node fluxes are used to calculate the flow volume, V, in each node as (Downer et al., 2002b): where q lat (m 2 /s) is the lateral inflow from the overland flow cells adjacent to the node, and q recharge (m 2 /s) is the exchange between the groundwater and channel. These new volumes are used to compute nodal values of A, d, and wetted perimeter at the n + 1 time level. Calculations proceed from the upstream boundary to the downstream boundary.

2-D routing scheme in GSSHA
Overland flow in 2-D GSSHA model is the same methods described for 1-D channel routing, except that the where p ij and q ij are discharges in x and y direction. Depths in each cell are calculated for the (n + 1)th time step based on the flows for each cell (Julien and Saghafian 1991): By this original formulation, two additional methods of solving the equations have been proposed, an Alternating Direction Explicit scheme (ADE) and an ADE scheme with an extra Predictor-Corrector step (ADEPC) (Downer et al., 2000, b).
In ADE scheme, inter-cell flows are first calculated for x-direction according to Eq. (6). Depths in each row are updated based on the flows for x-direction (Downer et al., 2002b): Inter-cell flows in the y-direction are computed using the updated depths (Downer et al., 2002a, b):  Depths in each column are updated based on the flows in the y-direction (Downer et al., 2002b): In the ADEPC method, additional steps are added to improve the accuracy and stability. As before, during each sweep, either by rows or by columns, the head is estimated based on the calculated flows [Eqs. (9) through (11)]. Next, using the updated depths, the flow is calculated for the n + 1 time step (Downer et al., 2002b): The discharge for time step n + 1/2 is the average discharges obtained from n and n + 1 time steps, defined by Eq. (13) (Downer et al., 2002b): The discharges are then used to update the original depths [in Eqs. (9) and (11)].
The flowchart for the current study is presented in Fig. 2.

Innovative approach of hydraulic routing in GSSHA
In this study, as an innovative approach at the level of large GSSHA distribution model, the cross-sections of all 414 waterways of the Karavandar River hydrographic network are uniquely extracted from DEM and used in three states of modeling the flood routing calculations by Diffusive Wave method (Fig. 3): 1. Modeling main Karvandar River and major tributaries only. 2. Modeling the streams with drainage areas greater than 5 km 2 . 3. Modeling non-stream at the entire watershed area.
Some of the selected sample cross-sections used in 1-D hydraulic routing by diffusive Wave method, moving from upstream to downstream, are shown in Fig. 4.

Study watershed and available data
The study watershed, Karvandar River watershed at Daman hydrometric station, southeast of Iran, has 3450 km 2 drainage area. The average slope of watershed is 16.6%, and its average elevation is 1300 m above mean sea level.
Daman hydrometric station has 45 years of daily and peak discharges but no flood hydrographs. There are five  The location of rain gauge stations and Thiessen polygons on the physiographic model of the Karvandar River watershed is illustrated in Fig. 9.

Results and discussion
Daman hydrometric station situated at the watershed outlet records the discharges corresponding to the rainfall and storm events of the watershed. Five storm events and related floods are analyzed to calibrate the GSSHA Distributed Hydrological Model parameters and validate the model. The solution grid is performed using the inverse distance weighted (IDW) method based on the spatial distribution     of five rain gauge stations. Due to the lack of flood hydrographs, calibration and validation are performed for the peak discharge of floods at the watershed outlet. The base GSSHA distributed hydrological model of the watershed used in this research includes a square grid containing about 86,200 square cells of 200 m. The Green-Ampt model is used to simulate the infiltration process. Flow routing on the watershed surface and in the main stem of river or in major tributaries is done by a 2-D and 1-D diffusive model, respectively.
The hydrological parameters of the GSSHA model, such as the rates and depths of infiltration and surface roughness coefficients, were determined using geological maps and cluster analysis of satellite images of the study area. The topography of GSSHA was produced using the DEM of the Alos Palsar satellite of Japan with a spatial resolution of 12.5 m.
Based on overall data obtained from field surveys, observations, and satellite images, Manning's n value for all streams used in 1-D diffusive wave routing of model calibration is defined as 0.035. It is very difficult to access all waterways to determine Manning's roughness coefficient, due to the large size of the study watershed, difficult topographical conditions and the lack of access roads; therefore, it is not possible to determine the Manning roughness coefficient separately in each waterway, especially in the scenarios of modeling all waterways with a drainage area greater than 5 km 2 , for 414 waterways (Figs. 10, 11, 12).
Distributed models of surface roughness, soil type, and combined model in study watershed is shown in Figs. 10 through 12, and corresponding mapping tables used in the GSSHA model are presented in Tables 1 and 2. Thus, the base model was calibrated for the two flood events and validated for three other major flood events. The hydrological simulation time was 1800 min (30 h), and the optimal time step in simulating the GSSHA distributed hydrological model was defined as 10 s by trial and error process.  After calibration and validation of the base model, five proposed storms, including 24-h storms uniformly distributed on the watershed with the cumulative depths of 60, 65, 70, 75, and 80 mm, were applied and the flood hydrographs were extracted. The main goal in choosing the uniform spatial distribution of design storms is to compare the topographic, geometric, and physical changes of the base model and perform sensitivity analyses in identical meteorological conditions independent of the storms' spatial variations.
Preliminary results showed that unlike the lumped and semi-distributed hydrological models, the concentration-time of the watershed in flood events with different peak discharges is fully variable, depending on the physiographic and hydrological characteristics of the watershed and the magnitude of design storms. By increasing the design storm's depth, the time to peak discharge may decrease or increase, depending on the watershed's topographic, geometric, and hydrographic complexity.
Sensitivity analysis was performed using the GSSHA calibrated baseline model for the following four scenarios.

Sensitivity analysis based on accuracy of DEM
To evaluate the accuracy of the topographic model used in the GSSHA base model, the Alos Palsar DEM of Japan, with a spatial resolution of 12.5 m, was replaced by the USA Aster GDEM with a spatial resolution of 30 m. The new GSSHA model was performed for five aforementioned proposed storm events, and the results were analyzed and compared with the results of the base model. Results showed that peak discharge of flood events in the base GSSHA model is slightly lower than the corresponding values in the GSSHA model with a spatial resolution of 30 m (Tables 3,  4; Figs. 13, 14). The much smaller volume of the numerical model based on topography with a spatial resolution of 30 m compared to the base model (about five times), these results confirm that using a topographic model with a lower spatial resolution does not have practical errors in the accuracy of the GSSHA model.

Sensitivity analysis based on density of streams
The influence of stream density on determining the peak discharge by GSSHA model was investigated. For this purpose, firstly, a 1-D hydraulic routing of streams was provided for Karvandar River and its major tributaries by the diffusive wave method. Another GSSHA model was provided for the streams and waterways containing the drainage area greater than 5 km 2 . Five storm events were simulated by these two models and their results were compared in Table 5 and Figs. 15 and 16. As deduced from the table and figures, by increasing the density of waterways in the GSSHA model, the peak flood increased significantly, emphasizing the high sensitivity of the GSSHA model to the streams and waterways density.
In another scenario, the simulation was performed without 1-D diffusion wave hydraulic routing, and the results were compared with the results of the base model. The results for design storms with different cumulative depths showed that peak discharge and flood volume for the 1-D hydraulic model in waterways were significantly higher than those performed without 1-D simulations (Table 6; Fig. 17). The differences in results were more significant for the storms with low cumulative depths. For a 60 mm storm depth, the ratio of peak discharge obtained from a model by ignoring hydraulic routing to a model by considering hydraulic routing was 9%, while this ratio was 37% for an 80 mm storm depth.

Sensitivity analysis based on grid size
To evaluate the effect of cell size on the accuracy of results, five different cell sizes were adapted to 100,150,200,250, and 300 m in modeling the geometry of the watershed. Five storm events were performed by each model proper to the corresponding cell size. Results showed that the peak discharges in the three models with grid cell sizes of 100, 150, and 200 m were very close and had considerable difference from the peak discharges with 250 and 300 cell size. These results indicate that using a solution grid with a cell size of more than 200 m   in GSSHA models can lead to significant errors in the results of hydrological distributed simulations. Therefore, it is suggested to limit the solution grids with cell sizes less than 200 m (Tables 7, 8 ,9,10,11,12;Figs. 18,19,20,21,22,23). Sensitivity analysis for simultaneous effects of cell size and stream density was also investigated. A scenario with a100 m cell size was also performed. Results confirmed that by reducing the cell size, the sensitivity of the model to the stream density decreases (Table 13; Fig. 24a, b).
Therefore, using the GSSHA model with a solution grid of smaller cell sizes can improve model calibration. However, due to the increasing number of cells in the network, the simulation time of the distributed hydrological numerical model increases significantly.

Conclusion
In this study, using an innovative approach in geometric modeling of all small and large streams and waterways from upstream to downstream of the study watershed, effects of topographic, geometric, and physical main factors on the hydrologic response of the GSSHA model were analyzed.
Results showed that the peak discharge of flood in the base GSSHA model developed based on the DEM with a spatial resolution of 12.5 m is slightly lower (< 4%) than the corresponding values in the GSSHA model developed based on the DEM with a spatial resolution of 30 m, which confirms that using a topographic model with a lower spatial resolution has no substantial effect on the accuracy of GSSHA model. Also, by increasing the density of waterways in the GSSHA model, the peak discharge of flood increases significantly (44-57%) compared with the corresponding values in the base GSSHA model. Furthermore,   results showed that peak discharge values in the three models with grid cell sizes of 100, 150, and 200 m (base model), especially for intense floods, are close (< 15 and 0.7% difference for 60 and 80 mm storm depths, respectively) and compared with the results of the two models with grid cell sizes of 250 and 300 m show significant and meaningful differences (about 20 and 12% difference for 60-80 mm storm depths, respectively). These results indicate that a solution grid with a cell size of bigger than 200 m in GSSHA model can lead to significant errors in the results. It is recommended to use satellite rain models in future studies and compare their results with each other. Also, compare the results of GSSHA model with the results of other hydrological models.
Funding The research has not been supported through any funds.
Data availability All data generated or used during the study are applicable if requested.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval Not applicable.

Consent to participate
The authors have made a significant contribution to this manuscript, have seen and approved the final manuscript.

Consent to publish
The authors have agreed to publish the study in Applied Water Science Journal.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.