Flood risk mapping and crop-water loss modeling using water footprint analysis in agricultural watershed, northern Iran

Spatial information on flood risk and flood-related crop losses is important in flood mitigation and risk management in agricultural watersheds. In this study, loss of water bound in agricultural products following damage by flooding was calculated using water footprint and agricultural statistics, using the Talar watershed, northern Iran, as a case. The main conditioning factors on flood risk (flow accumulation, slope, land use, rainfall intensity, geology, and elevation) were rated and combined in GIS, and a flood risk map classified into five risk classes (very low to very high) was created. Using average crop yield per hectare, the amount of rice and wheat products under flood risk was calculated for the watershed. Finally, the spatial relationships between agricultural land uses (rice and wheat) and flood risk areas were evaluated using geographically weighted regression (GWR) in terms of local R2 at sub-watershed scale. The results showed that elevation was the most critical factor for flood risk. GWR results indicated that local R2 between rice farms and flood risk decreased gradually from north to south in the watershed, while no pattern was detected for wheat farms. Potential production of rice and wheat in very high flood risk zones was estimated to be 7972 and 18,860 tons, on an area of 822 ha and 7218 ha, respectively. Loss of these crops to flooding meant that approximately 34.04 and 12.10 million m3 water used for production of wheat and rice, respectively, were lost. These findings can help managers, policymakers, and watershed stakeholders achieve better crop management and flood damage reduction.


Introduction
Flooding is one of the most devastating and costly natural hazards. It has severe socioeconomic and environmental consequences, including destroying farmland, reducing crop yield, and causing regional freshwater shortages (Mind'je et al. 2019). Flooding can only occur not only in lowland areas but also in mountainous environments. Analysis of flooding and its relationships with explanatory variables can help water managers identify the most effective variable in flooding (Hosseini et al. 2020). Therefore, identification of floodprone areas and of the most influential conditioning factors is an essential tool in mitigating the effects of flooding (Khosravi et al. 2016). Various models have been developed to estimate potential flood risk areas and simulate flow Rahmati et al. 2019). However, some of these models cannot simulate and evaluate flood risk under different scenarios, while in some heterogeneity of input data, e.g., on land use and geology, makes determining the threshold flow more difficult (Zhao et al. 2018). Other models have limitations especially in ungauged and extensive areas. As a solution, multi-criteria analysis (MCA) is a useful first step in mapping and evaluating flood risk areas. In MCA, it is possible to select, evaluate, and combine relevant factors to map the final flood risk (Santos et al. 2019). Understanding the particular interrelationships between water and food can enhance the resilience of water-food systems, since food security is highly associated with water and both are affected by a changing climate. Therefore, quantifying the effects of floods on crop production, and consequently on food security, are important (Pacetti et al. 2017). These detrimental effects can influence food availability, identified by the Food and Agriculture Organization of the United Nations (FAO) as a pillar of food security (FAO 2015). Therefore, it is important to manage irrigated and rainfed agricultural systems in light of the relationship between agricultural land uses and flood-prone areas, in order to maintain the ecosystem service of food production. In an MCA approach to identify the effect of flooding on irrigated and rainfed agricultural crops in the present study, the "water footprint" (WF) concept introduced by Hoekstra (2003) was used. The WF of an agricultural crop is defined as the total volume of water (rainfall or irrigation) used to produce the product (Hoekstra 2009). The damage caused by floods to crop production can thus be converted to WF as a complementary indicator. Some previous case studies have evaluated the effects of flood events on food availability, e.g., Pacetti et al. (2017) calculated the damage of flooding to agricultural areas in Bangladesh and Pakistan in terms of lost calories. However, previous studies have estimated flood damage using conventional statistical analyses, which produce average parameter estimates, and thus spatial variations in flood damage to crops and their associated WF have been ignored. A more sophisticated analytical technique was needed to overcome this limitation. In response, geographical weighted regression (GWR), a statistical model, was developed to investigate spatial correlation and heterogeneity (Xia et al. 2018). GWR examines spatially non-stationary parameters, so model performance is improved by reducing spatial autocorrelations. In recent years, GWR has been widely employed to good effect in different fields, to analyze, e.g., groundwater quantity (Taghipour Javi et al. 2014;Almeida et al. 2018), rainfall and environmental indices (Georganos et al. 2017;Ahmadi et al. 2018b;Salimi et al. 2018), land surface temperature (Kalota 2017;Zhao et al. 2018), urban and regional differences (Dadashpoor et al. 2019;Duncan et al. 2019) and ecology and human geography (Tu 2011;Li et al. 2018).
The staple foods in Iran are rice-based and wheat-based foods (Nejad et al. 2011;Karizaki 2016) and the Talar watershed in Mazandaran province, northern Iran, is one of the most important regions for domestic cultivation of these crops. Therefore, these two major crops and the Talar watershed were selected as a case for calculating the risk of potential crop-water losses due to flooding using the WF approach. Floods cause severe damage to agricultural land and residential areas in Mazandaran province , with 70% of available credit going to repair the damage caused by flooding (Sadeghi-Pouya et al. 2017). In this MCA study, the crop damage and associated water loss caused by flooding were calculated by integrating satellite images, field data, and agricultural for northern Iran. Specific objectives of the study were to: (i) create a flood risk map using conditioning factors; (ii) determine the water losses associated with flood-related losses of wheat and rice, through the WF approach; and (iii) explore the spatial relationships between area of flood risk classes and areas of main crops (e.g., wheat and rice) in the study watershed using the GWR technique.

Study area
The Talar watershed in northern Iran (36° 36′-36° 46′ N; 55° 23′-54° 31′ E) extends along the coast of the Caspian Sea, to which the Talar river drains in a south-north direction (Fig. 1). This mountainous watershed (mean altitude ~ 1800 masl) covers 2055 km 2 and has a Mediterranean rainfall regime, with mean annual precipitation of 552.7 mm, the majority falling in spring. Rice and wheat, which together account for about 23% of agricultural land, are the main crops (Maghsood et al. 2019). The terrain is mountainous, with steep slopes covering more than 60% of the area. The main land uses are rangeland, forest, Fig. 1 Maps showing the location of the Talar watershed in northern Iran rainfed agriculture, irrigated agriculture, and residential area (Shooshtari and Gholamalifard 2015). The watersheds in northern Iran have experienced drastic changes due to mining activities, road construction, and residential development. Large areas of forest and rangelands have been converted to cropland, orchards, and residential areas in the past four decades (Kavian et al. 2018;Pirnia et al. 2019). These changes have affected the hydrological response, ultimately altering runoff volume and flow regime in watersheds ). In the Talar watershed, recent decades due to intensive deforestation and land use change, the potential of runoff generation and peak flow have been increased 12.38% and 41.8%, respectively (Khaleghi 2017). In Table 1, characteristics of rainfall station are presented.

Data sets
For the purposes of this study, Talar watershed border and sub-watersheds of the Talar River were delineated using digital elevation model (DEM) with resolution 30 m and Arc-SWAT extension 2012.10_4.21v. A land use map of the Talar watershed was extracted using Landsat 8/Operational Land Imager (OLI) image (21.06.2016) and supervised classification algorithm in the ENVI 5.3 (Torabi Haghighi et al. 2018). Five land use types were identified: forest, rangeland, irrigated agriculture, rainfed agriculture, and residential area.
The spatial variability in flood risk was determined based on six factors ( Fig. 2): slope (S), elevation (E), flow accumulation (F), geology (G), land use (L), and rainfall intensity (R), see Eq. (1), which have direct important impacts on flood risk Karatzas 2011, 2016;Kazakis et al. 2015). These factors were prepared as six thematic maps, which were then combined into one final flood risk map using linear algebraic function and their weights in the GIS environment. Factors F and G are qualitative, whereas factors S, E, and R are considered quantitative. The factors were classified into five flood risk zones (FRZ): very low, low, moderate, high, and very high. Jenk's natural breaks method was applied to classify quantitative (numerical) factors, whereas qualitative factors were classified based on their effect in flood recharging. For example for factor L (land use), forest was classified as very low flood risk, but residential area as very high flood risk. For each classified factor's flood risk rating, a numerical value was allocated (very low (1), low (2), moderate (5), high (8), and very high (10)) (Kourgialas and Karatzas 2016).
Since all factors do not have the same effect on flooding condition, two types of effects were considered: minor effect, where a change in one factor has an indirect effect on another factor (allocated 0.5 points), and major effect, where a change in one factor has a direct effect 1 3 on another factor (allocated 1 point). The final rating of each factor was calculated by adding together the points allocated to minor and major effects ( Table 2). All factors, FRZs, and points scores used were based on previous studies (Kourgialas and Karatzas 2011;Kazakis et al. 2015). Evaluation of the literature indicated that the six selected factors provide necessary and useful information for improving flood risk modeling (Zerger 2002). The final score for each factor (Table 2) was obtained by multiplying proposed weight of effect (RL) by factor rate (FR). Finally, a flood risk map of the Talar watershed was prepared using linear summation method of factors (FREGLS) in GIS as (Kourgialas and Karatzas 2011): where x i is the thematic map of each factor i, w i is the weight of each factor i and F, S, L, R, G, E are flow accumulation factor, slope factor, land use factor, rainfall intensity factor, geology factor, and elevation factor, respectively. Using DEM with spatial resolution 30 m, the elevation (E) was classified, and a slope (S) map was produced using the 3D Analyst tool in GIS. Land use and land cover affect the volume of runoff, with a direct impact on both the time for which soil receives rainfall and the amount (Kazakis et al. 2015).
The geology (G) and land formation-based classification of the study area were produced based on the Iranian geological map (Geological Survey of Iran, 1997). Permeable geological formations with high porosity and fractures reduce runoff coefficient and runoff volume, whereas fine-grained and impermeable formations increase runoff volume and flood risk (Table 2).
Flow (F) accumulation was generated using DEM and the spatial analyst tool in the GIS environment. In this map, pixels with higher values have more hydrological connections with other pixels, and thus make a higher contribution on flood risk.
In order to calculate rainfall intensity factor (R), precipitation data from seven rainfall stations in the Talar watershed were used. The map of rainfall intensity was created using Modified Fournier Index (MFI) (Morgan 2005), calculated as: where p is the mean monthly rainfall for month i (1 ≤ I ≤ 12) and P is the mean annual rainfall. MFI shows the average monthly rainfall for stations, calculated using interpolated spline method. This method is a useful way to indicate spatial variation in, e.g., rainfall, especially in data-limited situations (Lloyd 2005). MFI can be used for the Mediterranean rainfall regime, in which flashy floods and overflow of stream banks are common (Belmonte and Beltrán 2001; Kourgialas and Karatzas 2016). All thematic maps were prepared with spatial resolution of 30 m × 30 m and summed with the Raster calculator tool in the GIS environment considering Eq. 1.

Flood damage to agricultural products and associated water loss
The effect of flood damage on agricultural products in terms of crop and water (irrigation water) loss was estimated using spatial data such as land use maps and agricultural statistics (Brémond and Grelot 2013;Giang et al. 2020). Determining the damage to agriculture due to flooding involved estimating loss of agricultural products, taking into account the topography, land use, and characteristics of cultivated crops in the study area (Pacetti et al. 2017). The effect of flooding on crop production includes loss of the crop itself, reduced food security, loss of energy, and loss of water used for production of the crop. Thus, the percentage of different land uses in each sub-watershed in the Talar watershed was calculated. According to the annual agricultural organization report (Ahmadi et al. 2018a), irrigated agriculture in this region is mostly rice, grown near the main river, while rainfed agricultural lands is mainly devoted to wheat production. Thus, rice and wheat land were assumed for irrigated and rainfed agriculture in the study area, respectively. In order to estimate agricultural production, total amount of rice and wheat crops was calculated based on potential for agricultural production per unit area in sub-watersheds (Ahmadi et al. 2018a, b). For estimating agricultural product loss, all rice and wheat farms located in high and very high flood risk areas were considered potential agriculture crop losses (Kourgialas and Karatzas 2016).
The amount of water loss associated with crop losses was estimated using WF (Eq. 3). This indicator represents direct or indirect use of water to produce goods or services and includes three types of WF: blue (using surface or groundwater), green (using rainwater), and grey (the amount of water to assimilate pollutants) (Hoekstra 2017). Flood damage to agricultural land represents loss of water, especially blue water used for, e.g., irrigation, which directly impacts future agricultural production. Therefore, in the present study, the amount of water footprint loss (WFL) was calculated based on water requirement of the crop (irrigation or rainfall) as: where PAP is the potential of agriculture production, A is the area (ha) under rice and wheat production, P is the potential production of each crop (ton/ha), A i is the wheat or rice area (ha) on land in high and very high flood risk classes, and WF i is the water footprint of wheat or rice (m 3 kg −1 ), as suggested by Mekonnen and Hoekstra, (2011).

Geographically weighted regression
Geographically weighted regression (GWR) was used to assess relationships between agricultural land uses (rice and wheat) and flood risk areas at sub-watershed scale. This model can explore the spatial relationship between dependent and independent variables considering no-stationarity properties of targeted phenomena (Stewart Fotheringham et al. 1996). The GWR equation is (Fotheringham et al. 1998): where ( u j , v j ) is the coordinates for location j, βi ( u j , v j ) is the local regression coefficient for independent variables χ i at location j, 0 (u j , v j ) is the intercept, ε j is an error term, and y j is the value of the dependent variable for the j th sample. Local R 2 in the GWR model is an indicator of how dependent and independent variables are fitted together. The higher the value of local R 2 , the lower the residual square sum, indicating a higher correlation (Wu et al. 2017). Table 3 shows the conditioning factors and sub-classification of flood risk zones and their proposed rating. The effects of the different conditioning factors were as follows:

Results and discussion
Flow accumulation A key factor in flood risk, as high value of this layer represents concentrated flow, and consequently higher flood risk.
Slope An important spatial factor for identification of flood risk through surface runoff velocity and vertical percolation (Rahmati et al. 2016). The highest slope is assigned the highest rating (Table 3). In this study, the slope was extracted in percent based on the digital elevation model. The slope factor influence on the water velocity and plays a major  1 3 role in flooding in highlands and reflect a constant threat in lowlands due to gentle slopes. The slope map was created in ArcGIS 10.5 to quantify topographic controls on flood conditions.
Land use Influences flooding and inundation level through altering infiltration, changing the relationship between groundwater and surface water, and debris flow (Kazakis et al. 2015;Areu-Rangel et al. 2019). Runoff and consequently flood conditions vary considerably under different LULC patterns. In the Talar watershed area, the main land uses were found to be rangeland and forest (Fig. 1).
Rainfall intensity Expressed using MFI, it ranged from 21.27 to 48.36, with higher values in the north of the watershed and lower values in central and southern areas (Fig. 2).
Geology Can amplify and extenuate the frequency and magnitude of floods (Kazakis et al. 2015). Some formations, such as Karst, significantly affect flood generation, so a lower rating is associated with higher infiltration capacity and fewer flood events.
Elevation Plays a vital role in controlling surface flow movement and flood depth ). In the Talar watershed, elevation varies between 213 to 4003 m asl, with the highest values in the south and the lowest in the north. Elevation was identified as the most important factor in flood risk in the study area (Table 3), as in other study areas (Tehrany et al. 2015;Ozkan and Tarhan 2016;Wang et al. 2017). The weighting for elevation was higher than for rainfall intensity and inferred intensification of rainfall with increasing elevation, confirming findings in previous studies (Kourgialas and Karatzas 2011).
In order to identify areas at risk of flooding in the Talar watershed, conditioning maps (Fig. 2) were combined in a single map (Fig. 3). This revealed areas with very high flood risk (VHFR) and high flood risk (HFR) in northern sub-watersheds. Overall, flood risk appeared to be more intense in the north, and around the river network in the south of the watershed. This part of the watershed mostly consists of forest, agriculture, and residential areas. Green areas and forests are reported to reduce flood risk, because of their infiltration Fig. 3 Flood risk map for the Talar watershed at sub-watershed scale capacity and evapotranspiration potential (Gaston et al. 2005;Smith et al. 2011;Warhurst et al. 2014). However, the actual effect can differ, as the results of this study showed. First, the design of vegetation cover might be inappropriate from the flood mitigation viewpoint, so that green lands are convex and higher than adjacent roads, increasing the flood risk, as reported for Shanghai, China (Wang et al. 2017). Second, the elevation factor, which was identified as the most influential conditioning factor in the Talar watershed, is low in forested and northern parts of the watershed. Some eastern, western, and southern parts of the watershed, which were found to have low and very low flood risk (LFR-VLFR) (Fig. 3), are mostly covered by rangeland and have high elevation.

Validation
In order to verify the final flood risk map, 134 historical flood points were used (Fig. 3). The positions of observed flooding points were determined through document obtained from regional water company of Mazandaran and field observation, interviews with local people, and water tail evaluation, and recorded by global positioning system (GPS) conducted in 2017. The overall accuracy of the flood risk map was assessed for two classes, HFR and VHFR. The location of historical floodplains in different flood level risk classes was determined by overlaying historical flood risk points on the flood risk map and using the spatial analyst tool (to extract values to points). A total of 111 out of 134 flooded points was found to be located in the high and very high flood risk classes in the final flood risk map, which means that the overall accuracy of the flood risk map was 83%.

Flood risk and damage
The area of agricultural land in flood risk zones was calculated for each sub-watershed. In order to compute the potential water loss, the area and yield of the rice and wheat crops had to be determined. Based on the land use map, an area of 2516.44 ha of the Talar watershed was occupied by rice and 31,260.71 ha by wheat (corresponding to 17% of watershed area). As the basin-based information on yield was not available in the study area, so we considered the average of province for our study. The average yield of wheat (1.73 ton ha −1 ) and rice (4.79 ton ha −1 ) in Mazandaran province (Ahmadi et al. 2018a) was assumed.
The area of rice and wheat farmland in high and very high flood risk areas (HFR-VHFR) in each sub-watershed is mapped in Fig. 4, and corresponding data are summarized in Table 4. For rice farmland, the highest flood risk level was observed for sub-watersheds 1 and 2, with area 369.79 ha and 310.69 ha, respectively. For wheat, the highest flood risk level was observed for sub-watersheds near the outlet. In general, a large proportion of land under rice (613.67 ha) and wheat (12,750.45 ha) was found to be exposed to VHFR flood risk damage and inundation in the 500 m buffer zone of the Talar river.
Assuming the flooded agricultural areas lost all cultivated crops and that all agricultural areas were productive, the potential losses of rice and wheat, derived by the intersection of land use and flood risk areas were calculated to be 7972.05 tons and 18,860.65 ton, respectively. The associated irrigation water loss through rice and wheat crops, calculated as WF, was 12.10 and 34.04 million m 3 (MCM), respectively, or over 46 MCM in total (Table 5 and Fig. 5).

Spatial relationship between areas of flood risk zones and crops
Spatial statistics can play a key role in environmental studies. Our GWR model results indicated that the value of local R 2 varied from 0.11 to 0.82 in the study area (Fig. 6). In the case of rice-HFR, the classes with high values of local R 2 were shown to be located in northern and central areas of the watershed, where the agricultural land use and residential areas are mostly located. For rice-VHFR, classes with medium and low values were shown to be located in the south of the watershed. In addition, the spatial distribution of local R 2 revealed a gradual decrease from north to south. In the case of wheat-HFR, the classes with high values of local R 2 were located in eastern and south-western parts of the watershed, while other parts of the watershed had low value of local R 2 (Fig. 6). Classes with high and medium values were located in the south of the watershed. These findings indicate that the flood potential map is sensitive to location, which is consistent with other studies (Sim et al. 2014;Wang et al. 2017;Purwaningsih et al. 2018). Flood risk zoning in agricultural watershed with intensive rainfall and floods has received attention from agricultural policymakers and planners because understanding of flooding conditions enhance their responses to flood risk. In current study, the efficiency of the water footprint analysis for spatial flood risk mapping and crop-water loss modeling was investigated and the resulting map was interpreted with the rice and wheat farms agricultural productions. Results demonstrated that the employed model had excellent performance in agricultural watershed.

Conclusions
Due to the special topographic and climatic conditions of Mazandaran province, agricultural lands are always in danger of annual flood damages. Rice farms are located near the river to provide cost-effective water supply, so these lands are faced with flood risk in comparison to other areas. As a result, assessment of spatial flood risk to cropland is important for better flood risk management and damage reduction in agricultural watersheds. In the present study, a flood risk map was prepared for the Talar watershed based on MCA method, one of the most important irrigated and rain-fed agricultural areas in Iran, where flooding occurs periodically. The map was created using six conditioning factors (flow accumulation, slope, land use, rainfall intensity, geology, and elevation). Of these,  5 Potential production and water footprint loss at sub-watershed scale elevation best explained frequent flood occurrence. Verification of the flood risk map using historical data on flood events revealed 83% accuracy for the map. The results showed that 31% of the total area of Talar watershed is under very high and high flood risk, and that 38% of its cropland (rice, wheat) is under high-very high flood risk. The irrigation water loss associated with crop losses to flooding was estimated to be 46.15 MCM, based on WF analysis. In this paper, in order to prepare flood risk map and to estimate agricultural watercrop losses, we integrated MCA modeling with water footprint method and agricultural statistics. Since the flood as a natural-human destructive phenomenon depends on various environmental factors, the use of multi-criteria analysis modeling reduces this complexity. Also, using the water footprint method, water loss can be calculated for agricultural crop losses caused by flood damage. In addition to agricultural products, water is also wasted, and this is an important issue in water resources management that has received less attention in flood damage studies. Also, in order to better understand the relation between flood