Comprehensive analysis of the impacts of climate and land use changes on flood damage (Case study: AqQala, Iran)

Industrial development and population growth have accelerated deforestation and environmental degradation with a significant role in rising greenhouse gas emissions over the past decades. Climate change impacts have been intensified in low latitude regions, i.e., Iran, on variables such as temperature and precipitation (Hajarizadeh et al. 2013; Thompson et al. 2021). Broadly speaking, the characteristics and attributes of floods are highly dependent on the location and climatic conditions of the study area, and as a result, the approaches and techniques used by different researchers differ accordingly, upon their research objectives. Climate change can alter the spatiotemporal distribution of precipitation and thus affect the quantitative and qualitative characteristics of water resources (Behboudian et al. 2023; Uniyal et al. 2023). In recent years, land use change has also greatly affected the scale of floods and damage caused by them (Liu et al. 2023; Zhao et al. 2023). Many researchers have emphasized that land use change has a fundamental impact on the scale of flood damage as important as climate change (Rogger et al. 2017; Ngondo et al. 2021; Sugianto et al. 2022; Sun et al. 2023; Xiao et al. 2023). Some others concluded that land use change has a more powerful effect on flood intensity than climate change (Ward et al. 2008; Agarwal et al. 2022). Furthermore, Iran is facing deforestation which exacerbates the impacts of climate change by reducing carbon sequestration capacity and increasing the vulnerability of landscapes to environmental hazards such as soil erosion, floods, and droughts (Gholoubi et al. 2019).

The influences of climate and land use change have been widely investigated from water resources aspects including surface and underground water. In terms of surface water, these impacts are studied independently or through synergistic approaches. The results of most studies indicate an increase in flood characteristics. To exemplify, the impact of climate change exclusively is studied by Bai et al. (2019), and the higher emission factors result in a more dramatic increase in all future precipitation and flood risks compared with historical observations. For land use change impact, the hydrologic model developed by Hussein et al. (2020) indicated that increasing flooding potential depends on the rate of change in urbanization, along with the magnitude of rainfall events and watershed characteristics. That said, there are some limited investigations that report decreasing or insignificant changes in flood characteristics (Mfwango et al. 2022). Insignificant impacts on the extent of floods are reported by Kuntiyawichai et al. (2020), when both influences of climate and land use change are investigated together. The result differences under representative concentration pathways (RCP) 4.5 and RCP 8.5 were concluded to be negative and positive for inundation compared with baseline, respectively. For climate model setups, a handful of investigations have used CMIP6 for anticipation processes while most research is conducted in this area via CMIP5. The CMIP6 models provide an improved dynamical process with finer resolution, and future climate change simulations are created using the new SSP-based emission scenarios (Eyring et al. 2016; O'Neill et al. 2016). In such studies, a wide variety of hydrologic and hydraulic models are used to achieve different research objectives. Several models have been used for rainfall-runoff modeling (R-R). The HEC-RAS hydraulic model has been mainly used to prepare flood inundation maps. There are studies available with the main focus on the estimation of flood damage, most of which have concluded an increase in future (Hooshyaripor et al. 2020; Kuntiyawichai et al. 2020; Pariartha et al. 2023). For instance, Kuntiyawichai et al. (2020) demonstrated that considering both the maximum flood depth and flood duration in the calculation of total direct flood damage resulted in a clear increase in direct flood damage for higher return periods (RPs) under both baseline and future conditions. The study highlighted the importance of incorporating these factors to better assess the potential impacts of flooding in the region. The results indicated that, under future conditions, the total direct damage is projected to be higher compared to the baseline scenario, with the magnitude of increase varying depending on the RP and the climate scenario (RCP 4.5 or RCP 8.5). These findings emphasize the significance of accounting for various flood characteristics to improve flood risk assessments and planning for future scenarios.

On 26 March 2019, due to heavy rains, floods occurred in the northern provinces of Iran. These floods caused dramatic financial losses and casualties by these floods with hundreds of people injured. In this natural disaster, 300 mm of rainfall in two days lashed part of these regions, which is equal to the average annual precipitation for Golestan province itself. Reports indicate that adjacent regions have also suffered rainfall of 50–70% of the total annual precipitation in less than 5 days. This record has been the highest reported in a period of nearly 70 years of meteorological data available for this region (Karami 2019).

This disaster seems to be the result of simultaneous influences of climate and land use change, in which their prediction-uncertainties make the situation more complicated to study. Such conditions are considered as ‘Deep Uncertainty’. This type of uncertainty will lead to inappropriate and even more conflicting decisions. Therefore, in order to define an appropriate and effective decision-making process for such challenges caused by uncertainties at this level, it is vital to develop methods that help planners and policy makers in providing long-term master plans. The main findings of the present study are to facilitate a systematic approach to properly formulate adaptation strategies and manage the vulnerabilities imposed by future flood risks and damage within deep uncertainty conditions.

The present investigation has been conducted in five main steps which are explained in detail in each related section. Briefly, these steps are as follows: (1) Anticipation of temperature extremes derived from the climate model; (2) Projection of land use changes in 2050 and 2080; (3) Hydrological modeling of the study area to obtain the hydrograph in each of the defined scenarios; (4) Hydraulic modeling to extract flood inundation maps; and (5) Direct damage estimation in the worst-case scenario while comparing with the status quo.

Daily precipitation is collected from nine synoptic stations including Gorgan, Marave, Gonbad, Kalale, Aliabad, Hashem, Torkaman, Inche, and Minudasht for the past years (1975–2020). The precipitation data and temperature extremes are obtained from the Iranian Meteorological Organization. The anticipation of future precipitation and temperature changes (minimum and maximum) is applied for the anticipation of future discharge computed by the SSA model.

In this research, the CMCC-ESM2 model's output has been employed to investigate three SSPs provided by IPCC (CMIP6). Different GCM models including HadGEM3, AccESSESM1, MRI-ESM2, and CMCC-ESM2 were investigated in this region. Our findings suggest that CMCC-ESM2 yielded better results for climate change studies. Specifically, the study focuses on SSP126, SSP370, and SSP585 scenarios. The CMCC-ESM2 model is utilized because of its remarkable performance with positive implications for future climate change scenarios (Cherchi et al. 2019). SSP126, SSP370, and SSP585 scenarios are selected to candidate optimistic option, average mode, and the most critical scenarios, respectively.

In the given context, H_i represents the observed values, while F_i denotes the forecasted values. H_avg and F_avg correspond to the average values of the observed and forecasted data, respectively. The parameter N represents the total number of observations.

Finally, precipitation, maximum temperature, and minimum temperature data have been produced corresponding to three climate change scenarios.

For the preparation of land use maps, a combination of three software packages is utilized, namely ENVI 5.3 software, Landsat satellite imagery acquired from the Enhanced Thematic Mapper Plus (ETM+) sensor for the 1991 dataset, and Landsat 8 imagery obtained from the Operational Land Imager (OLI) sensor for the 2002 and 2018 datasets. Prior to processing, geometric and radiometric corrections are considered on images using the ENVI software similar to Birhanu et al. (2019). It is critical to prepare high-resolution satellite images to obtain an unambiguous description of land uses in the study area. Next, validation of the present status for existing land uses has been carried out categorized by forest, agriculture, urban, bare land, pasture, and water. Then, samples are chosen from satellite images for each category, and they are then split into training (70%) and test (30%) sets to facilitate validation. The acquisition of these samples from satellite images is performed in a random manner, taking into account visual interpretation with the assistance of topographic maps and images from Google Earth (similar method to Tarawally et al. (2019)). The training samples make up a portion of the total image pixels, with their contribution capped at 10% (as recommended by Hanifehlou et al. (2022)). The supervised classification method of the maximum likelihood algorithm is employed using the ENVI software to generate land use maps for the years 1991, 2002, and 2018. This classification approach relies on the variance and covariance of classes to evaluate the process (as explained in Islam et al. (2018)). This method provides a robust and reliable approach for generating accurate land use maps based on remotely sensed imagery, making it a suitable choice for our study, and it has been used in some previous studies in this region (Hanifehlou et al. 2022). To achieve accurate land use classification in an area, it's essential to take into consideration the characteristics of the region, existing land cover types, their distribution, and spatial patterns. Additionally, some covers may have similar spectral reflection, which can affect soil reflection. Due to these factors, an automatic analysis method may not provide reliable results. As a result, the maximum likelihood method was employed as a supervised classification method to achieve high-accuracy classification of the mentioned area. Afterward, the results extracted from prepared land use maps are evaluated and validated based on the Kappa coefficient applied for 30% of the gathered samples. The Kappa coefficient for 1991 is 74.33. For 2002, it is 4.92% higher than 1991, and the coefficient for 2018 is even a bit more pronounced. The coefficient's higher value in 2018 is attributed to the improved radiometric resolution of the images compared with the 1991 and 2002 cases.

The Land Change Modeler (LCM) program of TerrSet software imports the land use maps for 1991 and 2002, as well as the digital elevation model (DEM) and slope maps of the study area to forecast future land use changes. Logistic and artificial neural networks (ANNs) accessible in the LCM are applied to model transition potential maps within the TerrSet software. The land use map for 2018 is anticipated by the Markov chain method after preparing sub-models. This method creates a stochastic model characterized by a sequence of potential outcomes, where the probability of each outcome is contingent on the outcomes in proximity and described through a transition matrix (Jokar Arsanjani et al. 2013; Noszczyk 2018). To evaluate the effectiveness of the developed model, the anticipated map of 2018 is compared with the current map of this year and validated using the Kappa coefficient. A satisfactory performance of the model is indicated by the Kappa coefficient of 0.92. Subsequently, the validated model is applied to forecast the land use maps for 2050 and 2080. The spatial data preparation is carried out by applying ArcGIS 10.8 software.

The estimation of streamflow in the Gorganrud basin is conducted by applying the SSA model, which is based on HEC-1 principles and specifically designed for rainfall-runoff (R-R) processes. The classification of soil hydrological groups is extracted from the map provided by NASA (HYSOGs250m). In ArcGIS 10.8 software, the Gorganrud basin is delineated and subdivided into multiple sub-basins using Arc Hydro Tools, incorporating a 30-meter DEM and stream network. The study area is partitioned into 25 sub-basins for analysis. The SSA model's hydrologic elements, such as sub-basin, source, junction, reservoir, sink, reach, and diversion, are interconnected in the river network, providing insights into the runoff processes and their impact on the drainage system (Scharffenberg 2016). Additionally, three reservoir elements are set to simulate hydrograph detention and attenuation.

In this study, the Soil Conservation Service (SCS) Curve Number (CN) loss method is utilized. The SCS-CN method applies the curve number approach for incremental losses. Originally designed to calculate total infiltration during a storm, the methodology is adapted here to compute incremental precipitation by reevaluating the infiltration volume at the end of each time interval. Infiltration for each interval is derived from the difference in volume between the end of two adjacent time intervals.

In the context of the study: T_c refers to the time of concentration in hours, L represents the hydraulic length of the watershed, which is the longest flow path in feet, CN denotes the SCS runoff curve number, and S stands for the average land slope of the watershed in percentage.

Statistical tests are performed based on the long-term data of maximum annual 24-h rainfall which are extracted from the climate model. The frequency analysis is conducted via HyfranPlus software. Following the selection of the best statistical distribution, the 24- and 36-h precipitation height values in the RPs of 2 to 100 years are considered in different scenarios. Universal methods such as p-value, χ2, … are used for comparing statistical tests. In this research, we used the ‘Ghahraman method’ (Ghahraman & Abkhezr 2004) to convert daily rainfall into hourly values. Considering the fact that the T_c of the basin is more than 24 h, and the duration of rainfall should be more than the T_c, 36-h rainfall is applied for R-R modeling.

Having three SSP scenarios (SSP126, SSP370, and SSP585) and three land use maps for 2018, 2050, and 2080, there is a need to define combined scenarios to include the simultaneous impact of climate and land use change. To accomplish this task, the characteristics of the scenarios for SSA are defined in Table 1 based on climate scenarios, different land use maps, and the target years for anticipation (2050 and 2080) in this research. Scenario 1 serves as the baseline, representing the present situation without any modifications.

Hydraulic modeling is used to obtain the depth and velocity of flood discharge in each scenario to prepare a corresponding inundation map in this research. The HEC-RAS hydraulic model is applied to calculate water surface profiles for steady, gradually varied flow in natural or constructed channels in the RPs of 25, 50, and 100 years. The water surface profiles are calculated using the standard step method where each cross-section is built based on the previous one by iteratively solving the energy equation (similar to Brunner (2016)). The geometric data editor in HEC-RAS is utilized to build a schematic of the connected river system. Then additional hydraulic elements including junctions and structures are depicted. In the HEC-RAS models, all cross-section attributes are determined based on the terrain profile description of the river, considering the potential for flooding on both banks. These attributes include overbanks, stations, elevations, Manning's roughness coefficients, reach length, left and right banks, as well as energy loss coefficients, encompassing friction, expansion, and contraction losses. To initiate water surface calculations in HEC-RAS, boundary conditions are defined for each reach, both upstream and downstream. Moreover, internal boundary conditions are considered for connections to junctions. The flow data from each sub-basin, computed by the SSA model, along with the released flow data from dams, extracted from the HEC Data Storage System (DSS), are brought into the HEC-RAS models.

The so-called Normal Depth boundary condition is applied for the downstream section in the HEC-RAS setup model. It assumes normal flow (uniform flow) conditions at the downstream boundary of the river. This assumption allows for the calculation of the normal depth or stage for each calculated flow using Manning's equation, with the main channel's bottom slope set at 0.002%. To evaluate the models, a validation process is conducted to investigate whether they can accurately represent physical processes and make precise predictions under various conditions. Maximum daily discharge and annual peak flow rates during the years 2004–2020 are used as observed data for validating the SSA peak runoff rates as simulated data for the AqQala basin. R² and NSE are used as reliable criteria to assess the model performance.

Here are the variables and their meanings: DFD stands for direct flood damage in US dollars (US$). DDL (j, L, D) represents direct flood damage per land use type j at L, D, measured in US$ per land use type. AAL (j) denotes the average area per land use type j per unit in square meters (m²). PL (i, j) represents the percentage of land use type j in cell i, dimensionless (–). A (i) refers to the area of cell i in square meters (m²). i is the cell number. j represents the land use types 1, 2, and 3 (agriculture, bare land, and urban). L represents the maximum flood depth in meters (m). D is the flood duration in days.

The calculation provided above does not include direct flood damage to infrastructure; the total direct flood damage (TDM) can be estimated considering that infrastructure damage cost can be as high as 65% of DFD (as suggested by Reinsurance World Map of Natural Hazards (Berz et al. 2001)). In summary, the decision to exclude direct flood damage to infrastructure in Equation (4) was influenced by a combination of data availability, complexity, scope, and methodological considerations.

It is worth clarifying that the calculation of agricultural land damage is based on the global price of the dominant crop of the study area (wheat). The average efficiency of the farmlands in a crop year is considered to be 4.2 tons per hectare, worth roughly 522 USD per ton. Also, the calculation of damage to the urban area is set to be 8.24 USD per square meter based on the flood damage calculation issued by the National Flood Report for 2019 flooding of Golestan province.

The average annual precipitation of historical data for the study area is 495 mm. The averages derived from the analysis of the climate model data for the same area are 614, 723, and 634 mm under the climate scenarios SSP126, SSP370, and SSP585, respectively, until 2080. The maximum and minimum annual temperatures of historical data for the AqQala area are 23.13 and 12.67 °C, respectively. The maximum annual temperature is anticipated to be 25.34, 25.72, and 26.64 °C under the climate scenarios SSP126, SSP370, and SSP585, respectively, until 2080. Similarly, the minimum annual temperature is anticipated to be 14.65, 15.14, and 15.94 °C under the climate scenarios SSP126, SSP370, and SSP585, respectively, until 2080.

Accordingly, the climate model data analysis for the AqQala area shows an increase in average annual precipitation by 24, 46, and 28% under the climate scenarios SSP126, SSP370, and SSP585, respectively, until 2080. Also, the maximum annual temperature has increased by 2.21, 2.59, and 3.50 °C under the climate scenarios SSP126, SSP370, and SSP585, respectively, until 2080. Similarly, the rise of minimum annual temperature is anticipated to be 1.98, 2.46, and 3.27 °C under the climate scenarios SSP126, SSP370, and SSP585, respectively, until 2080.

The future climate change was anticipated for the period of 2020–2080 in comparison to the baseline period of 1975–2020, under a set of shared socioeconomic pathways (SSPs) greenhouse gas scenarios (SSP126, SSP370, and SSP585) from CMIP6 model simulations (CMCC-ESM2 for this study). Following the high-resolution downscaling (1° grid spacing), the input is used for the calibrated SSA model to generate a future hydrograph, in which the flood inundation maps were extracted by the calibrated HEC-RAS model.

The areas of forest, agriculture, urban, bare land, and pasture lands in the study area in 2018 are 2,170, 2,972, 500, 2,429, and 1,989 km². For 2050, the area of the land uses with the same abovementioned order will be 1,961, 2,177, 1,052, 3,029, and 1,847 km². And, for 2080, the areas of these land uses are 1,563, 948, 2,582, 3,381, and 1,594 km², respectively.

Briefly speaking, the anticipated results for land use changes of the AqQala area in 2080 compared with 2018 indicate an increase in urban and bare lands by 416 and 39%, respectively, where, at the same timeline, a decrease of 28, 68, and 20% is anticipated for forest, agricultural, and pasture lands, respectively. Generally speaking, the land use for forest, agriculture, and pasture lands will be decreased while the land use for urban and bare land will be increased in the future.

From a practical standpoint, accurately projecting land use conditions for 2050 and 2080 may present significant challenges. However, it is essential to recognize that aforementioned models possess the capability to forecast future conditions based on available statuesque inputs. It is also evident that the forecast for land use cover leans toward a pessimistic outlook. An important observation pertains to the increase in bare land despite projected future increments in precipitation. It is noteworthy that the used land use models operate autonomously of anticipated precipitation augmentations due to their inherent limitations in predictive capacity; they are reliant on learning from the past for future projections. On the one hand, rainfall is crucial for vegetation growth, but on the other hand, alterations in rainfall patterns, intensity, and distribution can have intricate and sometimes harmful effects on soil erosion and land degradation. In addition, temperature increase during the projection period should be considered along with participation as it may contribute to the expansion of barren land areas.

Evidently, provided that the trend of land use changes in the study area continues, with the increase of urban and bare lands and the decrease of forest, agriculture, and pasture lands, the permeability of the land will decrease, and it will become surface runoff. Subsequently, as confirmed by the results, much larger floods will cause more damage.

The examination of rainfall variations is performed at various RPs. The baseline period encompasses 46 years of data from historical simulation runs (1975–2020), while the future climate is projected using a 60-year period (2021–2080) based on the described scenarios. The 25-, 50-, and 100-year RPs of annual 24-h maximum rainfall are estimated while applying the different distributions under three SSPs of SSP126, SSP370, and SSP585 in the AqQala basin. Annual 24-h maximum rainfall characteristics and precipitation values in different RPs for each SSP scenario are provided in Table 2.

By observing the results related to 15 scenarios (output of the R-R model), the flood peak discharge values for scenarios 1–15 for 100-year RP are 637, 499, 579, 531, 780, 615, 711, 654, 454, 502, 464, 876, 636, 700, and 649 m³/s, respectively. The top three, in terms of flood intensity, are S12, S05, and S07, from high to low. Among 15 scenarios, the rank of the S01-PY2018-PrBl-CN2018 (present situation) scenario is seven.

By observing the results related to 15 scenarios (output of the R-R model), despite the increase of precipitation and temperature amounts in three climate scenarios SSP370, SSP585, and SSP126, assuming that other variables are constant, the amount of flood peak decreases. Even though the increase in precipitation amounts has been observed in the three climate scenarios, the fluctuation of changes in the average annual forecasted precipitation values by climate models in the coming years is much less than the fluctuation of precipitation values in the present situation.

However, with the change of land use maps in 2018, 2050, and 2080, the flood peak discharge values will increase chronologically. Therefore, it can be stated with certainty that the effect of land use change on floods in the future is greater than climate change.

Hydraulic modeling has been specifically concentrated on a 14.5-km segment of the Gorganrud river as it flows through AqQala city. After validation and calibration of the flood discharges calculated in the previous section, the calculated values of two scenarios S01-PY2018-PrBl-CN2018 and S12-PY2080-PrBl-CN2080 are selected for hydraulic modeling using HEC-RAS software. The S01-PY2018-PrBl-CN2018 scenario is chosen due to the representation of the present situation and the closeness of the flood values obtained to the average values of the 15 scenarios. The S12-PY2080-PrBl-CN2080 scenario is considered to be the worst-case scenario because of the maximum flood values. The map of flood depth and velocity changes in the present situation and maximum flood scenarios in the 25-, 50-, and 100-year RPs are prepared for estimating flood damages.

The total flooded areas for the 25-, 50-, and 100-year RPs in the S12 scenario are 682, 755, and 769 ha, respectively. The total flooded areas for the 25-, 50-, and 100-year RPs in the S01 scenario are 422, 546, and 816 ha, respectively. The total coverage of flooded areas for the 25- and 50-year RPs in the S12 scenario, compared to the S01 scenario, has increased by 62 and 38%, respectively. The total coverage of flooded areas for the 100-year RP in the S12 scenario, compared to the S01 scenario, has decreased by 6%.

The simulation findings indicated that the potential flooded area in AqQala with 25-, 50-, and 100-year RPs are 422, 546, and 816 ha, in the present situation, respectively. Having future scenarios in scope, the extent of flood inundation of the worst-case scenario (S12) is larger than the baseline case for 25-, 50-, and 100-year RPs by 682, 755, and 769 ha, respectively.

In terms of comparing the existing situation (S01: average value) with the worst-case scenario (S12: maximum value), the total damages of flooded areas for the 25-, 50-, and 100-year RPs in the S12 scenario are 53.2, 64.5, and 66.3 million dollars, respectively. The total damages of flooded areas for the 25-, 50-, and 100-year RPs in the S01 scenario are 20, 39.3, and 51.4 million dollars, respectively. The scale of damage over flooded areas for the 25-, 50-, and 100-year RPs in the S12 scenario, compared to the S01 scenario, has increased by 166, 64, and 29%, respectively.

This study examines the simultaneous impact of climate and land use change on floods in the AqQala area. Considering such concurrent influences, models project future floods until 2080, compared to 1975–2020. Despite increased precipitation and temperature in the future, flood peaks decrease. Land use changes, particularly urban expansion, and forest decrease, contribute significantly to increased surface runoff and larger floods. Land use changes emerge as a more influential factor on future floods and damages than climate change. Simulation shows potential flooded areas and damages under different RPs, with the worst-case scenario projecting significant increases in inundation and damages compared to the baseline.

When developing a model to predict land use change, it is essential to take into account the strategic documents that have been created for the specific region, if any exist. These documents provide valuable information on the region's development plans and policies, which can help in identifying the patterns and factors that may influence land use change. Therefore, it is highly recommended to consider any available strategic documents in the region when developing a model for land use change. It may be challenging to accurately predict long-term land use conditions, but we acknowledge that models can project future conditions based on present inputs. This is an important factor to consider for practical purposes in this area.

To elaborate on the challenges regarding deep uncertainty and mitigate the conflicting decisions, this study can be a good starting phase to enhance the improvement plans on how to increase the water carrying capacity for this case in the future studies. Harmoniously, to reach more feasible low-cost eco-friendly plans, it is necessary to implement suitable smart corrective plans in order to choose the best possible decision. That said, in the continuation of this research, to develop flexible adaptation pathways, as a creative decision-making tool for policy makers, can be suggested to improve the existing situation.

The authors are grateful for the contributions by Dr Yashar Madjidi, and Dr Saeed Movahedi.

Conceptualization: A.M. and S.A.H.; methodology: A.M. and S.A.H.; software: A.M., S.A.H., and B.G.; validation: A.M., S.A.H., and B.G.; formal analysis: A.M. and S.A.H.; investigation: A.M. and S.A.H.; writing – original draft preparation: A.M.; writing – review and editing: A.M., S.A.H., and S.N.; supervision: S.A.H. and S.N.

No funding for this research.

The manuscript has been read and approved by all named authors and there are no other persons who satisfied the criteria for authorship but are not listed. The order of authors listed in the manuscript has been approved by all the authors.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare no competing interests.