Urban flood susceptibility evaluation and prediction during 2010–2030 in the southern watersheds of Mashhad city, Iran

Urban flood susceptibility evaluation (FSE) can utilize empirical and rational procedures to focus on the urban flood evaluation using physical coefficients and land-use change ratios. The main aim of the present paper was to evaluate a flood susceptibility model in the southern watersheds of Mashhad city, in Iran, for 2010, 2020, and 2030. The construction of the model depended on the utilization of some global datasets to estimate the runoff coefficients of the watersheds, peak flood discharges, and flood susceptibility evaluations. Based on the climatic precipitation and urban sprawl variation, our results revealed the mean values of the runoff coefficient (Cr) from 0.50 (2010) to 0.65 (2030), where the highest values of Cr (> 0.70) belonged to the watersheds with real estate cover, soil unit of the Mollisols, and the slope ranges over 5–15%. The averagely cumulative flood discharges were estimated from 2.04 m3/s (2010) to 5.76 m3/s (2030), revealing an increase of the flood susceptibility equal 3.2 times with at least requirement of an outlet cross-section by  > 46 m2 in 2030. The ROC curves for the model validity explained AUC values averagely over 0.8, exposing the very good performance of the model and excellent sensitivity.


Introduction
Recent scientific publications support the fact that the occurrences of climate-related disasters (e.g., flash floods and severe storms) have significantly increased under the abrupt changes through the hydro-meteorological conditions Ahmadi and Moradkhani 2019;Mihu-Pintilie et al. 2019). On the other side, many researchers have reported that increasing urbanization can induce higher flash flood occurrences in urban areas (Ahmad et al. 2017;Scionti et al. 2018;Shah et al. 2020).
According to the literature review, the flood events in urban areas notably have increased trends influenced by climate change and urban sprawl expansion (e.g., Ntelekos et al. 2010;Huong and Pathirana 2013;Li et al. 2013;Zhao et al. 2019). An important point is a balance between urban sprawl development and climate change effects through the future flash flood susceptibility (Dawson et al. 2011). Urban growth can be responsible for over 50% of the increase in flood risk, and urban areas face global challenges of flood susceptibility due to climate change and the development of residential areas and urban sprawl (Berndtsson et al. 2019).
In this regard, urbanization impacts on runoff measures and flood damages are important fields of research (Isidoro et al. 2012), especially in the non-western developing countries. For instance, the study of Abass et al. (2020) explained a direct relationship between urban sprawl and flood occurrence using geospatial techniques and field observations in Ghana. Similarly, the hydraulic modeling of urban floods by Devi et al. (2019) indicated the average increase in inundation extent considering the urban sprawl and extreme rainfall event in India.
Several models have been used to approximate the urban flood modeling on several local and national scales regarding the previous works. For instance, Tierolf et al. (2021) have used a flood assessment based on land-use change characteristics in the present and future periods, and they revealed an overall increase in urban flood risk in Southeast Asian countries. Using a hydrological model system (HMS), Feng et al. (2021) reported that areas influenced by a flash flood and floodplain increase due to the urbanization process in Canada. Also, in other researches, Dastorani et al. (2010) and Hoseini et al. (2017) have analyzed the flood-discharge using watershed modeling systems (WMS) and revealed suitable results adopting with local empirical findings. Furthermore, some studies using land-use classification and hydrologic models (e.g., Kulkarni et al. 2014;Zope et al. 2016;Wagner et al. 2016) revealed the increased runoff and flood susceptibility for Indian urban areas in the future.
It is not easy to quantify the urbanization impact on flooding since different urbanization scenarios result in different outcomes (Du et al. 2012;Fletcher et al. 2013). Overall, flood susceptibility increases significantly by transitioning from natural areas to urban sprawl fields (Prosdocimi et al. 2015;Miller and Hutchins 2017). In this regard, more statistical studies should be carried out to show the long-term quantitative effects of urban development in the flood susceptibility evaluation (Berndtsson et al. 2019). Urban flood susceptibility evaluation (FSE) and management is a proper way to mitigate the urban human losses and economic damages (Tingsanchali 2012), which utilizes the empirical and rational procedures to focus on the urban flood evaluation using physical coefficients and land-use change ratios (Su 2017).
Iran is highly disaster-prone, suffering from droughts, floods, earthquakes, landslides, and anthropogenic and technological disasters. In the semi-arid regions of Iran, flash floods are categorized as a dangerous type of hydroclimatic hazard. For example, in March 2019, Iran experienced several major waves of storms and heavy rainfall within some weeks, demonstrating anomalous events in the past 100 years. According to government official reports, the flash floods have burst about 140 watershed down-streams and inundated about 200 urban areas (UNCTI 2019). Accordingly, the severe flooding of 2019 in Iran was a great momentum for the urban decision-makers and dwellers to consider the forthcoming urban flood hazards in the pre-urban and urban watersheds. Now, the main aim of the present paper is to evaluate a flood susceptibility model in the southern watersheds of Mashhad city, in Iran, during three temporal windows of 2010, 2020, and 2030. The construction of the model depends on the utilization of some global remotely sensed (RS) and reanalyzed datasets through geographical information system (GIS) to evaluate spatial and temporal levels of runoff variations, peak flood discharges, and flood hazard susceptibilities toward the urban sprawl and land-cover changes. The integration of both GIS and RS presents an exceptional technique in monitoring and predicting urban growth and urban floods (e.g., Park et al. 2011;Al-Ghamdi et al. 2012).

Study area
The area under investigation includes about 22 pre-urban and urban watersheds in the southern Mashhad city, Iran, which located between 36°15′ and 36°20′ N latitudes and 59°26′-59°34′ E longitudes (Fig. 1). Based on the urbanized elements, the study area is located between two main orchards of the Mashhad, namely Vakil-Abad and Malek-Abad gardens. According to the municipal administrative boundaries in 2020, the study area with a population rate of 340,000 inhabitants as 10.8% from the total Mashhad city (3,128,000) and a surface area of 43.3 km 2 as 12.7% from the entire city (340 km 2 ), belongs to the Region 9. This region is located over the southern mountains of the Mashhad, including the natural relief and ecological and recreational landscapes. The study area has a semi-arid climate with a mean annual temperature of 14 °C and annual precipitation of 260 mm based on a long-term database (Yazd et al. 2019).
Some of the study watersheds (e.g., watersheds W11, W12, W13, and W14) dominantly are the natural upstream areas, which urban residential sprawls and real estates have not degraded. However, the downstream parts of the study area located on the watersheds W17, W21, and W22, wholly comprise the urban real estate. The focal outlets of all watersheds in the W10 are flowed across the Ferdowsi University zone and are effluent in the main channel into the Malek-Abad Garden, draining stream flows to the Kashafrud River in the northern part of the city.
During the several stages of urban sprawl development, the natural parts of the study area have suffered from urban sprawl growth, especially along the several main roads such as Piruzi, Fakuri, and Namaz highways. In recent years, the construction of the southern highway has destroyed the natural patterns of the study area, inducing the sprawl diffusion into the natural upstream watersheds and variation of the magnitude and frequency of the flooding occurrences. Hence, it may be anticipated that the expansion of urban real estate in the study area will influence the risk of future flooding events within the next decade.

Data preparation
Data preparation of this study follows the research model of the present study (Fig. 2) and needs to consider at least five global datasets. In this regard, the remotely sensed and reanalyzed data were gathered from spatial grid pixels or time series, focusing on the study area's geographical coordination (i.e., 36° N and 59° E).
The procedure for detecting and identifying the datasets presented here is based on Rabbani et al. (2021). Land characteristics of the watersheds such as area, length, height, and slope were obtained from georeferenced tagged image files of the global digital elevation model derived (DEM) via https:// lpdaac. usgs. gov/ produ cts/ astgt mv003 (from ASTER images version 003, with a horizontal resolution of 30 m, released from 2013). The soil units of the study area were extracted from the global soil-grid dataset (with a pixel resolution of 250 m, released from 2014) via https:// soilg rids. org, exposing the soil properties and classes (Hengl et al. 2014(Hengl et al. , 2017. Land-use types and the land-cover dataset were considered from the global land-use/ land-cover (LULC) database via https:// lpdaac. usgs. gov/ produ cts/ mcd12 q1v006 (from MODIS images with a pixel resolution of 500 m, version 006, yearly) retrieved from satellite products in 2010 (Li et al. 2017).
The spatial expansion of the urban sprawl was estimated in the past (2010), present (2020), and the probable future (2030) based on the global dataset of human builtup and settlement extent (HBASE) via https:// sedac. ciesin. colum bia. edu/ data/ set/ uland sat-hbase-v1/ datadownl oad, retrieved from Landsat imageries (available at 30-250 m resolution, version 001 based on the imageries in 2010) . Ultimately, the time-series of daily-precipitation rate (data with 0.1-degree resolution, from 1990 to 2021, FLDAS_NOAH01_C_GL_MA_v001) data was collected from the fourth version of the geospatial interactive online visualization and analysis infrastructure (GIOVANNI) program via https:// giova nni. gsfc. nasa. gov/ giova nni, maintained at NASA Goddard services center.

Research model
The FSE model, presented in this research (Fig. 2), includes some steps. First, we acquired several reanalyzed data (e.g., from DEM, soil grids, LULC, HBASE, and Giovanni) for producing the GIS-based spatial analysis or statistical time-series within 2010-2030 concerning the land and climatic characteristics (e.g., elevation, slope, soil units, land covers, urban sprawl, and precipitation rates). Second, the quantitative procedures were assumed to estimate the key hydrological factors of watersheds such as runoff coefficients (Cr), precipitation intensities (Pi), peak flood discharges (Qp), concentration times (Tc), and outlet cross-sections (Sc) as are described in following equations.
The runoff and flood models, such as rational methods and empirical equations defined by Thompson 2006, Devi et al. (2019, and Cheah et al. (2019), could aid in the estimation of the number of precipitation rates, land characteristics, and some various physical parameters of the watersheds. Hence, the runoff coefficient (Cr) is the ratio that depends on factors such as the intensity of soil units, the types of land-use/landcover, precipitation rates, and slope ranges (Mousavi et al. 2019). According to Table 1, the runoff coefficient of different watersheds can be determined by overlapping the aforementioned factors (e.g., Thompson 2006, Mousavi et al. 2019. Accurate determination of the peak flood discharge in each watershed has the main role in managing water resources (Bhadra et al. 2008;Hoseini et al. 2017). According to Parak and Pegram (2006), the peak flood discharge due to an anomalous rain fall event on a watershed can be determined from the rational formula as expressed below equation: where Qp is the flood peak discharge in cubic meters per second (m 3 /s), Cr is the runoff coefficient (unitless), Pi is the maximum precipitation intensity in millimeters per hour (mm/h) that is determined based on the time-series linear trend, and A is the watershed area in squared kilometers (km 2 ).
The time of concentration (Tc) for watershed flow is computed using Kirpich (1940) formula assigned by Thompson (2006), which has been calibrated using the  where Tc is the time of concentration for the watershed in the hour (h), L is the length of the longest watercourse (km), and S is the slope of the longest watercourse (km/ km). Furthermore, we can use an empirical equation to estimate the engineering dimensions of watershed outlets, floodways, and cross-sectional designs. On this basis and owing to the main factors of time of concentration (Tc), peak flood discharges (Qp), and the length of watersheds (L), the outlet cross-section (Sc) of each watershed is estimated to facilitate the hydraulic design and control using below equation: where Sc is the outlet cross-section (m 2 ), Qp is the flood peak discharge (m 3 /s), Tc is the time of concentration for the watershed (s), and L is the length of the longest watercourse (m). In this step, the field observations were carried out in the outlet points of the watersheds to control the estimated results compared with the quo status (2020). Third, the flood susceptibility evaluation of the watersheds was approximated as the final output of the FSE model. In this regard, FSE of the watersheds within three periods of 2010, 2020, and 2030 were defined as the cumulative flood discharges of each watershed after the concentration times considering the upstream watershed discharges and outlet cross-sections. Ultimately, the model validity was computed based on the receiver operating characteristic (ROC) procedure.
The ROC curve, a widely used informative area under the curve (AUC), was plotted to validate the flood susceptibility model in the present study regarding the urban sprawl effects. AUC ranges from 0.5 to 1, where the AUC values above 0.8 represent the very good performance of the model and excellent sensitivity (Yesilnacar and Topal 2005). The overall performance indicator of a model is drawn by plotting false positive (FP) rate or 1-specificity on the x-axis against true positive (TP) rate or sensitivity on the y-axis in a ROC curve (Liuzzo et al. 2019;Pirnia et al. 2019). The methodological details of the ROC metrics have been well described by Rahmati et al. (2019) and Rahman et al. (2021). (2)

Estimation runoff coefficients
The runoff coefficient is one of the important parameters considered in surface runoff and flood discharge estimation methods in various watershed management and flood control projects (Zeinali et al. 2019). Theoretically, the runoff coefficient is a fraction of rainfall volume retained by the land, soil, and vegetation characteristics, i.e., it relates to the transformation of the total rainfall in net rainfall influenced by infiltration, canopy interception, and surface detention (Del Giudice et al. 2014). In the rational/empirical method, Cr values can be calculated based on the land use/land cover classified map in addition to the slope layer and the hydrological group of the soil units (Mousavi et al. 2019). According to Table 1, we need to verify some physical factor layers of each watershed in the study area to determine the Cr values because the runoff coefficient of different watersheds is typically determined by overlapping some land cover and slope layers (Thompson 2006;Suharyanto et al. 2021). For instance, the Cr values for a region with a cover of rangeland, slope range of below 5%, and soil unit of Entisols can be calculated as 0.25. Accordingly, the dominant situation of land cover, slope range, and soil unit were extracted for each watershed in GIS using the Zonal Statistics extension to estimate a principal Cr value for each watershed. It should be noted that the small urban and pre-urban watersheds are often represented in single land use and soil unit classes in flood risk assessments (Tierolf et al., 2021). Hence, the Cr values can be fixed in the limited values (0-1) drawn in a table (Baiamonte 2020).

Analyzing DEM and soil data
For this purpose, the elevation and slope maps, derived from DEM data, in addition to soil units map, acquired from the global soil grid dataset, were produced in GIS (Fig. 3). The dominant characteristics of each watershed were extracted using the zonal statistics in GIS in Table 2. On this basis, the surface areas of the study watersheds vary from 0.45 to 4.43 km 2 . Total 22 watersheds are located in the altitude range of 1020-1520 m above sea level. The variation of averagely length and ∆ height of the watersheds obtained equal 1.1-4.7 km and 50-400 m, respectively. About 16 watersheds (55% of the total study area) have critical slope ranges over 15%, and dominant soil units in 14 cases out of 22 watersheds (48% of the total study area) are observed as Mollisols. The soil unit of the Mollisols could be considered as the most intensive soil loss-prone zone in the study area around Mashhad city and due to its relation with the barren land and brown-fields (Ebrahimi et al. 2021).

Attributing land-cover data
According to the global land-use/land-cover (LULC) as well as HBASE dataset, derived from satellite imageries, the main types of land-cover of the study area in 2010 dominantly are categorized as rangeland area (60%) and residential real estate (40%). The percentage of land-cover types in the quo status (2020) is changed as 45% and 55% for rangeland and real estate areas, respectively. Moreover, the contribution of the rangeland and the sprawl expansion of real estate area to the land-cover types can be anticipated equal to 20% and 85%, respectively, in 2030 (Fig. 4). A dominant trend of sprawl growth in the study period (2010-2030) depends on the leapfrog expansion of the new settlements around the main roads and highways. The pattern of the new settlement in this part of the Mashhad city is contributed to suburb real estates and occupied brown-fields, which absorbed increasing urban sprawl (Yazd et al. 2019). The zonal statistics of land covers are presented for each watershed in Table 3. The urban sprawl effect and real estate cover have influenced the study area by 7 watersheds (45% of total area) and 12 watersheds (65% of total area) from 2010 to 2020, and it is predicted for 18 watersheds (> 85% of total area) up to 2030.
In Fig. 5, the satellite imageries (Google Earth) and general perspectives for a sample of the sprawl expansion are shown for two time-windows of 2010 (past) and  Figure, we can observe the creeper expansion of urban real estate around the Namaz highway (in the southern part of two isolated hills). This fact reveals that the construction of any main roads (such as the Southern highway) can trigger the sprawl expansion and impervious surfaces over the watersheds during the next time-periods.

Estimating runoff coefficients (Cr)
The Cr values in each watershed were calculated in Table 3 using the overlapping some data layers (i.e., land cover map, slope range, and soil units) and considering Table 1. In this regard, the mean values of Cr are changed from 0.50 (in 2010) to 0.65 (in 2030). The highest values of Cr (> 0.70) belong to the watersheds with real estate cover, soil unit of Mollisols, and the slope ranges over 5-15%, which will be predicted in 2030 for 7 watersheds (W1, W3, W5, W6, W13, W14, and W18) in the upstream location of the study area (Fig. 6). The highest values of Cr (> 0.70) in the present time are estimated only in a watershed (W6) in the central part of the study area. Several scholars have demonstrated that the urbanized watersheds with much more runoff coefficient (Cr) will produce much more flood discharge and susceptibility (Hung et al. 2018;Zeinali et al. 2019). Hence, predicting watersheds with high Cr values should make the urban decision-makers aware of the critical effects of sprawl expansion on the increasing flood hazards in the pre-urban areas. In the next section, we will estimate the accurate quantity of flood discharges in the study area.

Estimation peak flood discharges
Equation 1 revealed that in addition to runoff coefficients (Cr), we need to determine the precipitation rate to estimate each watershed's peak flood discharges (Qp).

Analyzing precipitation data
Precipitation rates describe the water supply to a watershed, a portion of which reaches the watershed outlet as surface runoff. The amount and duration of the precipitation rates are the essential characteristics of a storm for hydro logic analysis and can be combined to describe the intensity and frequency of the precipitation (Hayes and Young 2006). During the flash flood events in 2019, as anomalous events within the last 100 years (UNCTI 2019), some western and northeastern Iran watersheds have experienced the maximum daily-precipitation rates over 120-160 mm (Mahmoudzadeh and Daneshvar 2019). In this study, we assumed the mean daily-precipitation rates between 2001 and 2020 for the Mashhad city shown in Fig. 7. Among these data, the maximum values of precipitation rate per hour were analyzed in Fig. 8 to reveal the linear trend equation as below (same as the method used by Shanableh et al. 2018), demonstrating the increasing trend of anomalous daily-precipitation intensity (Pi) in the study area from 2001 to 2020: According to the precipitation data from 2001 to 2020 and its linear trend equation, the regression coefficient and coefficient of determination (R-squared) were obtained as 0.55 and 0.30, representing 30% of the (4) Pi = 1.24 × (Year − 2001) + 9.74 dependent variable (Pi) can be predicted by the independent variable (Year). Due to anomalous precipitation rates in recent years (2018)(2019)(2020), the goodness of fit in this equation is not maxima but is good enough to determine the future chaotic precipitation events in the study area.
The estimated maximum values of daily precipitation rate for the study area can be assumed as 20, 34, and 47 mm/h for time windows of 2010, 2020, and 2030, respectively. These estimations revealed an increasing trend for precipitation extremes ranging from 1.7 to 2.4 times in 2020 and 2030 regarding 2010. Although this range of variations in precipitation ranges depends on the peak of rainfall anomalies in the given study area and may not represent the actual climate change

Calculating peak flood discharges (Qp)
The estimation of Qp is shown in Table 4. The highest values of Qp in 2010, 2020, and 2030 were estimated equal to 4.102, 6.973, and 9.639 m 3 /s for the watershed W10. The watershed W10 is the ending outlet of all 22 watersheds, which is flowed across the zone of Ferdowsi University and is effluent in the main channel into the Malek-Abad Garden. The lowest values Qp in 2010, 2020, and 2030 were obtained equal to <0.5, <1.0, and <1.5 m 3 /s for the watersheds W8 and W14, which have the lowest surface area in the upstream location of the study area. In the entire watersheds, the average flood peak discharges increased from 1.3 to 2.4 m 3 /s during the period from 2010 to 2020 and will increase up to 3.6 m 3 /s in 2030.

Estimation outlet cross-sections Estimating time of concentration (Tc)
Based on Equation 2, the Tc values were estimated for all watersheds, as shown in Table 5. The largest Tc (1.12 h) belongs to the watershed W22 in downstream of the study area, which has the highest values of surface area (44.3 km 2 ) and river length (4.7 km). Contrarily, the shortest Tc (0.15 h) belongs to watersheds W6, W14, and W19 in the central part of the study area, which have the lowest river lengths (<1.2 km).

Determining outlet cross-sections (Sc)
Based on Equation 3, the outlet cross-sections were formulated to design flow ways and engineering control of the flood discharges downstream (Table 5). The broadest cross-section required for designing outlet channels was assumed for Watershed W10 with estimations of 17.56, 32.01, and 46.87 m 2 in 2010, 2020, and 2030. It should be noted that the peak flood flows were estimated based on the outlet of each watershed downstream, not the upstream. However, the cumulative flood flows in the outlet of a watershed were considered based on the effluents of the upstream watersheds to determine the downstream outlet cross-sections. For instance, the effluent of the entire study area needs to design with concerning a cross-section of 46.87 m 2 (e.g., with dimensions 3 × 15 m) in the future for accurate and invulnerable ejection of the cumulative flood discharges.
After the field observations, we controlled the most outlet cross-sections (Fig. 9) and found an insufficient design for cross-sections regarding the present (2020) and the future (2030) flood effluents. For example, the outlet of watershed W10 with a width of 10 m and height of 2 m is not appropriate to eject the peak flood effluents.

Flood susceptibility evaluation (FSE)
FSE of the watersheds within three time-periods of 2010, 2020, and 2030 were defined as the cumulative flood discharges of each watershed after the concentration times considering the upstream watershed discharges and outlet cross-sections ( In the entire watersheds, the average FSE increased from 2.04 to 3.84 m 3 /s from 2010 to 2020 and will increase up to 5.76 m 3 /s in 2030. In other words, the flood susceptibility will increase an average of 3.2 times in the future. Recent research indicated that the effects of the urbanized areas on flooding will increase approximately 2.7 times by the year 2030 (Güneralp et al. 2015;Park and Lee 2019). Hence, the research result in our study area revealed similar effects on the increasing flood susceptibility regarding the effects of the urban sprawl and climatic variations in the future. Globally, the previewed scenarios have demonstrated an increase in the frequency and magnitude of urban floods due to the changes in precipitation patterns and urban sprawl expansion (Zope et al. 2014;Yin et al. 2015;Park and Lee 2019).

Model validation and implication
Our results revealed that urbanization growth could result in higher impervious land areas, higher surface runoff coefficients, and higher peak flood discharges, Fig. 9 The existing structure of outlet cross-sections in the field observation for some watersheds including a W9, b W10, c W20, and d W21 adopting Feng et al. (2021). In this regard, the FSE model could present the essential estimations to determine the flood hazard quantities (statistical dimension), the floodprone zones (spatial dimension), and the needs for future design (temporal dimension). We briefly concluded that the highest flood susceptibility (> 10 m 3 /s) would increase over six watersheds (46% of the total study area) in the future. In this regard, we seriously pointed out that the prediction of the watersheds with high values of Cr should make aware the urban decision-makers about the critical effects of sprawl expansion on the increasing flood hazards in the pre-urban areas. In Fig. 11, the ROC curves for the model validity explained AUC values averagely over 0.8 for both periods of 2020 and 2030, exposing the very good performance of the model and excellent sensitivity. In detail, the ROC curve also revealed that the FSE model could be more sensitive in the prediction of future time.
The managerial implication of this research depends on prioritizing flood susceptible zones for land developers during the planning processes for urban and pre-urban regions. Urban planners, city officials, and policy-makers lacking appropriate hydrological understanding might implement equations and coefficients to create strategies against the risk of flooding in urban areas (Berndtsson et al. 2019). In this regard, the FSE model is constitutive of the current era of urbanization and must be linked to urban planning and management (Ahmad and Simonovic 2013;Su 2017). Meanwhile, the academic implication of the FSE model relates to the statistical and analytical procedures, which have the potential to supply dependent indicators concerning the runoff coefficient (Cr) and time of concentration (Tc). The localized GIS-based works are proposed to determine the Cr and Tc values in the urban hydrological integrations. Besides, comparative research is recommended to detect the triggering role of climatic parameters and urban sprawl effect on flood susceptible zones.

Conclusion
The severe flooding of 2019 in Iran was a great momentum for the urban decision-makers and dwellers to consider the forthcoming urban flood hazards in the pre-urban and urban watersheds. In the present paper, we evaluated a flood susceptibility model in the southern watersheds of Mashhad city, in Iran, during three temporal windows of 2010, 2020, and 2030. As mentioned by NRIDR (2005), to integrate the GIS and RS for detection of flood susceptible prone zones in the national flood risk programs of Iran, we constructed an FSE model based on the utilization of RS and GIS tools to evaluate spatial and temporal levels of runoff variations, peak flood discharges, and flood hazard susceptibilities.
The spatial and statistical analysis revealed that about 16 watersheds (55% of the total study area) have critical slope ranges over 15%, and dominant soil units in 14 cases out of 22 watersheds (48% of the total study area) are observed as Mollisols, the most intensive soil units regarding erosion and weathering. Furthermore, the urban sprawl and real estate cover have influenced the study area by 7 watersheds (45% of total area) and 12 watersheds (65% of total area) from 2010 to 2020, and it is predicted for 18 watersheds (> 85% of total area) up to 2030. Our results revealed the mean values of runoff coefficient (Cr) from 0.50 (2010) to 0.65 (2030), where the highest values of Cr (> 0.70) belonged to the watersheds with real estate cover, soil unit of Mollisols, and the slope ranges over 5-15%. According to the precipitation data from 2001 to 2020 and its linear trend equation, the estimated maximum values of daily precipitation rate were assumed as 20, 34, and 47 mm/h for time windows of 2010, 2020, and 2030, respectively. Consequently, the averagely cumulative flood discharges increased from 2.04 to 3.84 m 3 /s during the period from 2010 to 2020 and will increase up to 5.76 m 3 /s in 2030. The flood susceptibility evaluation (FSE) will increase an average of 3.2 times in the future. Although the climatic variations in Iran exposed the increasing extreme events such as urban floods (Daneshvar et al. 2019a), the structure and design of urban areas can influence the intensity of local climate disasters (Ren 2015). In the study area, i.e., Mashhad city, Daneshvar et al. (2019b) have already noted the relationships between urban sprawl and local climatic factors. Similarly, our findings revealed both effects of climatic anomalies and sprawl expansions on the intensity of flood disaster and susceptibility in the study area.