The Performance of Physically Based and Conceptual Hydrologic Models: A Case Study for Makkah Watershed, Saudi Arabia

Abstract

Population growth and land use modification in urban areas require the use of accurate tools for rainfall-runoff modeling, especially where the topography is complex. The recent improvement in the quality and resolution of remotely sensed precipitation satisfies a major need for such tools. A physically-based, fully distributed hydrologic model and a conceptual semi-distributed model, forced by satellite rainfall estimates, were used to simulate flooding events in a very arid, rapidly urbanizing watershed in Saudi Arabia. Observed peak discharge for two flood events was used to compare hydrographs simulated by the two models, one for calibration and one for validation. To further explore the effect of watershed heterogeneity, the hydrographs produced by three implementations of the conceptual were compared against each other and against the output of the physically-based model. The results showed the ability of the distributed models to capture the effect of the complex topography and variability of land use and soils of the watershed. In general, the GSSHA model required less calibration and performed better than HEC-HMS. This study confirms that the semi-distributed HEC-HMS model cannot be used without calibration, while the GSSHA model can be the best option in the case of a lack of data. Although the two models showed good agreement at the calibration point, there were significant differences in the runoff, discharge, and infiltration values at interior points of the watershed.

1. Introduction

Hydrologic models are used to solve a range of specific problems in the management and development of water and land resources, including flood simulation and prediction, aquifer recharge management, runoff estimation, and drainage network design (e.g., [1,2,3]). The uncertainties inherent in these models can be reduced with the availability of accurate input and watershed data [4]. Hydrologic models include lumped models, semi-distributed models, or fully distributed models [5]. These models can also be classified based on the model formulation as empirical or physically-based models with conceptual formulations typically used in lumped models and physical equations in the other two types [6]. The parameters of the lumped models are not directly related to the watershed’s physical characteristics, which is an important shortcoming of these models [7]. Generally, the application of lumped models is limited to gauged watersheds not undergoing significant change in their conditions, and they have to be calibrated using significant observational data [7].

In a semi-distributed model, the basin to be simulated is divided into smaller sub-basins to capture the spatial variability of the inputs and to represent the heterogeneity within the basin [8]. The data required for a semi-distributed model is typically less than that needed for a fully distributed model. In addition, the spatial representation and mathematical formulations of the hydrologic processes are simpler [9]. Although a semi-distributed model requires a limited number of parameters compared to a fully distributed model, it requires more calibration data [10]. In general, the natural spatial and temporal variations of the hydrological properties such as soil type, land use, and rainfall are not accurately represented in the semi-distributed hydrologic models [11].

The fully distributed models try to simulate the watershed conditions accurately while the physical processes and/or the spatial heterogeneity are mathematically treated at the grid cell scale [12,13]. The distributed models can be either empirically distributed models or physically-based distributed models. The difference between them is the transformation of the surface runoff, which is empirically simulated in empirical distributed models and through physical formulations in physically-based distributed models [14]. The model grid size is selected based on the required level of detail of runoff simulations and the spatial resolution of inputs and watershed data [15,16].

Several studies have compared different types of hydrologic models. Paudel et al. [17] compared the performances of lumped and distributed models and concluded that the simulation of land use changes using distributed models is significantly better than using lumped models. Sith & Nadaoka [18] found that the physically-based, fully distributes model GSSHA (Gridded Surface/Subsurface Hydrologic Analysis) performed slightly better than the semi-distributed, semi-physically-based model SWAT (the Soil and Water Assessment Tool) for short-term streamflow simulations, while, for long-term simulations, both models had acceptable performance. Hui-lan [19] concluded that the hydrographs predicted by a distributed model compared with field observations were much better than those predicted by a semi-distributed model. El-Nasr et al. [20] found that the semi-distributed SWAT model and the fully distributed MIKE SHE model performed well when used to simulate the hydrology of a large catchment, but the MIKE SHE model was slightly better. Meselhe et al. [21] investigated the impact of spatial and temporal sampling of rainfall on runoff predictions using a conceptual and physically-based hydrologic model over a small watershed. They found that the physically-based hydrologic model was more sensitive to both the temporal and spatial samplings of rainfall than the conceptual model. Other researchers evaluated the performance of HEC-HMS (the Hydrologic Engineering Center-Hydrologic Modelling System) and GSSHA models and reported a much better performance by GSSHA (e.g., [22,23]).

Several hydrologic studies have been conducted in the Arab Peninsula using different modeling approaches. Al Abdouli et al. [24] used HEC-HMS to quantify coastal runoff and examine its temporal and spatial distribution in several watersheds in the United Arab Emirates (UAE) to assess the potential for recharge of depleted aquifers and flood mitigation in urban areas. Sharif et al. [25] simulated an extreme flood event in Hafr Al Batin city, Saudi Arabia using the physically-based, fully-Distributed Hydrologic model (GSSHA) driven by Global Precipitation Measurement (GPM) satellite rainfall data. Their results demonstrated that the urban portion, which is 6.8% of the watershed area, produced about 85% of the generated runoff. Embaby et al. [26] developed an integrated approach combining geotechnical investigation and hydrologic modeling to generate hazard zones for Tabuk city, Saudi Arabia. The approach was successful in identifying areas prone to flooding. Al-Zahrani et al. [6] developed a flood hazard map for Hafr Al-Batin city by integrating the results of three models (hydrologic, hydraulic, and flood delineation models). They recommended that their result can be used for designing flood control structures and other planning purposes. Fathy et al. [27] evaluated hydrologic models over Wadi Sudr, Sinia, Egypt, and found that the runoff hydrograph can be estimated more accurately using the distributed model. To the authors’ knowledge, no study tried to investigate the differences between conceptual semi-distributed and physically-based fully distributed models in simulating extreme flood events in the region.

In recent years, the Makkah region in western Saudi Arabia has been suffering from an increase in the frequency of flash floods due to its natural complex topography, which is dominated by mountains with steep slopes and very shallow soil layers. These events are triggered by the occurrence of infrequent high convective short-lived rainfall events. A few hydrologic studies have been performed in the Makkah region. However, all studies used lumped or semi-distributed models and were mostly focused on simulating peak discharge. Abdelkarim & Gaber [28] assessed the impact of flash flooding in the Makkah area by constructing a flood risk map for a watershed in Makkah City using HEC-HMS and HEC-RAS. Their findings suggested a flood control policy based on using the existing hydraulic facilities and a range of new structures that can help protect lives and the urban infrastructure. Elfeki et al. [29] evaluated the flood hazard on another watershed in Makkah, also using HEC-HMS and HEC-RAS, and recommended some structural measures to mitigate floods. Bastawesy et al. [30] assessed the potential of flash flooding in a third watershed in Makkah using GIS and recommended a range of mitigation measures to contain the 50-year design flood. Dawod et al. [31] applied simple hydrologic models and four national regression models to estimate the flood peak discharge in several ungauged small watersheds in Makkah. They concluded that the Curve Number (CN) approach is the best methodology for flood estimation in case of the availability of soil type, land use, and metrological data. Al-ghamdi et al. [32] used the Curve Number method to estimate flood hazards in Makkah City for the period 1990 and 2010. They found that the expansion of the residential areas was behind the significant rise in significant flood hazards. Al Saud [33] analyzed the morphometry of Wadi Aurnah, Makkah, using GIS and introduced a simple approach for rapid assessment of the flooding potential.

The main goal of this study is to simulate recent floods in the Makkah region using a physically-based fully distributed hydrologic model and compare its performance to that of the semi-distributed HEC-HMS model. Two flood events that occurred in 2010 and 2018 were used in this study. One event was used for model calibration and the other for validation. The difference in model output was discussed and was further investigated by comparing the different implementations of HEC-HMS to explore the impact of watershed heterogeneity on the model results.

2. Study Area

Makkah province is one of thirteen provinces in Saudi Arabia and is considered one of the most important regions due to its location in the historic Hejaz region and its extended coastline on the Red Sea. The total area of Makkah region is approximately 153,200 km² occupied by around 8,557,000 people [34], making it the third-largest and most populous province. The city of Makkah (Makkah Al-Mukarramah), the Makkah Province’s capital, is located 70 km from Jeddah on the Red Sea with an average elevation of 277 m above sea level. The region is characterized by a dry climate with sparse, inconsistent rainfall events [35]. Plains and mountains are the main geomorphological units of the study area. In the past few decades, the city of Makkah has witnessed unusually rapid urbanization, making it the third most populated metropolitan area in Saudi Arabia. The increase in the urban areas in Makkah province between 1992 and 2013 was threefold, as reported by Alahmadi & Atkinson [36]. This growth is mainly due to the significant increase in the area covered by Makkah City from about 12% of the province in 1992 to 22% in 2016 [37]. Recently, Abdelkarim & Gaber [28] reported that the urban areas exposed to flooding have increased by approximately 25 fold between 1988 and 2019. In addition, a high percentage of Makkah’s road network (more than 50%) is subjected to inundation during major flood events [36]. Because of the rapid urban growth, municipal authorities now require flood impact studies being conducted before any land development activity within the city limits. The Makkah watershed used in this study was delineated using a high-resolution digital elevation model (DEM) of 10 m obtained from King Abdulaziz City of Science and Technology (KACST). The delineated watershed includes most of Makkah city with a drainage area of 1725 km² (Figure 1).

3. Rainfall Events

The calibration and validation run of the two models were driven by two storm events that occurred recently in the Makkah watershed. Due to the very sparse rain gauge network in the study area, capturing the features of the extreme precipitation through ground observations was not possible. Moreover, the very limited rain gauge data are available only at a daily time scale. Therefore, the two hydrologic models were forced using the satellite rainfall data, which would have the ability to capture the spatial and temporal variability of the storm events at much higher resolutions.

A previous study by the authors using the Integrated Multi-satellitE Retrievals for GPM (IMERG) satellite rainfall products showed reasonable performance of the IMERG Early run product in the Makkah region compared to the other IMERG products. Accordingly, the two simulated models were forced by the IMERG Early run product in this study. The IMERG Early run data with spatial and temporal resolutions of 0.1° × 0.1° and 30 min were downloaded using FileZilla software from the Precipitation Measurement Missions (PMM) website (http://pmm.nasa.gov/data-access/downloads/gpm (accessed on 7 March 2020)). The data were converted from HDF5 into a gridded ASCII format using an R script developed by the authors and was then reformatted as input to the hydrologic model. The latest version of the IMERG product, Version 6 (IMERG-V06), was used. The latest product includes significant updates such as a higher maximum rainfall threshold from 50 to 200 mm/h, full inter-calibration to the GPM combined instrument dataset, and the use of an updated rain retrieval algorithm [38].

The first storm that occurred on 13 February 2010 with an average watershed total of 54.7 mm based on the Early IMERG product (Figure 2) was used for calibration. The IMERG Early product showed that the intense rainfall was concentrated over the northern portion for the 13 February 2010 event and the southeastern part of the watershed for the 3 November 2018 event. The second storm that occurred on 3 November 2018 was used for validation. Three rain gauge stations recorded the total rainfall accumulation for this event: Makkah (J114), Arafah (9004), and Muntasaf-Huda (J205), with rainfall totals of 7.8, 22.0, and 20.0 mm, respectively. The watershed-averaged total rainfall observed by the three gauges was 18.36 mm based on the Inverse Distance Squared Weighting method, while it was 42 mm based on Version 6 (Figure 2). The Early IMERG product significantly overestimated the rainfall amount for the 3 November 2018 storm event compared to observations by the ground rain gages. However, the average rainfall estimated by the Early IMERGV06 is closer to what can be inferred from a study that has been done on Wadi Al-Nu’man watershed in Makkah [29].

Rainfall frequency analysis was also used to estimate the rainfall for different return periods to assess the models’ performance. The recorded data from 2006 to 2018 was used to estimate the return period for maximum daily rainfall over the study area using Log-Pearson Type III distribution (LPT III) [39]. Figure 3 shows the frequency distribution as fitted to the LPT III distribution.

Table 1. Maximum 24-h rainfall (mm) for different return periods for the study area based on LPT III distribution.

Return Period	2	5	10	25	50	100
Rainfall (mm)	28	48	64	88	107	129

shows the maximum 24-h rainfall values for the return periods of 2, 5, 10, 25, 50, and 100 years.

4. Models Description

4.1. Semi-Distributed Modeling (HEC-HMS Model)

A semi-distributed model such as HEC-HMS is a compromise between a lumped model and a fully distributed one. In HEC-HMS, the watershed is divided into sub-basins to represent the watershed’s spatial heterogeneity, such as the variation of land cover/use and soil. The outputs (hydrograph, runoff volume) of the model are calculated at each sub-basin outlet and then routed downstream to estimate the hydrograph and runoff volume at the whole watershed outlet. HEC-HMS is arguably the most used model in the United States due to its flexibility, ease of use, and acceptance by almost all governmental entities. It is also one of the most popular models across the globe. The HEC-HMS model has the advantages of high-efficiency simulation, short computation times, and comparatively low calibration needs. However, it lacks the ability to consider localized effects of variation of rainfall, soil, land use, and topography [40]. In addition, the model is incapable of simulating the overland flow even at the subbasin level. Although HEC-HMS requires relatively fewer parameters for calibration, these parameters are not directly related to the watershed’s physical characteristics.

The earliest version of HEC-HMS was originally developed to simulate hydrologic processes of branched (or dendritic) watersheds [41]. The model’s software environment includes a database, computation engine, data entry utilities, and results’ reporting tools. The mathematical model used to represent the integrated watershed processes is based on a user-specified methodology (Figure 4). For example, there are eight precipitation methods to choose from, including gauge weighted average, specified hyetograph, and frequency storm. The most common method used to simulate runoff in HEC-HMS is the Soil Conservation Service (SCS) Curve Number method [11,42]. However, there are also 11 methods to choose from to determine infiltration and runoff volume, including the Curve Number and Green-Ampt.

In this study, the Curve Number (CN) loss model [43] was used to compute surface runoff. According to this method, the excess rainfall is computed according to

$$Q_r=\frac{{\left(P-0.2S\right)}^2}{P+0.8 S}$$

(1)

where Q_r is the direct runoff (excess rainfall); P is the rainfall, and S (storage in the watershed) is related to the Curve Number (CN), which is computed by the following equation:

$$S=\frac{25400}{CN}-254$$

(2)

The CN values in Equation (2) range from 30 for high infiltration soils to 100 for impervious surfaces where the runoff is practically equal to rainfall. The SCS Unit Hydrograph method [44] was used to transform the excess precipitation into a runoff hydrograph using the following equations, which estimate the peak discharge (Q_p) and the time to the peak discharge (t_p):

$$Q_p=2.08 \frac{A}{Tp}$$

(3)

$$T_p = \frac{\Delta t}{2}+t_{lag}$$

(4)

where Q_p = peak discharge (m³/s), T_p = the rise time of the hydrograph (hours), A = watershed area (km²), Δt = the excess duration of precipitation (hours), and t_lag = lag time (hours).

The Muskingum routing method was selected in this study, and the following equation is solved by this method:

$$Q_t=\left(\frac{\Delta t-2KX}{\Delta t+2K\left(1-X\right)}\right)I_t+\left(\frac{\Delta t-2KX}{\Delta t+2K\left(1-X\right)}\right)I_{t-1}+\left(\frac{2K\left(1-X\right)-\Delta t}{\Delta t+2K\left(1-X\right)}\right)Q_{t-1}$$

(5)

where O_t = the outflow; t = time; I = the inflow; K is the storage time constant for the reach, and X is weighing factor.

4.2. The Physically-Based, Fully Distributed Hydrologic Model (GSSHA)

Figure 5 shows the GSSHA model set up. The GSSHA model accepts spatially variable inputs, making it very useful for evaluating land-use change effects, modeling ungauged watersheds when there is no data for calibration, and simulating the water quality and sediment yield [17,45]. In the GSSHA model, the basin is divided into grid cells considering the spatial variability of the watershed properties. The flow in every grid point can be estimated using physical equations. Accordingly, detailed hydrologic predictions at the grid-cell level across a basin can be provided by this approach. The GSSHA is a process-based model that is used to simulate channel flow (one dimension) and overland flow (two dimensions) using finite-volume solutions of the St. Venant equations of motion [46,47], and all the processes are simulated over each grid cell. The major components of the GSSHA model are: temporally and spatially varying precipitation, evapotranspiration, simplified lake storage and routing, surface runoff routing, an interception by plants, infiltration, wetland peat layer hydraulics, saturated groundwater flow, unsaturated zone soil moisture accounting, overland sediment erosion transport and deposition, overland contaminant transport and uptake, and in-stream sediment transport [48].

In GSSHA, the overland flow routing can be computed using one of three numerical techniques: ADE-Prediction Correction (PC), Explicit, and Alternative Direction Explicit (ADE). The characteristics or type of watershed will control the selection of the appropriate numerical technique. The following equation represents the ADE scheme used in this study. The inter-cell flows are calculated using Equation (6) along the x-direction:

$$P_{ij}^N=\frac{1}{n}{\left(h_{ij}^N\right)}^{\frac{5}{3}}{\left(S_{fx}^N\right)}^{\frac{1}{2}}$$

(6)

Equation (7) is used to calculate the flow depths in each cell at the n + 1 time level based on the flows along the x-direction:

$$h_{ij}^{\frac{N+1}{2}}=h_{ij}^N+\frac{\Delta t}{\Delta x}$$

(7)

Interflows can be calculated from each cell using Equation (8) along the y-direction:

$$q_{ij}^{\frac{N+1}{2}}=\frac{1}{n}{\left(h_{ij}^{\frac{N+1}{2}}\right)}^{\frac{5}{3}}{\left(S_{fy}^{\frac{N+1}{2}}\right)}^{\frac{1}{2}}$$

(8)

The updated column depths are calculated using Equation (9) along the y-direction based on the interflows:

$$h_{ij}^{N+1}=h_{ij}^{\frac{N+1}{2}}+\frac{\Delta t}{\Delta x} \left(q_{i,j-1}^{\frac{N+1}{2}}-q_{ij}^{\frac{N+1}{2}}\right)$$

(9)

where: $p_{ij} and q_{ij}$—overland flows from cell ij in the x and y directions, respectively, $h_{ij}$—the water depth in cell ij at the nth time level, n—roughness coefficient, $S_{fx} and S_{fy}$—the slopes of the water surface in the x and y directions, respectively, and $\Delta x=\Delta y$—cell’s dimensions.

The head discharge is computed based on Manning’s equation to route the channel flow Equation (10):

$$Q_{i+1}^n=\frac{1}{n} A_i^n {\left(R_i^n\right)}^{\frac{2}{3}} {\left(S_{f_{i+\frac{1}{2}}}^n\right)}^{\frac{1}{2}}$$

(10)

where $Q_{i+1}^n$ = intercell flow, n = channel roughness, S_f = the friction slope in the x-direction, A = cross section area (m²), and R = hydraulic radius. Equation (11) is used to calculate S_f.

$$S_{f_{i+\frac{1}{2}}}^n=S_{o_{i+\frac{1}{2}}}^n-\frac{h_{i+1}^n-h_i^n}{\Delta x}$$

(11)

If the effects of the backward flow are considered, the flow is calculated using the hydraulic head in the downstream cell as follows:

$$Q_{i+\frac{1}{2}}^n=\frac{1}{n} A_{i+1}^n {\left(R_{i+1}^n\right)}^{\frac{2}{3}} {\left(S_{f_{i+\frac{1}{2}}}^n\right)}^{\frac{1}{2}}$$

(12)

The volume of flow is calculated at each node using Equation (13):

$$V_{ij}^{1+N}=V_i^N+\Delta t\left(q_{lat}^{N+1}\Delta x+q_{rec}^{N+1}\Delta x-Q_{i-\frac{1}{2}}^N-Q_{i+\frac{1}{2}}^N\right)$$

(13)

where:

• $q_{lat}$ = lateral flow (m²/s)
• $q_{lat}$ = flow exchange between the channel and groundwater (m²/s)
• $Q_{i-\frac{1}{2} }^{N}and Q_{i+\frac{1}{2}}^N$ = the computed inter-cell flows from depths (h) in the longitudinal (x), at the Nth time level.

The infiltration is calculated using the GAR method as follows:

$$f (t)=K\left(\frac{S_f \left({\theta }_s-{\theta }_i\right)}{F \left(t\right)}+1\right)$$

(14)

where:

• f (t) = Infiltration Rate (cm/hr)
• S_f = wetting front suction head (cm)
• K = effective hydraulic conductivity (cm/hr)
• θ_i = initial soil moisture content
• θ_s = water content of the soil at natural saturation.

5. Watershed Data

The primary input data required for the GSSHA and HEC-HMS models include rainfall data, digital elevation models (DEM), soil types, and land use/cover, which need some processing and preparation in the right input format. The DEM data, at 10 m resolution, were obtained from King Abdulaziz City of Science and Technology (KACST). The soil type data were downloaded from the SoilGrids^TM global digital soil mapping system (www.SoilGrids.org (accessed on 20 March 2020)) as mass fractions of clay, silt and sand in percentages. The final soil type map (Figure 6) was developed for the Makkah watershed using ArcGIS10.6 [49], and the soil texture classification based on the clay, silt, and sand contents and then adjusted after comparison with aerial photographs and field information. The land use/cover data (Figure 6) were obtained from the OpenLandMap data portal (www.openlandmap.org (accessed on 25 March 2020)), and adjustments were made after comparison with aerial photographs. The hydrologic soil group data were obtained from the HYSOGs250m product ([50] https://daac.ornl.gov/SOILS/guides/Global_Hydrologic_Soil_Group.html (accessed on 26 March 2020)) at 250 m resolution as shown in Figure 7a. This data were used to estimate the CN values (Figure 7b) needed for the HEC-HMS model.

6. Methodology

The methodology adopted in this study is described through the flowchart shown in Figure 8. ArcGIS 10.6 was used to prepare the raster data for WMS [51] to be processed and exported to HEC-HMS. The process in WMS includes delineation of the watershed, flow accumulations, flow directions, and computing CN values based on the hydrologic groups in Figure 7. The watershed was subdivided into 55 sub-basins for HEC-HMS simulations, as shown in Figure 9a. The SCS Curve Number and SCS dimensionless transformation methods were applied to compute the runoff volumes and to model direct runoff, respectively, while the Muskingum method was used for channel routing. The inverse distance method was used to generate the rainfall hyetograph for the selected storm events before being used as input to the HEC-HMS model.

The channel roughness, channel cross-sections, and other parameters such as hydraulic conductivity were required for the GSSHA model. ArcGIS and WMS were also used to prepare the raster data for the GSSHA model. All streams were smoothed to avoid the effect of the adverse and flat slopes resulting from inaccuracies in the DEM. The watershed was divided into grid cells at 150 × 150 m resolution (Figure 9b). The land use/cover and soil types were imported by WMS to prepare the grid index maps which then were used to estimate the required parameters such as infiltration parameters, roughness, initial moisture, overland retention depth, and impervious fraction. The WMS was also used to extract channel cross-sections from the 10 m DEM (Figure 10). The no-flow condition was adopted for all first-order streams boundary conditions in GSSHA as there is no baseflow in the basin [47].

7. Results

7.1. Models Calibration

The storm event of 13 February 2010 with total rainfall accumulation of 55 mm resulting in a peak discharge of 431 m³/s and a flow depth of 1.6m measured at a point of intersection between Wadi Uranah (610.8 km²) and Mashar Arafat [52] was used to calibrate both the GSSHA and HEC-HMS models. Since HEC-HMS is not a fully distributed model, each IMERG rainfall grid value was input to the HEC-HMS as point observation at the grid center, and this process was applied over the whole watershed. The rainfall hyetograph for each sub-basin for the selected storm event was generated using the Thiessen Polygons method. On the other side, GSSHA can accept temporally and spatially varying precipitation from a satellite. A time step of 15-min was used for the discharge output for the two models. Figure 11 shows the simulated hydrographs from HEC-HMS and GSSHA simulations before calibration. It can be noticed that HEC-HMS significantly overestimated the observed peak discharge, while the GSSHA model slightly underestimated the observed peak discharge. Accordingly, the HEC-MS was calibrated by adjusting the CN values and Muskingum parameters and the GSSHA model by adjusting the hydraulic conductivity and roughness values.

Both models produced acceptable hydrographs after calibration, as shown in Figure 12. However, the runoff volume produced by HEC-HMS after calibration was smaller by 38% than that produced by GSSHA. In other words, the percentage of losses due to infiltration that resulted from HEC-HMS (73.6%) is more than that for GSSHA (57.5%). We could not verify the two models’ runoff volume since the observed runoff volume was not available. It can be seen from Figure 12 that the base time of the hydrograph simulated by GSSHA is much greater than that for HEC-HMS. This could be attributed to the significant effect of the grid cells in the distributed model, representing the spatial variability of the watershed properties more accurately. On the other hand, the basin is divided into sub-basins in the HEC-HMS model, and the surface slope is taken as the average for each sub-basin. This cannot accurately represent the real conditions, especially in complex terrain dominated by mountains with steep slopes such as Makkah watershed, where the average slope masks the effect of relief variability.

7.2. Models’ Validation

The HEC-HMS and GSSHA models were validated using the storm event of 3 November 2018 when a peak discharge of 640–680 m³/s was estimated at three points in the main channel near the outlet of Wadi Al-Nu’man, which has a drainage area of 683.4 km² (Figure 1). These estimates were based on a maximum flow depth of about 3.5 to 4 m measured by observers from the University of Umm Al-Qura and a consulting company [28]. The observed runoff hydrograph could not be reconstructed because the flow depth and peak discharge values were the only available observed data. Moreover, there is no information about the peak discharge timing and hydrograph duration. Accordingly, we used one statistic measure in model validation, which is the error in peak discharge (ε_p) using the following equation:

$${\epsilon }_p (\%)=\frac{\left(P_o-P_s\right)}{P_o}\times 100$$

(15)

where: P_o is the observed peak discharge; P_s is the simulated peak discharge.

When using the IMERG Early run product as input, HEC-HMS overestimated the peak by 4.5%, while GSSHA overestimated the peak by just under 1% (

Table 2. Comparison between GSSHA and HEC-HMS simulation results.

Model	Peak Discharge (m3/s)			Runoff Volume (m3)			Error in Peak Discharge (%)
	13 February 2010		3 November 2018	13 February 2010		3 November 2018	13 February 2010		3 November 2018
	Before Calibration	After Calibration	Validation	Before Calibration	After Calibration	Validation	Before Calibration	After Calibration	Validation
GSSHA	393.2	427.6	663.6	12,866,921	14,311,943.5	19,455,825.4	8.8	−0.8	0.6
HEC-HMS	847.9	429.5	689.8	17,150,300	8,882,200	11,933,600	96.7	0.3	4.5

and Figure 13). However, the estimated runoff volume and the percentage of infiltration losses estimated by the two models were quite different. Again, the runoff volume estimated by HEC-HMS is much lower than the GSSHA estimate. This may indicate that HEC-HMS was not able to represent the complex topography of the study area.

7.3. Models’ Outputs Comparison

Figure 14 and Figure 15 show the hydrographs simulated by the two models at the entire Makkah watershed outlet. The outlet peak discharge values produced by HEC-HMS for calibration and validation events were lower than those produced by the GSSHA model by 16.3% and 36.7%, respectively. For the 2010 and 2018 storm events, the HEC-HMS model estimated runoff ratios of 21.2% and 18.9%, respectively, while the runoff ratios calculated by GSSHA were 39% and 41.5%, respectively.

Unlike HEC-HMS, GSSHA provides distributed outputs of major hydrologic states and fluxes such as infiltration, runoff at high temporal and spatial resolutions, and soil moisture as maps. However, each sub-basin’s infiltration values can be estimated from HEC-HMS sub-basin results. These values can be used to develop average infiltration depth maps using GIS. Figure 16 shows maps of the cumulative infiltration depths (m) produced by HEC-HMS and GSSHA. The figure highlights the significant difference between the estimates of the two models. In HEC-HMS, the average infiltration is calculated by considering the whole sub-basin as one unit, while in GSSHA, the infiltration is calculated for each grid cell. As a result, the uniform infiltration over each sub-basin is not realistic. The GSSHA infiltration map shows the infiltration in limited areas, around the main channels, and in the flat areas with the soil of high hydraulic conductivity, which could be closer to the real conditions. In addition, the GSSHA infiltration map shows very low infiltration depth values over the mountains because of the low hydraulic conductivity and steep slopes. In contrast, the HEC-HMS map shows higher infiltration depths over these areas.

The GSSHA model has the ability to estimate all model outputs at any location inside the basin. To compare the two models’ outputs at the sub-basin scale, we allow GSSHA to estimate output hydrographs at points that coincide with some selected HEC-HMS sub-basins outlets. The locations of the points selected for comparison are shown in Figure 16. The comparison statistics are shown in

Table 3. Comparison between GSSHA as a reference and HEC-HMS outputs at selected locations.

Basin	Model	Peak Discharge (m3)	Runoff Volume (m3)	Differences (%)
				Peak Discharge	Runoff Volume
80B	HEC-HMS	8.1	78,300	78.5	92.9
	GSSHA	37.6	1,097,316
81B	HEC-HMS	9.1	90,100	78	90.8
	GSSHA	41.3	978,141
127B	HEC-HMS	23.4	335,400	59.2	77.7
	GSSHA	57.3	1,505,590
120B	HEC-HMS	135.4	1,882,000	26.9	59.5
	GSSHA	185.1	4,646,771
122B	HEC-HMS	131	1,907,500	41.8	63
	GSSHA	225.1	5,152,439

. The table shows that the peak discharges estimated by HEC-HMS are 26.9 to 78.5% lower than the GSSHA estimates, while the runoff volumes are lower by 59.5 to 92.9%. This significant difference between the two models can be attributed to several reasons. For example, HEC-HMS produced high infiltration values for the sub-basins 80B and 81B, located in the upper right portion of the basin, whereas the GSSHA model produced very low infiltration resulting in huge differences in the peak discharge and runoff volume, as shown in Table 3. Sub-basin 120B includes one of the widest channels underlain by highly permeable soils. For this sub-basin, GSSHA estimated high infiltration values mostly in and around the channels, as shown in Figure 16, while the whole sub-basin was subjected to high infiltration as shown in the HEC-HMS map. This comparison confirms that the calibration and validation at one location do not lead to an agreement of the results of the two models over the entire watershed, especially for large watersheds. As a result, several locations for model calibration and validation are required. However, in case of the difficulties of getting the required data, the calibration/validation at the location of the existing/suggested hydraulic structure will still be valuable.

A major advantage of the GSSHA model is its capability to simulate infiltration for small subsegments of the watershed. This can help to identify areas within the watershed with high infiltration and consequently can help to determine locations within the watershed where excess runoff for managed aquifer recharge and prevent flooding without costly structural measures can be utilized. Figure 17 shows the distributed runoff depths. The detailed information on runoff depths produced by the GSSHA model can be utilized to guide innovative flood mitigation and countermeasures in the urbanized portions of the watershed and even traffic planning and rescue efforts during extreme events.

7.4. Effect of Topography on HEC-HMS Simulations

To gain more understanding of the effect of complex topography and spatial heterogeneity of the study area on the model predictions, GSSHA and three models of HEC-HMS were used to estimate the outlet hydrograph for the whole Makkah watershed. The ModClark approach, one of many rainfall-runoff methods in HEC-HMS, provides an advanced feature for gridded runoff simulation. The ModClark approach is based on Clark’s Unit Hydrograph [53]. It employs a grid representation that makes it convenient for applying a gridded CN method. This helps examine the effect of the study area’s complex topography on the output hydrograph using the same input data of the semi-distributed models. Two different versions of the semi-distributed model were also used. One is based on the SCS unit hydrograph (SCS-UH), and the other one is based on the Clark Unit Hydrograph (CU-H). The basin lag time is the main parameter of SCS-UH, while the time of concentration (t_c) and a storage coefficient (R) are the main parameters of C-UH. The Muskingum method was used for channel routing for all HEC-HMS models. The three models (semi-distributed (SCS-UH), semi-distributed (CUH), and ModClark) were prepared with the same input data, while GSSHA was prepared with the same minimum required input data. Since the objective was only to show the effect of watershed heterogeneity, Type II storm distribution with a return period of 10 years was used as input to the four rainfall-runoff models. Figure 18 shows the gridded CN map developed for the ModClark model, which is significantly different from the sub-basins based CN map (see Figure 7).

Table 4. Modeling approach.

Model	Direct Runoff/Overland Flow Routing	Runoff Volume Computation/Infiltration	Channel Flow Routing
HEC-HMS semi-distributed (SCS-UH)	SCS unit hydrograph	SCS curve number	Muskingum
HEC-HMS semi-distributed(C-UH)	Clark unit hydrograph	SCS curve number	Muskingum
HEC-HMS ModClark	Clark unit hydrograph	Gridded SCS curve number	Muskingum
GSSHA	Two-dimensional diffusive wave routing	Green and Ampt with redistribution (GAR)	One-dimensional diffusive wave

summarizes the modeling approach used.

It can be noticed from Figure 19 that the hydrograph base time produced by GSSHA and HEC-HMS ModClark are much greater than those produced by the semi-distributed models with the (SCS-UH) model estimating the shortest base time. In addition, it can be seen in Figure 19 and

Table 5. Comparison between GSSHA and HEC-HMS (semi-distributed and ModClark) models.

Model Type	Rainfall (m3)	Runoff Volume (m3)	Runoff (%)	Losses (m3)	Losses (%)	Peak Discharge (m3/s)
GSSHA	94,788,583	37,622,799	39.7	53,142,119	56.1	858.6
HEC-HMS ModClark	94,963,200	40,485,600	42.6	54,241,500	57.1	521.5
HEC-HMS (Semi-distributed using SCS transform method)	94,894,200	39,013,400	41.1	55,880,800	58.9	1536.1
HEC-HMS (Semi-distributed using Clark transform method)	94,894,200	39,013,400	41.1	55,880,800	58.9	1026.3

that the peak discharge produced by GSSHA is greater than that for HEC-HMS ModClark and less than that produced by other HEC-HMS models. The difference in the peak discharge between the semi-distributed (SCS-UH) and (C-UH) models can be attributed to the difference in the unit hydrograph assumptions.

With the same input data, loss method, transform method, and channel routing, the semi-distributed (C-UH) and ModClark models are almost identical except in that ModClark is a distributed model [54]. Thus, the difference in results can be primarily attributed to the study area’s complex topography. The difference in the peak discharge between the semi-distributed (C-UH) and ModClark models (Figure 19 and Table 5) is due to ModClark discretization of the watershed into smaller grid cells, resulting in an increase in the travel time reflecting the effect of the complex topography of the study area. This demonstrates the impact of the study area’s complex topography that can not be represented very well by the lumped or semi-distributed HEC-HMS models. However, the GSSHA model produced higher peak discharge than the ModClark model, which could be attributed to the difference in the modeling methods (Table 5). Halwatura & Najim [55] found that the performance of the SCS-CN was not satisfactory in many cases. Walega et al. [56] found that the prediction of direct runoff by modified Sahu-Mishra-Eldho—curve number (SME-CN) is much more accurate than the original SCS-CN method implemented in HEC-HMS.

8. Conclusions

In this study, a physically-based, fully distributed model (GSSHA) and a conceptual semi-distributed model (HEC-HMS) were used to simulate flood events over an arid watershed in Makkah, Saudi Arabia. The steep slopes and large fraction of impervious surfaces of the study area increase drainages efficiency and contribute to the occasional flash floods in the city. Hydrographs simulated by the two models were compared to observed peak discharge for one flood event using satellite rainfall as input. One more event was used for model validation. The two models were prepared with the minimum required data so that the comparison between the two models would be meaningful. The performance of the models was measured by the peaks’ discharge. The effort required to build the two models was comparable since the same input data were used, and a similar number of parameters were calibrated. However, additional effort was needed to define channel cross-sections for the GSSHA model. Moreover, cell-to-cell surface runoff computations in GSSHA required more time than did the lower resolution sub-basin runoff computations of HEC-HMS.

Although the GSSHA model performed better than HEC-HMS, which significantly overestimated the peak discharge before calibration, the two models’ performance after calibration was comparable. However, the runoff volumes produced by HEC-HMS were much smaller than GSSHA for all calibration and validation events. The failure of HEC-HMS to capture details of the complex terrain and land surface heterogeneity was likely behind this discrepancy. The cell-to-cell routing of flow that was only possible with GSSHA. The GSSHA model can make use of higher-resolution DEMs and rainfall products that are becoming more available and are more suitable in urban areas, where the natural landscape is significantly modified.

The difference in the output of GSSHA and HEC-HMS was investigated further by comparing GSSHA hydrograph to those produced using three implementations of HEC-HMS (semi-distributed (SCS-UH), semi-distributed (C-UH), and ModClark). A uniform design storm was used as input. The two semi-distributed (SCS-UH and C-UH) models were implemented with the same input data, losses method, and routing. The results showed that the semi-distributed (C-UH) model produced less peak discharge than the semi-distributed (SCS-UH) one. This could be related to the effect of the storage coefficient in the C-UH model, which is geomorphological data and unique for each watershed. The HEC-HMS ModClark model hydrograph was comparable to that of GSSHA, indicating that it was able to capture the effect of the complex topography of the watershed. The comparison results were intended only for situations when a quick and simple hydrologic model simulation is to be performed. All models could have been built more deliberately by incorporating more of the watershed features, such as the detention/retention structure in residential areas, if observed runoff data were available. This can be investigated in more detail in future studies.

A major advantage of the GSSHA model is that the detailed information on high infiltration areas produced by the model can help identify the locations where the excess runoff can be utilized in aquifer recharge. This can help reduce flooding without the need for costly structural measures. Moreover, the detailed information on the runoff depths produced by the GSSHA can be used for designing flood control structures and other planning purposes.