Assessment of the Sustainability of the Territories Affected by Gully Head Advancements through Aerial Photography and Modeling Estimations: A Case Study on Samal Watershed, Iran

Gully erosion is considered one of the major issues of land sustainability because it can remove considerable volumes of sediment and productive soils. Once started, gullies can continue to move by headcut retreat, or slumping of the side walls. Studies of gully development require constant monitoring activities which are not possible in not-well-explored areas, such as the arduous region of Iran, due to costs and a lack of geoinformation. Thus, the present research attempts to assess gully evolution using only two digital aerial photographs of different periods (1968 and 1994) and field assessment (2009) to estimate the gully head advancement based on frames geometry and rigorous procedure in southwestern Iran. Also, the gully head advancement was estimated and compared among them by different empirical equations. The results indicated that the mean of gully head advancement was 1.4 m year−1 and 1.2 m year−1 during 1968–1994 and 1994–2009, respectively, and the annual average of sediment mobilization was 26.8 m3 ha−1 in 2009. The model assessment indexes indicated that SCS (Soil Conservation Service) II was the best model for gully head advancement estimations in this study area. The main reasons for this can be associated with the Rp factor (previous gully head advancement) and the local environmental conditions. We conclude that the sustainability of the territory has been greatly affected due to this advancement. We also hypothesize that gully head changes could be related to the susceptibility of geological formations, climate, soil properties, and the coincidence of other gullies’ formation with common drainage networks in the study area. Based on the obtained results, land managers can use the results to distinguish the gullies in this region with a higher environmental risk, and to decide an effective implementation of soil conservation measures in order to include them in the land management plans.


Introduction
Soil erosion directly affects the sustainability of the territory because it enhances soil, water, and nutrient losses, decreasing fertility and damaging ecosystems [1,2]. The gullies are one of the most complicated linear erosion features, due to their negative impacts generated on agricultural, rural, or urban areas, such as flash floods or landslides [3,4].
The prediction of a gully's evolution is complicated, because they show different multistages, and are controlled by a large number of natural and anthropogenic factors [5][6][7]. Gully erosion is the most advanced stage of water erosion that also affects soil properties and depths, restricts specific land uses, and can threaten constructions [8,9]. Moreover, the gullies can reflect intermittent watercourses depending on previous soil water conditions and rainfall intensities [10], where the processes of channel incision due to erosion can be very intense and unpredictable [11]. Thus, due to these facts, the effects on sediment and water supply into the watersheds in the form of high and rapid peaks should be better understood [12][13][14]. This may explain, in spite of numerous efforts made in different parts of the world, an inclusive model to predict the initiation, growth, and evolution of gully erosion that has not been developed until now. Therefore, more research would be made to find possible solutions for the prediction and mitigation of these abovementioned negative effects using, for example, vegetation cover, mulches, or developing sustainable land management plans [15,16].
Fortunately, as a first step, the negative role that gully erosion plays in land degradation and sediment yield are being well-documented in different parts of the world in the last decades. In Western Europe, the sediment contribution of gully erosion has been measured to be from 30 to 80% [12]. In Australia, the contribution of gully erosion to the total soil erosion was estimated to be approximately about 50%. Also, there are some other study cases that have shown gully and bank erosion are the major sources of sediment and water losses, and the negative off-site effects [17]. For example, in Romania, gully head retreat has been measured up to 0.92 m year −1 [18]. In southern Poland on steep slopes, the mean gully advance rate was 0.21-0.52 m year −1 , and on the valley bottoms, from 0.18 to 1.98 m year −1 [19]. Some researchers report gully head retreat in the northeast of China from 6.2 m year −1 to 8.6 [20], Also, in India, changes were estimated to be 4-5 m year −1 [21] or, in Ethiopia, up to 0.34 m year −1 [22]. Regardless, it is important to indicate that there are different kinds of gully head retreat assessments, but they must be adapted depending on each specific environmental conditions and human impacts as it was recently confirmed in a review paper where the variations can oscillate between 0.01 m year −1 and 135.2 m year −1 , with median rates of 0.89 m year −1 , and an average of 5.0 m year −1 [17].
It is important to remark that soil erosion due to gullies can be found in different kinds of landscapes and under distinct climate conditions, which indicates that it is a very complex geomorphological process. However, in semi-and arid environments and, specifically, in non-developing or rapidly developing countries, which register an elevated concentration of extreme rainfall events and big extensions of bare soils due to the intensive agricultural practices and grazing, the negative impacts of gullies are much more severe [23].
Nowadays, the use of aerial photography and remote sensing techniques, in combination with geographic information systems (GIS) capabilities, are useful tools in the studies carried out in very specific environments, such as semi-and arid areas [24][25][26], and with high spatial variations [27,28]. The collection of historical data of gully's changes and advancements are indispensable to achieve gully studies; however, it requires annual monitoring, which is much too expensive for non-developing or rapidly developing countries, such as Iran. Therefore, the use of modelling techniques with GIS could support the smaller amount of aerial photography available in the country.
There are two general category models for gully growth; the first has an empirical basis, such as Thompson [29], Seginer [30], and SCS (Soil Conservation Service) [31] models, and the second a conceptual one, as well as physically based processes, such as EGEM (the Ephemeral Gully Erosion Model) [32], GULTEM (the model to Predict Gully Thermoerosion and Erosion) [33] and CHILD (the Channel-Hillslope Integrated Landscape Development) [34] models. According to the use of the empirical basis model in this research, it should be noted that the Thompson model was developed with data collected from some areas in the United States and, subsequently, in this model, there are two important factors to be considered: (1) rainfall, and (2) the clay content. The Soil Conservation Service of USDA during 1960s, based on 210 gullies in six areas of the Rocky Mountains in the United States, developed the first model of SCS for prediction of headcut retreat. In this model, unlike the Thompson model, soil properties are not considered for predictors of gully growth, because they are based on the total area and rainfall. Finally, the Soil Conservation Service added the R p factor (past average annual rate of gully head advance) and other two previous factors developed by the SCS II.
The pioneer study carried out by [35] can be applied for gully head advancement calculation in future years. This author, using two aerial photos of 1957 and 1968, showed a solution for the best accurate calibration of the SCS I model in Zahan-Ghaen watershed in Iran, where there is no dense vegetation. After that, [36] also calculated the average advancement of one gully of 0.73-0.1 m year −1 , based on digital information and comparing three different series of aerial photos in two regions of New Zealand. They stated that if the aerial photos have an appropriate resolution, they are able to be used for gully head advancement assessment. Decades after, many researchers used digital aerial photos with a high spatial accuracy to assess gully head retreat [37] and also in combination with SfM (structure from motion) techniques [38,39]. However, these studies did not insist on the sustainability of the territory. We consider that main goal of these kinds of results should be highlighted tools to develop sustainable land management plans in non-developing and rapidly developing countries.
One example where gully erosion is a big concern for rural and urban areas is Iran [40][41][42]. Thus, in this study, we tried to estimate gully head advancement through digitizing the minimum number of aerial photos between 1968-1994 and 2009, in combination with field observations in the Samal watershed in Iran. Moreover, we compared our results with some empirical models, such as Thompson, SCS I, and SCS II. Finally, the correlation between some environmental factors (topography, morphology factors of gully, and soil attributes) and gully head advancement was assessed to find the key factors on gully advancement. We hypothesize that this research allows us creating a useful database for land management planners in order to conserve the sustainability of this watershed.

Study Area
In order to conduct the purposed research, a part of the Samal watershed near the Samal village in Boushehr province was selected. Samal watershed is located between 51 • 8 to 51 • 25 E and 28 • 59 to 29 • 11 N ( Figure 1a) in Southwest of Iran. The total area of the Samal watershed has 29,750 ha, while that the selected study for gully measurement approximately covers 4500 ha in the northwest of Samal watershed (15.1%).
The geological map ( Figure 1b) shows that the catchment is located in the Zagros Mountain zone, which includes the Fars lithological groups (Aghajari and Mishan). In the highest areas, Bakhtiari conglomerates and, in the plains and steppes, quaternary alluvial materials and marls are the results of millennial years of erosion and sedimentation processes over the Cenozoic era, and mostly quaternary period [43]. These lithological groups are characterized by a high erodibility, therefore, the development of rills and gullies are very common in non-vegetated areas. In Figure 2, the two selected gullies can be observed, which are characterized by a claw-shape form ( Figure 2a) and U-shape morphology ( Figure 2b).
Based on the De Martonne aridity index, the climate patterns are characterized by desertic conditions [44]. In Figure 3, the mean annual temperature ( Figure 3a) and precipitation (Figure 3b) values during the last 30 years registered by the climate station of Boushehr (28 • 58 N, 50 • 49 E) are represented. We can observe a clear increase of the mean temperatures, which influences the rapid increase in weathering processes and evapotranspiration, and the decrease in vegetation cover and agricultural and livestock practices of this territory [45]. The annual average rainfall is about 180 mm.

Research Methodology
In Figure 4, a flow chart with the general context of the research methodology is included. a b Figure 3. Temperature (a) and precipitation (b) changes during the last 30 years in the studied area.

Research Methodology
In Figure 4, a flow chart with the general context of the research methodology is included.

Research Methodology
In Figure 4, a flow chart with the general context of the research methodology is included. a b Figure 4. Flowchart of the methodology used in this study.

Gully Mapping by Using Aerial Photos and Soil Sampling Procedures
Based on the visual interpretation of the aerial photos and the assistance of Google Earth data, the main gully was mapped using the software PCI Geomatica 1.9 (PCI Geomatics, Gatineau, QC, Canada) and ArcMap 10.5 (ESRI, Redlands, CA, USA). Mapping was carried out on 1.12 × 1.12 m orthophotos from 1968 to 2009. Linear features, the presence of headcuts, and steep side walls were the main criteria to be distinguished. Moreover, field observations were needed in order to confirm that the delimitation was correctly carried out. In the largest permanent gully selected for this research, soil samples into the headcuts and side walls were taken. It is important to remark the limitation of this kind of analysis based on the difficulties collecting a high number and frame of aerial photos, and the costs and time invested to make field campaigns in a not-well-connected arduous region in the southwest of Iran. They were transported to the laboratory, and soil texture (sand, silt, and clay), SAR (sodium adsorption ratio), pH, CEC (cation exchange capacity), OM (organic matter), and limes were analyzed following standardized methods [46,47]. We used the standard hydrometrics procedures for soil texture estimations [48] and a LECO (Laboratory Equipment Corporation) furnace to determine soil organic carbon by combustion [48]. We also determined CEC when washing the soil samples to remove the excess of salts, and using an index ion to estimate the total positive charge in relation to the original soil mass. This involves bringing the soil to predetermined pH values prior to conducting the analysis [48]. The pH value is a measure in a water solution (1/5) [48]. Finally, we used following formula [48] for SAR obtention: in which all cation concentrations are expressed in mmol/L.

Mapping Geomorphological Features
Digital topographic layers, such as contour lines, point elevations, rock outcrops, ephemeral water bodies, and streamlines (1:25,000), were used to extract the digital elevation model (DEM) through the Iranian National Cartographic Center (NCC) and DGN (cartographic map of Iran) data.
A DEM was used for the implementation of morphometric characteristics of this study area, such as (1) A t (gully outlet upslope area); (2) A 1 (headcut upslope area); (3) L (distance from the gully head to the boundary of the drainage basin), and (4) H (height of gully head).The spatial resolution used for this study was a 7 × 7 m cell size, and all the maps were obtained using the PCI Geomatica 9.1 (PCI Geomatics, Gatineau, QC, Canada). The image analyses were performed as follows: (i) data entry (images and ground data); (ii) collection of ground control points; (iii) solving the mathematical models; and, (iv) orthophotography analysis generation.

Orthophotography Preparation for Image Analysis in 1968 and 1994
A scanned diapositive of multitemporal aerial photos (1968 and 1994; Figure 5) was digitized through a rigorous method [49]. It is characterized by a mathematical approach using the software PCI Geomatica 9.1 (PCI Geomatics, Gatineau, QC, Canada) based on the detection of fiducial marks, a camera calibration, and geometry of the photo frame. This valuable dataset comes from the Iranian National Cartographic Center (NCC). Then, a mosaic of aerial photos (1 × 1 m with RMSE = 0.5 pixel) was created and, finally, the spatial distribution of the gully was mapped. Moreover, their associated contribution catchment boundary on each gully head was extracted through a high-resolution DEM (flow accumulation tool), and was rechecked and corrected based on the high accuracy mosaic orthophoto.

Measurement of Gully Head Advancement in 2009
After selecting the gullies (Figure 1a), detectable with the aerial photos of 1968 and 1994, the gully advancement between 1994 and 2009 was directly calculated performing measurement points along the borders with a GPS (Etrex 30, GARMIN, Olathe, KS, USA). The data positions were summarized using the ArcMap 10.5 (ESRI, Redlands, CA, USA) software. Also employed was the highly accurate GPS using a UTM (Universal Transverse Mercator) coordinate system to locate the headcuts for 2009. However, for the other previous periods, we had to use the aerial photos, paying attention to the coincidence with the roads and the topographic map (DGN).

Procedures of Gully Modeling to Estimate Head Advancement
The gully head advancement was estimated by using the Thompson [29], SCS I [31], and SCS II [31] models: (1) Thompson. R: Total gully head advancement during period of P calculation (m). A: drainage area above the gully head (ha). S: slope of approach channel above the gully head (%). P: total rainfall for 24 h that is equal or greater than 12.7 mm. E: percentage of clay content of soil profile through which the gully head is advancing.
(2) SCS I. R: rate of gully head advancement (m year −1 ), A: drainage area above the headcut (ha), P: total of 24-h rainfall events of 12.7 mm or greater, occurred during the gully's development (mm).
(3) SCS II. R: rate of gully head advancement (m year −1 ), Rp: past average annual rate of gully head advancement (m year −1 ), A1 and A2: ratio of the drainage area of a given upstream which can reach to the drainage area through which the gully has moved, P2 and P1: ratio of the expected long-term average annual inches of rain from 24-h rainfall events of 12.7 mm or greater to the average annual inches of rain from 24-h rainfall events of 12.7 (mm) or even greater for the studied period. * Since the Rp index is needed in this method, the model was applied only for the last period (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)

Measurement of Gully Head Advancement in 2009
After selecting the gullies (Figure 1a), detectable with the aerial photos of 1968 and 1994, the gully advancement between 1994 and 2009 was directly calculated performing measurement points along the borders with a GPS (Etrex 30, GARMIN, Olathe, KS, USA). The data positions were summarized using the ArcMap 10.5 (ESRI, Redlands, CA, USA) software. Also employed was the highly accurate GPS using a UTM (Universal Transverse Mercator) coordinate system to locate the headcuts for 2009. However, for the other previous periods, we had to use the aerial photos, paying attention to the coincidence with the roads and the topographic map (DGN).

Procedures of Gully Modeling to Estimate Head Advancement
The gully head advancement was estimated by using the Thompson [29], SCS I [31], and SCS II [31] models: Thompson. R: Total gully head advancement during period of P calculation (m). A: drainage area above the gully head (ha). S: slope of approach channel above the gully head (%). P: total rainfall for 24 h that is equal or greater than 12.7 mm. E: percentage of clay content of soil profile through which the gully head is advancing.
SCS I. R: rate of gully head advancement (m year −1 ), A: drainage area above the headcut (ha), P: total of 24-h rainfall events of 12.7 mm or greater, occurred during the gully's development (mm).
SCS II. R: rate of gully head advancement (m year −1 ), R p : past average annual rate of gully head advancement (m year −1 ), A 1 and A 2 : ratio of the drainage area of a given upstream which can reach to the drainage area through which the gully has moved, P 2 and P 1 : ratio of the expected long-term average annual inches of rain from 24-h rainfall events of 12.7 mm or greater to the average annual inches of rain from 24-h rainfall events of 12.7 (mm) or even greater for the studied period. * Since the R p index is needed in this method, the model was applied only for the last period (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)).
Finally, the results obtained from each model were compared to the real measured gully head advancements in each aerial photography (1968 and 1994) and GPS points (2009). In this regard, the coefficient of determination (R 2 ) and the mean square error (MSE) indexes were used to check the validation of each model.
where n is the number of data points, f i the value returned by the model and y i the actual value for data point i.

Determining Which Key Factors Affect the Gully Head Advancement
We used the regression methods because they are the most widely used statistical tools for data analysis. In this case, several response variables are studied simultaneously, thus, we are in the sphere of multivariate regression. The usual description of the multivariate regression model, that relates the set of m multiple responses to a set of n predictor variables, assumes, implicitly, that the m × n regression coefficient matrix is of full rank. Therefore, a multivariate linear regression procedure was performed to examine the relationship between the mean of gully head advancement, topographical variables, and soil properties, to determine the major factors affecting the gully head advancement. Only the best adjustments were added to the final results. The software SPSS 16 (IBM, Armonk, NY, USA) was used to carry out these analyses.

Soil Properties and Headcut Advancement Measurements
In Table 1, soil properties and vertical height were presented. The soils are characterized by a high content of sands (≈70%) with maximum values of 91%. The average clay content is 10.9%. As several previous research studies indicated, sand and clay fractions, f, are one of the main soil physical attributes for gully and rill development. These are the cases in the well-studied badlands [50,51], where the highest soil erosion rates have been quantified. However, our findings did not show a clear relationship between them, which can be addressed into the effect of other specific dominant factors, such as the runoff discharge, seepage, and mass movement capacity. The mean SAR values are 2.3, showing a maximum of 11.3. The mean value of CEC is 4.8 cmolc kg −1 , with 7.8 cmolc kg −1 being the maximum. pH values register values for non-acidic soils (7.7), and the percentage of limes reaches 64.6%. Moreover, it is important to highlight the extremely low content of organic matter, which reaches 0.22%. In southern Italy, one research [52] carried out a geomorphological analysis, and analyzed the relationship between topography and organic matter content, concluding that a clear positive correlation exists between the low organic matter content and high soil erosion rates. As we can observe in the aerial photography, the high soil erosion rates could have affected the presence of vegetation cover, which is a clear signal of the registered extreme low content of OM. Obviously, the absence of vegetation and the low organic matter content are highly correlated among them, which was also confirmed in Spanish areas with gullies [4]. Finally, we summarized mean height values of 1.37 m, reaching a maximum of 2.7 m and minimum of 0.5 m.
In Table 2, the results of measuring gully headcut retreat over the different periods are shown. Between 1968 and 1994, a higher gully's growth than from 1994 to 2009 is registered, reaching total average values of 1.36 and 1.23 m year −1 , respectively. However, a higher variability during the most recent studied period is noted, obtaining a higher standard deviation (0.64 vs 0.67 m year −1 ). The elevated variability of each gully allows us to confirm the statement made by some studies [28,53], which concluded that gully development, for some years, has been highly unpredictable. During the first studied period

Gullies' Head Advancement Measurements Using the Modeling Techniques
The estimated growth using the three selected empirical models is presented in Table 3. In the Supplementary Material, some parameters used for each estimation were also added.
In general, all the models underestimated the mean values, except the SCS II model. SCS I model registers a total average change of 0. 47   Other researchers in different parts of the world also applied these indexes to different gullies, obtaining a diverse rate of changes. For example, it obtained 0.206 m year −1 in Habl-e-Rood (Iran) [54]. Also, it registered changes from 0.73 to 0.1 m year −1 in New Zealand [36]. Another example can be found in the research of [19] in Poland, where a rate of 0.63 m year −1 was obtained. Also, Martinez-Casasnovas [55] recorded some changes in a gully headcut in Spain. General values showed a rate of change of 0.2 m year −1 . However, further research carried out at the same place demonstrated maximum values from 0.7 to 0.8 m year −1 . These researchers also used aerial photographs to measure gully retreat, and we can conclude that, after comparing to the model results, the gully head advancement was up to 3 times higher in some cases. The authors mentioned that the main differences may be explained because of the sensitive soil properties, that the models obviated the 24-h rainfall events with a high intensity and the interconnection drainage network, which is not included in the model.

3.3.Comparison between the Model with Real Observed Data
In order to select the best model, the real observed data on the aerial photos are compared to the model estimations using assessment indexes (Table 4). Comparing the calculation and model prediction, we demonstrate that among the used models, the SCS II model obtains the highest accuracy (R 2 = 0.19, R = 0.43) and the lowest error rate (MES = 0.16), being the best model to predict the gully head advancement in our study area. The main reasons to explain the high accuracy for SCS II model may be related to the importance of the R p factor (gully head advancement in previous years) that was not considered in the other models. This factor is able to involve the regional environmental effects of soil properties, land uses and geomorphological changes, which are neglected by the other variables in the rest of the tested models. These results are similar to the findings of [54] in Habl-e-Rood watershed (MAE, Mean Absolute Error= 37.3%) and [35] in Zahan-Ghaen watershed (RMSE, Root Mean Square Error= 0.895) which also considered the SCS II model as the best estimator of gully head advancement. Moreover, the comparison of the R2 during the two studied periods reveals that the accuracy of both SCS I and Thompson models for long-and short-scales is not similar. In fact, over a longer time span, more extreme rainfall events, and consequently, a higher gully retreat would be expected. Therefore, none of that mentioned models can be suitable for estimating short-term gully head advancements. The lower R 2 for SCS I and Thompson may be related to their considered factors. The correlation analysis revealed that SAR, contribution area, and headcut depth can be confirmed as the main effective factors for predicting gully head advancement. However, in the mentioned models, just the contribution area is considered and, therefore, the ability of these models (R 2 ) has been subjected to uncertainty. One of the main factors concerned, regarding the low accuracy of these models, may relate to the effect of soil chemical and physical attributes of gully erosion. None of the models has considered the effects of soil chemical attributes directly, while the results of regression analysis will reveal the effect of other factors than components of mentioned models.

Key Factors That Hide the Sustainability of the Territory by Gully Development
The most adjusted results of the multivariate linear regression between the mean of gully head advancement, topographical factors, and soil properties, are presented in Table 5.
The performed analysis shows that there are significant positive correlations between the gully upstream watershed area (A t and A 1 ), the distance from the gully head to the boundary of the drainage basin (L), the height of the gully head (H), the SAR, and the gully head advancement.
Related to the soil properties, we can confirm that sodium causes the dispersion of fine particles, especially of clay in the soil, and plays an important role in the sensitivity of the loss of resistance of the soil against the water impact [56].
The assessment of the role of various topographical changes shows that there is a significant correlation between these variables and gully erosion advancement. The influence of topography on soil properties and erodibility is coincident with other findings highlighted by previous researchers. For example, one researcher demonstrated [57], in a river close to the Three Gorges Reservoir in China, that the topography can change water quality characteristics and, subsequently, the intensity of the hydrological processes, such as the erodibility. One of these changes due to the topography is the increase in sodium and acidification of the soils, also coincident with our findings. Moreover, these impacts can be more extreme if both factors interact with land use changes due to agricultural practices and grazing [58]. The gully upstream watershed area, the distance from the gully head to the watershed ridge, the height of the gully head, and the gully head advancement, can explain about 85% of the total variance. These small positive correlations between the most mentioned variables (Table 5) indicate that the most important factors on gully extension are also related to the runoff generation (e.g., A t , A 1 , and L) and the association with the kinetic energy (H) received by gully's headcut. Also, this interconnected dynamic between runoff and soil loss was clearly confirmed using rainfall simulations and erosion plots by several authors in other arid environments [51,59,60]. Thus, the conservation of the vegetation cover, the application of mulches, or the regulation of agriculture and livestock activities should be a mandatory decision to preserve the sustainability [61,62].
Therefore, as the watershed area and distance from the watershed ridges to the gully headcut increase, consequently, more runoff and hydrological responses will be available on the gully's headcut and the gully growth can also increase [12,54]. Therefore, for future model developments of headcut retreats, anthropogenic factors should be considered, such as land use changes, grazing impact, or the elimination of vegetation cover. Also, for each time span (temporal scale), spatial factors should be considered. In this regard, it seems that use of the upstream watershed area of the gully is the major factor mentioned by previous works [17,37,63]. It was also found that gully erosion volume was highly correlated with these factors, which is similar to our findings.

Conclusions
Gully head advancement is a serious problem that must be considered as a big concern by the Iranian government and land management planners. Gullies can destroy infrastructure, decrease agricultural activities by loss of soil fertility, and hide transport facilities. Therefore, having knowledge about the rate of changes of gully headcuts could be considered as a useful tool to predict future gully development, in order to manage soil conservation policies. Our results indicated that gully advancement rates in arid lands, such as the Samal watershed, are dramatically higher than in other climate conditions, which can be related to sparse vegetation, soil susceptibility, and rainfall patterns. Our findings also revealed that over the next year-as the gully head advancement will continue-more sediment amounts will be removed by the gully headcut encroachment. In addition, the effects of future global warming over the following decades may change the hydrologic responses, and should be considered more in further research. Furthermore, although the SCS II model showed the highest coefficient (MES = 0.16), none of the Thompson, SCS I, and SCS II models are able to accurately predict the longitudinal gullies' advancements in this area. We conclude that this accuracy would be related to the R p factor. It seems that this coefficient is a good replacement for considering soil erodibility, slope, and land use changes on the upslope contribution area. Most of predicting models for gully developments are based on environmental characteristics of the regions in which they were generated, and additional data under various environmental (geo and climate) conditions could improve our understanding. Also, the comparison of the model accuracy in two different time period scales shows that none of the models are suitable for gully head prediction in short-time periods. Further multiple correlation analysis between gully head and environmental factors could reveal that five variables (A t , A 1 , L, H, SAR) can explain 85% of the longitudinal advancement of gully head changes.
For the future, to conserve the sustainability, soil conservation measures developed by land managers should be mainly focused on the abovementioned criteria. Our findings revealed that for precise prediction of the future geomorphological and hydrological dynamics of a gully, specific topographic and soil attributes should be considered. Moreover, in using this combined method of high-resolution aerial photos and field work, time and economic resources could be also saved and reused to develop efficient nature-based solutions. Finally, we hypothesize that this research allows us to create a useful tool for land management plan developments, in order to conserve the sustainability of semi-arid and arid territories.