Geomorphic threshold conditions for gully erosion in Southwestern Iran (Boushehr-Samal watershed)

Journal of Asian Earth Sciences 35 (2009) 180–189

Contents lists available at ScienceDirect

Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jaes

Geomorphic threshold conditions for gully erosion in Southwestern Iran (Boushehr-Samal watershed) Aliakbar Nazari Samani a,*, Hassan Ahmadi b, Mohammad Jafari b, Guy Boggs c, Jamal Ghoddousi d, Arash Malekian e a

Faculty of Natural Resources, University of Tehran, Karaj, P.O. Box 31585-3314, Iran Faculty of Natural Resources, University of Tehran, Karaj, Iran c School of Science and Primary Industries, Charles Darwin University, Darwin, NT 0909, Australia d Soil Conservation and Watershed Management Research Institute, Tehran, Iran e International Desert Research Center, University of Tehran, Tehran, Iran b

a r t i c l e

i n f o

Article history: Received 31 December 2007 Received in revised form 25 January 2009 Accepted 2 February 2009

Keywords: Gullying processes Erosion Area–slope Landuse impact Threshold Topography Soil attributes Iran

a b s t r a c t Globally, a large amount of research has been dedicated to furthering our understanding of the factors and mechanisms affecting gully erosion. However, despite the importance of gully erosion in arid and semi arid regions of Iran there has been no comprehensive study of the geomorphic threshold conditions and factors influencing gully initiation. The aim of this article is to investigate the gullying processes and threshold conditions of permanent gullies in an arid region of Iran based upon examination of the slope– area (S = aAb) relationship. The data were collected through field and laboratory studies as well as Digital Elevation Model (DEM) analyses. In total, 97 active headcuts were identified across the three study sites and classified based on dominant initiation process including piping, landsliding and overland flow. Soil properties, including EC, SAR and soil texture, as well as landuse practices were found to be the major factors initiating piping and bank gullies. All gullies initiated by landsliding and seepage processes were found to be located in steep areas (28–40% slope) with their distribution further influenced by the lithology and presence of a cohesionless sand layer within the soil profile. An inverse relationship between upslope area (A) and local slope (S), in which the a and b coefficients varied, was further investigated based on the dominant gullying process and land use. Gullies occurring in the rangelands that were dominated by overland flow had the strongest relationship while landsliding dominated gullies did not have a statistically significant S–A relationship. In comparison to theoretical and literature based relationships for gully initiation, relatively low values for b were obtained (0.182 to 0.266), possibly influenced by the presence of seepage and subsurface processes in many gullies. However, this is consistent with other studies in arid regions and may reflect greater potential for gullying in arid zones due to low vegetation cover and high variation in rainfall. In addition, the soil attributes together with land use practices influenced gully initiation thresholds. Application of the solved S–A relation for predicting vulnerable areas to gullying indicates that it is possible to predict the location of gullies with an acceptable level of accuracy; however other environmental factors should be integrated with the S–A relationship to more accurately identify the location of permanent gullies in arid regions. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Gully erosion is recognised as a major land degradation issue, causing both impacts on-site, through direct soil loss and off-site, through sediment deposition in downstream environments. Gullies are formed by water erosion and consist of several characteristics including a steep incised channel with an active headcut, unstable side wall and temporary water flow (Nordstrom, 1988; Poesen et al., 2003). * Corresponding author. Tel.: +98 261 2249313; fax: +98 261 2227765. E-mail address: [email protected] (A. Nazari Samani). 1367-9120/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jseaes.2009.02.004

A great deal of research has been undertaken on the contribution of gully erosion to overall soil loss and sediment production in a wide range of environmental and climatic conditions and at a variety of temporal and spatial scales. Although the relative importance of gully erosion has been well documented in the literature (Patton and Schumm, 1975; Vandaele et al., 1996; Poesen et al., 2003; Sidorchuk, 1999; Wasson et al., 2002), most soil loss equations and erosion models do not include the soil loss caused by gully erosion (Poesen et al., 2003). More than 60% of Iran’s area is located in arid and semi arid regions, with about 100 million ha at high risk of desertification (Ahmadi, 2004). The spatial and temporal variability of rainfall,

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

prolonged drought periods over past decades and intensive exploitation of agricultural lands as well as overgrazing of rangelands are the main causes of the desertification trend in Iran, leading to increased soil erosion and deterioration of ecosystems (Forest and Rangeland Organization, 2004). The vast area of Iran (1.6 million km2) combined with the presence of both wind and water as major erosion processes means that there is no particular and reliable soil erosion rate at a national level. However, based on estimation of 137 Cs measurements and suspended sediment gauging data for seven major dam watersheds of Iran, the water erosion rate in agricultural lands varies between 7.6 and 32 ton ha1 Y1 and 4.3– 22 ton ha1 Y1 in rangelands (cited by Hakimkhani, 2006). These large variations have been attributed to the wide range of environmental characteristics across Iran. Modelling of soil erosion in Iran has been based on empirical prediction models (e.g. PSIAC1, EPM2, and USLE3) considering the uncertainty subjected to such methods. Much more research is therefore required to understand the role of gully erosion in Iran. Gully erosion is particularly important in the Boushehr province located in Southwestern Iran (the study area of this research), which contains the biggest gas field of the world (Asalouyeh). In this region, gullies can undercut gas pipe lines and other related infrastructure, commonly causing them to subside and major damage to occur. Gully erosion is regarded as a threshold phenomenon controlled by a wide range of factors (Patton and Schumm 1975; Vandaele et al., 1996; Bull and Kirkby, 2002; Poesen et al., 2003; Valentin et al., 2005). The importance of threshold conditions is different in various climates and under different landuses, soil and land cover conditions. However, in different regions, researchers have found that there is a common inverse relationship between the upslope drainage area (A) and local slope gradient (S) of gully head locations ðS ¼ aAb Þ, as first defined by Patton and Schumm (1975). They stated that channel initiation is based on the kinetic energy of concentrated overland flow, which in turn is a function of runoff and slope (where upslope drainage area is used as a surrogate for runoff volume (Leopold et al., 1964)). This relationship has subsequently been applied to examine gully formation under a variety of conditions around the world, with regional differences being expressed in different values for the constant and exponent of the model (a and b, respectively) (Begin and Schumm, 1979; Montgomery and Dietrich, 1988, 1994; Prosser and Abernethy, 1996; Vandaele et al., 1996; Poesen et al., 1998; Vandekerckhove et al., 2000a,b; Desmet et al., 1999; Morgan and Mngomezulu, 2003). Comparison of the results of existing research indicates that variation of a and, b factors is independent, but depends on local environment characteristics and data collection methods. Furthermore, for most studies, values for the exponent b are more constant and vary between 0.2 and 0.4, while the constant a is more variable (0.0035–0.35). Variations in the value of b are commonly linked to the dominant process operating in the catchment. Values less than 0.2 are associated with overland flow while values greater than 0.2 commonly indicate that subsurface processes and mass movement are dominant (Morgan, 2005). Vandekerckhove et al. (2000b) showed that the trend of the threshold line is often clearer when data from mass movement dominated gullies is separated from hydraulic erosion dominance (changing b as slope of threshold line), while the constant ‘a’ varies by changing environmental attributes. Moreover, Montgomery and Dietrich (1994) consider AS2 as a suitable indicator for gully initiation and found a range in values of between 500 and 4000 m2, but Cheng et al. (2007) reported much smaller values (41–814 m2) and suggested that

1 2 3

Pacific Southwest Interagency Committee. Erosion Potential Method. Universal Soil Loss Equation.

181

AS2 does vary among regions and needs to be checked before application to an area. Although considerable research has been carried out on gullying processes and thresholds, much of this has been in humid to semi arid regions and has mostly focussed on ephemeral gullies in agricultural lands (Poesen et al., 2003; Vandekerckhove et al. 2000a,b). However, less effort has been applied to arid regions, particularly the rangelands of western and southern Asia. In Iran, Ghoddousi (2002) used environmental characteristics to predict the morphology of gullies through correlation analysis, while Adelpour et al. (2004) applied an experimental flume in the field to identify critical shear stress for channel initiation. However, no studies having been conducted to examine the geomorphic threshold conditions of gully initiation. This study therefore investigates gullying processes and threshold conditions for permanent gullies in a region of arid rangelands in Iran based upon examination of the slope– area (S = aAb) relationship, field study of dominant erosive processes and examination of associated environmental conditions, including soil characteristics and land use conditions. 2. Study area The study area is located in the Dareh-koreh watershed of Boushehr province, Southwestern Iran, at 29° 080 1800 N and 51° 130 1500 E (Fig. 1). The area consists of hills and alluvial plains with sedimentary geology including Gori limestone, Aghajari marl, Bakhtiyari conglomerate and Quaternary alluvium based on 1:100,000 geological maps. However, field observations indicated that some areas mapped as conglomerate, in which gullies were located, were actually covered by a hydro-aeolian Quaternary deposit. This is primarily due to the broad scale of the geological map (1:100,000, 1954). Also, these maps were generated by the National Petroleum Organization of Iran in 1954 and were therefore primarily developed for distinguishing formations containing crude oil resources. Soils are primarily Entisols (based on the American Soil Taxonomy) and are in the primary stage of development and limited in depth. Most profiles also identified as continuous sand layers approximately 20–35 cm thick within 60 cm of the surface. The study area is located in an arid environment with an annual average rainfall of 150 mm with mean annual temperature of 14.5°C, and relative air humidity of 52%. Generally about 80% of precipitation falls in 2–3 intense storm events, predominantly towards the end of autumn and winter, with a high temporal and spatial variability typical of such arid regions. For the purpose of this study three sites were selected for the analysis of gully head positions. The main characteristics of the sites are summarized in Table 1.

3. Materials and methods 3.1. Field measurements Topographic maps (1:20,000), recent aerial photographs (1:20,000) (1993) and geological maps were obtained from the National Cartographic Center (NCC) and Geological Survey Organization of Iran. Permanent gullies were identified by visual interpretation of the aerial photographs and examination in the field, with gullies greater than 0.7 m in depth being used for further analysis. The gully positions and all drainage channel(s) were then mapped using a portable GPS with a positional accuracy of approximately 3 m. To map gully locations and headcuts all sites were systematically walked from outlet to catchment divide and a series of points were marked using the GPS through all drainage networks. Data describing gully activity was also collected at this stage, based on pre-defined criteria as follows; abrupt and steep morphology, sharp edge, active plunge-pool and evidence of

182

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

Fig. 1. Location of the study area and sites of gully erosion.

Table 1 Main characteristics of the study area. Site number

Mean altitude a.s.l

Landuse

Lithology

Soil texture

EC ds/m

Vegetation cover

Geomorphology

Slope %

1 2 3

87 72 40

Rangeland Rangeland cereal crop Abounded land

Quaternary alluvial Quaternary alluvial Quaternary alluvial

Sandy loam Sandy loam Clay loam

2.6 1.9 40.8

Sparse annual grass Sparse annual grass Bare land

Hill Hill Alluvial pediment

23 20 1

tension cracks, flow marks as well as lack of vegetation cover on the bed, side wall and headcut of the gully (Dietrich and Dunne, 1993; Oostwoud Wijdenes and Bryan, 1994; Oostwoud Wijdenes et al., 2000). All field map data was then imported into a GIS environment and the gully networks were created and overlayed on the aerial photographs for validation purposes. Field observations were made regarding the characteristics of each gully’s upslope region, including the dominant land use (rangelands or abandoned land), and gully head slope, measured using an optical clinometer. A soil profile was also dug in each site and samples were taken from the different horizons to determine soil attributes. All of the soil samples were analysed in the laboratory to determine the chemical and physical attributes including: soil texture, Electrical Conductivity (EC), Na+, Ca++, Mg++, K+, Cl-, HCO 3 , lime, Sodium Absorption Ratio (SAR) and organic carbon through standard methods. Detailed rainfall records were not available for the area; however, because of the proximity of the sites we assumed that there is no significant difference in rainfall. For each gully, the dominant process of gully initiation was identified based on morphology and pre-defined criteria. Gullies dominated by landsliding were identified as having a heap of soil collapsed at the foot of the gully head, often with some ground vegetation cover (Fig. 2), while in gullies dominated by piping a hole and signs of sand boiling at the base of the gully head as well as in the side walls was observed (Fig. 3). Finally the presence of active plunge-pools was used to identify overland flow dominance (Fig. 4). Gullies were then separated based on the initiation process as these have been shown to exhibit different area–slope relationships (Montgomery and Dietrich, 1994; Dietrich et al., 1992; Vandaele et al., 1996; Vandekerckhove et al., 2000a,b).

to be calculated. Previous studies have used measurements based on aerial photographs, field optical clinometers, topographic maps and Digital Elevation Models (DEM) (Patton and Schumm, 1975; Vandaele et al., 1996; Desmet et al., 1999; Vandekerckhove et al., 2000a,b). However, the S–A relationship has been shown to be sensitive to the methods of data collection (Vandaele et al., 1996; Desmet et al., 1999). A DEM was constructed by interpolating 10 m interval digital contour lines (1:20,000, DGN) and elevation points provided by the NCC. The grid cell size of the DEM was 5 m. All processes were undertaken using the ESRI ArcGIS 9.1 software package and its extensions. A slope map was calculated from the DEM using the method introduced by Zevenbergen and Thorne (1987) and the

3.2. Area–slope analysis b The area–slope relationship ðS ¼ aA Þ requires the upslope contributing area (ha) and slope value (m/m) of each gully head

Fig. 2. Gully development associated with landsliding indicated by the heap of soil collapsed with the same ground vegetation cover at the lower part of the gully head.

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

Fig. 3. Headcut developed by piping. Note the hole in the gully wall associated with seepage at the lower part of the gully head.

Fig. 4. Abrupt gully head associated with overland flow, indicated by a plunge-pool located at the base of the gully head.

slope at each gully head extracted in the GIS. The upslope drainage area was assessed using the Arc Hydro extension of ArcGIS (Grise et al., 2002). Finally, all data describing the slope and contribution area for each gully was extracted in the GIS and transferred to SPSS for statistical analyses. To estimate values for the constant a and exponent b in the S–A relationship, the data were plotted on a log–log scale and according to power regression the best fitted line was selected. 4. Results and discussion 4.1. Gullying process and characteristics Totally, 97 active headcuts were identified across the three study sites, with 63 of these located in site one, 16 in site two and 24 in site three (Table 2). Although the number of gullies identified in each site is different, their density per unit area is comparatively similar. All identified gullies were found to be continuous and occurred in the drainage stream channel with sharp headcuts located at the starting point of tributaries. In site one and two overland flow and landsliding were found to be the dominant processes of gully initiation, while piping was identified as a major process in site three. This difference was found to be related to the topogra-

183

phy of the three watersheds, with site one and two dominated by higher slopes and site three characterised by low slopes (Table 2). That is, as slope increases the frequency of gully initiation by landsliding is generally increased. Although, channel evolution modelling by Montgomery and Dietrich (1994) indicates that the threshold slope for channel heads initiated by landslides is greater than 45%, we found many channel heads initiated by landsliding in areas with 28–40% slope. This result might be related to the presence of both seepage and landsliding processes in some gullies. The geology of sites one and two was found to be dominated by hydro-aeolian materials of quaternary formation within a hilly geomorphology (Fig. 5). These soils were characterised by silty/fine sand textured soils (57%), with bands of coarse sand of 20–35 cm thick occurring within 60 cm of the surface. The existence of such soil profiles in this area may also explain why landsliding has taken place in lower slopes when compared to Montgomery and Dietrich (1994). That is, the nature of materials at depth, particularly the presence of different soil strengths and cohesionless sand layers, together with drying and wetting cycles, creates tension or desiccation cracks during dry seasons consequently collapsing the top layers. Poesen et al. (2002) stated that in dryland areas two types of mass failure can be identified, one is continuous failure over long periods and the other is catastrophic shear failure in cohesive layers. In site one and two the latter, in combination with seepage, were identified as the most frequent processes of gully initiation and development. In site three the situation is different, with the gullies primarily being bank gullies. Bank gullies are frequent features in arid and semi arid regions (Poesen et al., 2002) and occur in river banks or terraces and retreat into gentle slope pediments, bank terraces or arable lands. In contrast to hillslope gullies that are formed due to critical flow shear stress at the soil surface, they are formed because of runoff overfalling from littoral banks and tunnel or piping processes caused by hydraulic gradients in dispersive materials (Poesen et al., 2002). As Table 2 indicates, all gullies in site three were identified as being initiated by piping processes. Table 3 shows the results of the mean value of measured soil parameters (EC, SAR and ESP4) relevant to piping for the samples taken in each soil profile. As shown in Table 3, the site three has considerably higher EC (40.86 dS/m) and SAR (48.2) values when compared with sites one and two, indicating that the soils are highly sodic (SAR 20–35) and saline (EC 10–26 dS/m). Sodic and saline soils are vulnerable to erosion, as the presence of high levels of sodium in the clay fraction causes dispersion and deflocculation of these materials, and consequently the subsurface flow creates tunnels and pipes (Gutierrez et al., 1988; Poesen et al., 2002; Romero Diaz et al., 2007). Studies by Harvey (1982) and Martın-Penela (1994) have also shown that the nature of materials at depth and presence of fragile silty-clay material containing cracks and discontinuities joints and faults filled by gypsum can cause differential infiltration in the soil profile and lead to piping. Although we were unable to find the intake area to study soil surface properties in the site three as all soils were disturbed by agricultural activities, samples from deeper layers in the soil profile (>70 cm) showed that the bottom material is clay loam with high levels of gypsum (41 meq/100 g of soil) and therefore increased vulnerability to piping. These results are consistent with previous research carried out by Gutierrez et al. (1988), Poesen et al. (2002), Romero Diaz et al. (2007), Faulkner et al. (2000). The presence of agricultural activities in site three is also likely to increase the risk of erosion caused by piping (Table 1). Tillage processes cause changes in bulk density and increase the

4

Exchangeable Sodium Percentage.

184

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

Table 2 Gully characteristics of the studied sites. Attributes

Site one

Gullying process No of gullies Percentage Area (ha) Gully density (No/ha) Mean slope (%)

Overland flow 41 65.1 31.5 2 13.3

Site two Landsliding 22 34.9

Site three

Overland flow 9 56.3 10.2 1.7 19.8

Landsliding 7 43.7

Piping 13 100 5 2.5 0.8

examination of the S–A relationship. The results presented by Vandekerckhove et al. (2000a) indicate that a negative power relationship exists between the local slope and upslope drainage area for bank gullies. However, they stated that because of different processes influencing gully initiation and sudden height drop at the rambla or barranco (arroyo) bank in one side and difficulties related to accurate measurement in the field the assumed S–A relationship could not be verified in their study area. But this does not mean that further research is not needed; on the contrary, calculation of S–A relationships in bank gullies remains essential to understanding the nature of dominant processes as well as recognizing vulnerable areas to gully erosion for future planning of land management. 4.2. Area–slope relationship

Fig. 5. Exposed soil profile showing a complex of Quaternary hydro-aeolian deposits typical of the hilly landscape in sites 1 and 2. Different layers are identified including those containing a mix of silt, clay and poor sand (Layer A) and layers of pure sandy materials (Layer B).

infiltration capacity of top soil and this, together with the differential porosity, texture, material, solubility and the rainfall regime, create different hydraulic gradients in soil horizons, which may lead to piping developing more quickly. This is particularly the case at the beginning of the growing season, when preparation of the land increases penetration of heavier rain into the soil, leading to increased dissolving of the soluble material at depth and an increase in piping risk. Piping in the area has led to crop fields becoming degraded and, although farmers attempt to surround affected land by an embankment levee to prevent runoff entering pipes and tunnels, the land is often abandoned and a shifting agriculture pattern is being seen. Therefore, to better understand bank gully initiation, more research is needed to elucidate changes and the main factors controlling piping in site three. In study site three we could not find recently developed gullies with the arroyo banks still intact in order to calculate an S–A relationship, consequently this site was not considered in the

4.2.1. Gully and catchment area–slope relationship In both sites one and two, the gullies merge with the downstream drainage network and consequently it was difficult to recognize gully networks from drainage channels. In other words, both extensions of stream channels and gully head locations in the study area are controlled by headward channel network extension rather than hillslope processes (Montgomery and Dietrich, 1989). The migration of headcuts leads to increases in the drainage density and consequently slope will increase (Leopold et al., 1964). This means that gully development in the area not only causes severe erosion and land degradation but also increases the drainage density and hillslope steepness by decreasing the length of a hillslope. Under such circumstances the entry time of runoff from hillslopes to streams declines which results in increasing the peak discharge at the outlet of the catchment (Bull and Kirkby, 2002). Therefore, the intensity of gully erosion is often an indicator of land degradation and catchment drainage evolution. The slope–area (S–A) characteristics were examined for gullies and the catchments of both site one and two. Fig. 6 shows the S–A relationship for catchments with the observed gullies S–A data overlaid, and clearly demonstrates that in all sites gully heads are distributed over the catchment, but have a minimum upslope contributing area of 50 m2 and slope of 0.0067 m/m. Previous research conducted by Willgoose et al. (1991a,1991b) and Hancock and Evans (2006) indicate that diffuse and fluvial processes together drive hillslope evolution in a catchment and these areas can be identified on an area–slope curve. Diffusive processes occupy that portion of the curve where slope generally increases with catchment area (the smaller catchment areas), while fluvial processes dominate where slope decreases with increasing catchment area. Furthermore, the study by Hancock and Evans (2006)

Table 3 Summary of measured soil parameters at each site. Site number

Soil Texture

OM%

EC dS/m

pH

Na meq/lit

K ppm

Ca meq/lit

Cl meq/lit

SAR

ESP%

Site1 Site2 site3

Sandy loam Sandy loam Sandy clay loam

0.22 0.27 0.55

2.63 3.64 40.86

7.7 7.5 6.8

9.2 12.6 566.7

28.0 24 433.3

27.8 26.2 149.5

4 18.2 154

2.3 2.4 48.2

2 2.2 40.5

185

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

used by Hancock and Evans (2006) possibly leading to larger threshold drainage area values. 4.2.2. Field measurement vs. DEM derived slope Table 4 shows the S–A relationship solved using the field and DEM derived gully slopes for the integrated (whole) dataset. As can be seen, the data derived from the field measurements has a clearer trend line (R2 = 0.49). The lower explained variance of DEM data reflecting the lower accuracy associated with derived slope values. Although the field and DEM derived relationships are different, the threshold straight lines through the lower-most points for both data sets are similar. The threshold line for data extracted from the DEM is Scr = 0.0121A0.2624, while it is Scr = 0.0146A0.2726 for measured data in the field (Fig. 8B and C). In fact both relations have relatively similar b and a coefficients. As the previous studies have illustrated, the field measurement data is more reliable and the change of coefficient often results from the methodology used (Vandaele et al., 1996; Desmet et al., 1999). The observed S–A data together with previous results taken from the literature (Montgomery and Dietrich, 1994; Vandaele et al., 1996; Vandekerckhove et al., 2000a,b) are shown in Fig. 7 (note: trend lines have been subsequently added by the author). Results in Fig. 7 show a clear variation in threshold conditions for gully incipience in different studies. This may be influenced by the type of gully, mechanisms of gully initiation and environmental characteristics controlling gully development, as well as the methodology used to assess A and S, as mentioned by other researchers (Vandaele et al., 1996; Desmet et al., 1999; Vandekerckhove et al., 2000a,b; Poesen et al., 2003).

Fig. 6. S–A relationship for site one (catchment W2, W3) and site two (bottom) with gully head data overlaid.

found that under natural vegetation cover all gullies initiated by landsliding and seepage processes were located at the convex region of the area–slope curve, with a drainage area approximately less than 1000 m2 (10 pixels). However, based on scatter plots of S–A in Fig. 6 our findings indicate that the convex (diffusive) region of the curve is at areas less than 250 m2. Two reasons may explain this difference between our results and those of Hancock and Evans (2006), including differences in: (1) the physical properties of the study area, particularly climate and geological structure and (2) the research methodology such as DEM size and algorithm used for generating slope. Our study area is located in a very susceptible geological structure with sparse vegetation cover (less than 10%) and under arid climatic conditions whereas their study area had a wet/dry tropical environment with an annual rainfall of 1389 mm and more dense vegetation cover. These conditions may therefore be responsible for the relatively low upslope area threshold for gully incision observed in our study. Furthermore, the size of the grid cells and method used for slope calculation will cause variations in results (Vandaele et al., 1996; Desmet et al., 1999), with the larger cell size

4.2.3. Coefficient and exponent comparison of the S–A relationship Comparison of the results with the other studies mentioned in Fig. 7 indicates that the results reported by Vandekerckhove et al. (2000b) are the most similar to our study. They calculated the average threshold equation (S = 0.101A0.267) for an arid area (Sierra de Gata) of the Mediterranean. In our study we found a similar exponent value (S = 0.0384A0.266) but with a different constant. Vandaele et al. (1996) reported that the exponent, b, varies between 0.26 and 0.6, but in most relations it is more or less constant and equals 0.4. Vandekerckhove et al. (2000b), however obtained the minimum value (0.104) for b. The value of the exponent b has been argued as indicating the erosion process and, possibly, dominant mechanisms of gully initiation, i.e. overland flow, landsliding and seepage whereas the influence of the methods used to assess A and S as well as characteristics of the study area are believed to be reflected in the value of constant a (Vandaele et al., 1996; Poesen et al., 2003). Compared to the exponent b the constant a ranges over several orders of magnitude, and may reflect the methods used to assess A and S as well as characteristics of the study area (Vandaele et al., 1996; Poesen et al., 2003). Montgomery and Dietrich (1994) reported that the relationship Acr / S2 is consistent across different landscape settings and gully initiation processes. According to their results the value

Table 4 Regression coefficient (intercept) a and exponent (slope) b of the fitted equation ðS ¼ aAb Þ for the DEM and field derived slope data. Coefficient of determination (R2), correlation coefficient (r) and significance level (P) of each dataset. Dataset

a

–b

R2

r

P

Integrated (1) Integrated (1a)

0.0384 0.0271

0.2661 0.3002

0.29 0.49

0.54 0.70

0.000* 0.000*

a *

Relationship is based on slope data measured at field. Significant correlation in 0.01 confidence level.

186

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

of AS2 varies between 500 and 4000. In this study we found that value for AS2 for gully heads vary between 0.9 and 425 and its average for gullies initiated by overland flow is 11 while for landsliding is approximately 48. This difference results from the importance of slope force and runoff to gully initiation, in steep areas (slope > 1 m/m) the exponent of S (2) causes an increase in the AS2 index whereas in low slope areas this exponent has a negative effect (i.e. when the slope is less than 0.5 m/m the exponent of S results in very small value for AS2). Because of sparse vegetation cover and crust on the soil surface (Poesen et al., 2002) and very intensive rainfall with a short period pattern (Mahdavi, 2004; Bull and Kirkby, 2002) Hortonian overland flow

10

Slope (m/m)

1

0.1

0.01

0.001 0.001

0.01

0.1

Area (ha)

Prosser&Abernethy,1996 Montgomery & Dietrich, Oregon Vandekerckhove et al,Spain, Gata, 2000 Integrated dataset

1

10

100

Montgomery & Dietrich ,Nevada,1994, Cheng et al.,China, 2007 Vandaele et al., Belgium,1996

Fig. 7. Relationship between upslope area and local slope of soil surface at the channel head in a variety of environments. Dashed line indicates the threshold for our integrated dataset. Dotted lines show threshold conditions for ephemeral gullies while solid lines indicate different dominant of gullying process.

in arid regions is the dominant gully initiation process and for a given area the generated runoff of a storm is more than humid regions. Gullies are predominantly initiated by extreme rainfall rather than average events (Zachar, 1982; Morgan, 2005), therefore in arid regions they can initiate more quickly than other areas. Consequently, it seems that in arid lands the threshold value for AS2 is less than other areas. 4.2.4. Factors affecting the S–A relationship (landuse, dominant process and climate) The shape of the cloud of points for the integrated dataset on Fig. 7 indicates a heterogeneous distribution among data, furthermore, the coefficient of determination (R2) (Table 4) shows only a relatively weak relationship (R2 = 0.29), although the confidence level indicates a high level of significance (P < 0.001). The data was further analysed based on the dominant initiation process. This was shown by Vandekerckhove et al. (2000a,b) and Montgomery and Dietrich (1994) to enhance the strength of the relationship and indicates a link between slope/area and the driving process of gully formation. Table 5 shows the slope–area relationship for gullies dominated by overland flow, landsliding and overland flow without farmland respectively. The relationship for overland flow dominated gullies strengthens (R2 = 0.33) relative to the integrated dataset, with the relationship for landsliding dominated gullies found to be not significant. However, the overland flow dataset included 5 outlier points which were found to be associated with agricultural lands (Fig. 8A). The effects of landuse and vegetation cover on gully initiation have been discussed in detail by several studies (e.g. Dietrich et al., 1992; Poesen et al., 1998; Vandekerckhove et al., 2000a,b). From the results of these researchers it can be drawn that relatively small changes in land cover in response to climate or landuse change, as well as soil disturbance, will dramatically decrease the upslope area required for gully initiation. These points were removed from the analysis, and the relationship was found to strengthen further (R2 = 0.52).

Fig. 8. Scattergram, fitted and threshold S–A lines separated based on dominant gully initiation mechanism including: (A) overland flow dominated gullies (DEM derived slope) – note the outlier points (shown as open circles) relate to gully heads located in agricultural lands; (B) overland flow dominated gullies with agricultural points removed (DEM derived slope), (C) overland flow dominated gullies with agricultural points removed (field derived slope) and; (D) landsliding and seepage dominated gullies.

187

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

Table 5 Regression coefficient (intercept) a and exponent (slope) b of the fitted equation ðS ¼ aAb Þ for the separated datasets. Coefficient of determination (R2), Correlation coefficient (r) and significance level (P) of the dataset separated based on dominant gully initiation mechanism. Dataset

a

–b

R2

r

P

Overland flow (2) Landsliding (3) Overland flow without farmland (4)

0.0369 0.0941 0.0438

0.1998 0.1278 0.1823

0.33 0.026 0.52

0.6 0.16 0.72

0.000* 0.252ns 0.000*

ns: Non significant. * Significant correlation in 0.01 confidence level.

Table 6 S–A threshold relationships reported in the literature. To compare the results, the relation is transformed in the far right column to S = aAb in which the area is expressed in hectares. Reference

Region

Vegetation

Geology

Annual rainfall (mm)

Relation

A

S

Montgomery and Dietrich, 1988

Open wood land

Granite

260

S = 0.35A0.6

ha

Montgomery and Dietrich, 1992

Sierra Nevada –

–

–

–

AS2 = 500–4000

m2

Vandaele et al., 1996

Belgium

Agriculture

Loess belt

700

S = 0.025A0.4

ha

Vandekerckhove et al. 2000b

Rangeland

Alluvial Fan Litic and Eutric soil Loess

182

S = 0.101A0.267

ha

900

S = 0.058A0.3

ha

Vandekerckhove et al. (2000a,b)

Sierra de Gata Loess Plateau Lesvos

405

S = 0.218A0.211

ha

Willgoose et al., (cited by Hancock and Evans, 2007)a

Northern Australia

Forest

Volcanic Rock and Eutric Leptisol ERA (Uranium)

1389

2.5A0.4S0.3 P 25

m2

m/ m m/ m m/ m m/ m m/ m m/ m m/ m

Cheng et al., 2007

Changed cropland to grassland Rangeland

Parameter units

Transformed relation

S = 0.35A0.6 S = 0.05A0.5 S = 0.025A0.4 S = 0.101A0.267 S = 0.058A0.3 S = 0.218A0.211 S P 0.01A1.33

a The source relation introduced by Willgoos et al. (1991a) is 2.5A0.4S0.3 P at and Hancock and Evans (2007) used it for evaluation gully position by taking into account the value of 25 as at.

Montgomery and Dietrich (1992) and Dietrich et al. (1992) have divided the landscape into overland flow, seepage and saturation erosion as well as landsliding based on the dominant process responsible for channel initiation. They have reported that the relation between A and critical S in all cases is a linear inverse one, but it is not true for saturated soil and seepage erosion often associated with landsliding (Vandekerckhove et al., 2000a,b). This is consistent with our findings, with no statistically significant S–A relationship found for gullies dominated by landsliding. This is exacerbated in our study by the lithology, and particularly the presence of a variable cohesionless-sand layer (Fig. 5). Table 6 shows S–A relationships and environmental characteristics of studies undertaken by other researchers. Comparison of Tables 5 and 6 illustrate that our findings are different from most previous research, except Vandekerckhove et al. (2000a,b). This is most likely attributable to the environmental characteristics of the study area, as well as methodologies used because none of the cited studies are located in arid conditions. Montgomery and Dietrich (1988) stated that in dry regions the thresholds line plots above wet regions and concluded that for a given local slope, the source area required to initiate a channel head should increase with increasing aridity to produce the same critical combination of runoff and local slope at the channel head (Montgomery and Dietrich, 1988). In contrast, Vandekerckhove et al. (2000a,b) indicated that this does not apply in their study area and consequently stated that the relative position of the threshold line is controlled by other factors beyond climatic condition. In our study, we also found a relatively low threshold line. In arid regions, this may indicate that sparse vegetation cover, and the intense, albeit sporadic, nature of the rainfall can generate sufficient runoff from small catchments for gully initiation.

The value of the exponent of b for the integrated dataset (0.266) is similar to that of Sierra de Gata (0.267) and Lesvos (0.211). Both studies are located in rangelands, although the Lesvos study area is considerably wetter. Although the study area of Montgomery and Dietrich (1988) is also located in an arid region (Sierra Nevada) the b is higher. This is likely to be because the area is more heavily vegetated (open oak woodland and grassland) and has a less erosive geology (old granite) than our study area. However, similar to their results we found that overland flow is the dominant process in gentle gradients and landsliding is responsible for channel initiation in steep regions. These researchers illustrated that overland flow and landsliding gullies can be separated based on slope (i.e. S > 50% for landsliding), but this was not possible in our study due to the combined effect of the processes in many gullies. 4.2.5. Application of threshold to predict areas vulnerable to gully initiation Identifying the potential for gully initiation in the study area can help land managers to better understand the potential damage of gullying and inform the development of soil conservation projects. The S–A relationship can be applied to distinguish the areas prone to gully erosion (Vandaele et al., 1996; Desmet et al., 1999), given similar geomorphology, climate and soil conditions. A small catchment (16.3 ha) located approximately 1800 m from the existing sites was identified for model application. All gully head locations were mapped (12 identified) and a DEM was used to extract the upslope drainage area and slope for each grid cell. Most previous studies (Patton and Schumm, 1975; Vandaele et al., 1996; Cheng et al., 2007) used a straight fitted line through the lower-most points of data. This line is the conservative

188

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

Fig. 9. Location of the observed gullies and predicted areas prone to gullying in a comparable catchment near the studied sites: (A) and (B) are derived from relation 4 and 2 (Table 5) respectively; (C) is based on a straight line fitted through the lowest points (from Fig. 8B) and; (D) is the area predicted using equation proposed by Hancock and Evans (2006) for at = 25.

boundary of gully and non-gully areas. In this study the average critical slope gradient was calculated through a threshold relationship obtained from the previous stage (Eqs. 2 and 4 in Table 5) and threshold line of Fig. 8B. In addition we used another equation obtained by Hancock and Evans (2006). Fig. 9A and B indicate that the solved S–A threshold relationship are able to capture 8 out of 12 gully positions observed in the field (68%). The map also illustrates that areas prone to gullying are distributed in all parts of the catchment, but upper areas are more susceptible than the other areas. Although a good consistency between observed gullies and predicted areas has been shown, the relative area prone to gullying in Fig. 9A and B are different. In other words, when we applied relation 2 (Table 5) about 6.5% of the catchment area was found to be prone to gullying while for relation 4 it increases to 12%. This is likely to indicate that areas subject to agriculture are more at risk of gully erosion. This reveals that by changing the landuse from rangeland to agriculture the area prone to gullying increases by up to two times. This finding is consistent with previous results reported by Montgomery and Dietrich (1994) and Vandekerckhove et al. (2000a,b). It is, therefore important that the future land management should ensure that areas vulnerable to gullying are not disturbed. Future research should also further examine this relationship under different landuses in the region to provide more information for effective management. The map resulting from solving the threshold line on the two lowest points Fig. 9C captures all of the observed gullies, but identifies 48% of catchment as prone to gullying and is therefore very conservative. Fig. 9D indicates that areas prone to gullying occur throughout the stream channel network, similar to the results reported by Hancock and Evans (2006). Although the area prone to gullying ob-

tained from their relation (5.4%) is near to results obtained from relation 2, the pattern is very different and only 42% of total observed gullies were identified. This is likely to be because while most research has concentrated on gully initiation points, the Willgoose et al. (1991a,b) study focuses on predicting the trajectories of channel initiation and hillslope evolution based on the Hortonian threshold concept (Bull and Kirkby, 2002). This relation appears to better predict gullying through channel network extension, however future research work is needed to calibrate this relationship to our study area. 5. Conclusion Throughout the arid regions of Iran gully erosion is a widespread process responsible for land degradation and has important implications for maintaining local infrastructure. Despite the significance of sediment contribution from gully erosion in arid and semi arid regions of Iran, no study has been conducted to understand the threshold conditions of gully initiation. Indeed, studies in arid regions throughout the world are relatively limited. This study has shown that the dominant process for incipient gullying is controlled by geomorphic characteristics and topography of the landscape. Generally, as the slope increases, the effect of landsliding on gully initiation becomes more significant. All bank gullies were found to be associated with highly sodic soil (EC = 40.86 dS/m and SAR = 48.2). In addition tillage processes increased the infiltration of top soil through changing bulk density and was found to develop piping more quickly. However, under natural conditions it is very difficult to identify headcut initiated by only one mechanism and normally a complex of mechanisms are responsible for gully initiation.

A. Nazari Samani et al. / Journal of Asian Earth Sciences 35 (2009) 180–189

Gullies were found to have formed throughout the catchment, with a large range in drainage area and slope. Grouping data describing upslope drainage area and local slope based on dominant gullying process indicated a significant relationship between A and S for gullies initiated by overland flow, but not for landsliding. Data discussed in this paper clearly indicated that landuse changes, that increase the effective rainfall intensity and runoff properties, decrease the topographic threshold and accelerate channel network expansion. Comparison of the threshold line found for overland flow gullies in this study with other studies suggests that it is relatively low. However, this is consistent with other studies in arid regions and indicates that arid regions may be more susceptible to gullying due to the low vegetation cover and high variation of rainfall. This study has also shown that the threshold S–A relationship can be used to predict vulnerable areas to gullying with an acceptable accuracy in order to examine the potential impact of different land cover. However, due to the complexity of gully erosion mechanisms and the limited amount of research undertaken in arid regions more research effort is needed to understand gullying processes and gully hazard zonation. Acknowledgments The research for this paper was funded by the Iran National Science Foundation (INSF) and this support is gratefully acknowledged. The authors would like to thank anonymous reviewers and Dr. Tim Kusky, associate editor of the Journal of Asian Earth Sciences for their helpful suggestions and comments. Special thanks are due to Professor Robert J. Wasson, Charles Darwin University, for his fruitful comments on an earlier draft and thorough review. The aerial photographs were provided by the National Cartographic Center of Iran. References Adelpour, A.A., Soufi, M., Behnia, A.K., 2004. Channel erosion thresholds for different land uses assessed by concentrated overland flow on a silty loam. In: Conserving Soil and Water for Society: Sharing Solutions. ISCO 200-13th International Soil Conservation Organisation Conference, Brisbane, Australia. Ahmadi, H., 2004. Introduction to Iranian Model of Desertification Assessment (IMDPA), Unpubl. Technical Report for presenting in Committee of Science and Technology, United Nation Convention to combat desertification, The 7th Session of Conference of Parties (COP7), Nairobi, Kenya, 56 pp. Begin, Z.B., Schumm, S.A., 1979. Instability of alluvial valley floors: a method for its assessment. Trans. Am. Soc. Agric. Eng. 22, 347–350. Bull, L.J., Kirkby, M.J., 2002. Channel heads and channel extension. In: Bull, L.J., Kirkby, M.J. (Eds.), Dryland Rivers: Hydrology and Geomorphology of Semi-arid Channels. Wiley, Chichester, UK, pp. 265–298. Cheng, H., Zou, X., Wu, Y., Zhang, C., Zheng, Q., Jiang, Zh., 2007. Morphology parameters of ephemeral gully in characteristics hillslopes on the Loess Plateau of China. Soil & Tillage Research 94, 4–14. Desmet, P.J.J., Poesen, J., Govers, G., Vandaele, K., 1999. Importance of slope gradient and contributing area for optimal prediction of the initiation and trajectory of ephemeral gullies. Catena 37, 377–392. Dietrich, W.E., Wilson, C.J., Montgomery, D.R., McKean, J., Bauer, R., 1992. Erosion thresholds and land surface morphology. Geology 20, 675–679. Dietrich, W.E., Dunne, T., 1993. The channel head. In: Beven, K., Kirkby, M.J. (Eds.), Channel Network Hydrology. Wiley, Chichester, pp. 175–219. Faulkner, H., Spivey, D., Alexander, R., 2000. The role of some site geochemical processes in the development and stabilisation of three badland sites in Almería, Southern Spain. Geomorphology 35, 87–99. Forest and Rangeland Organization, 2004. National Activities Program, Technical Report for United Nation Convention of Combat Desertification. I. R. Iran, 75 pp. Ghoddousi, J., 2002. Gully erosion morphology modelling and hazard zonation (Study area: Zanjanrood drainage basin), Unpubl. Ph.D thesis, Department of Arid and Mountainous Region Reclamation, University of Tehran, 368 pp, (abstract in English). Grise, S., Arctur, D., Booth, B., 2002. Implementing Arc Hydro. In: Maidment, D.R. (Ed.), Arc Hydro, GIS for Water Resources. ESRI Press, pp. 176–201. Gutierrez, M., Benito, G., Rodriguez, J., 1988. Piping in badland areas of the Middle Ebro Basin, Spain, In: Harvey, A.M., Sala, M. (Eds.), Geomorphic Processes in

189

Environments with Strong Seasonal Contrasts. Geomorphic Systems, Catena Supplement 13 II, pp. 49–60. Hakimkhani, Sh., 2006. An Investigation on using tracers in fluvial fine sediment sources fingerprinting (case study: the basin of Pouldasht flood spreading system, Makoo township), Unpubl. PhD thesis, Department of Arid and Mountainous Region Reclamation, University of Tehran, 255 pp, (abstract in English). Hancock, G.R., Evans, K.G., 2006. Gully position, characteristics and geomorphic thresholds in an undisturbed catchment in northern Australia. Hydrological Processes 20, 2935–2951. Harvey, A., 1982. The role of piping in the development of badlands and gully systems in Southeast Spain. In: Bryan, R., Yair, A. (Eds.), Badland Geomorphology and Piping. Geo Books Geo Abstracts, Norwich, pp. 317–355. Leopold, L.B., Wolman, M.G., Miller, J.P., 1964. Fluvial Processes in Geomorphology. W.H. Freeman, San Francisco. 525 pp. Mahdavi, M., 2004. Applied Hydrology, vol. 2, fifth ed. University of Tehran Press, Iran (In Persian), 420 pp. Martın-Penela, A.J., 1994. Pipe and gully systems development in the Almanzora basin (Southeast Spain). Zeitschrift fur Geomorphology 38, 207–222. Montgomery, D.R., Dietrich, W.E., 1988. Where do channels begin? Nature 336, 232–234. Montgomery, D.R., Dietrich, W.E., 1989. Source areas, drainage density, and channel initiation. Water Resources Research 25 (8), 1907–1918. Montgomery, D.R., Dietrich, W.E., 1992. Channels initiation and the problem of landscape scale. Science 255, 826–830. Montgomery, D.R., Dietrich, W.E., 1994. Landscape dissection and drainage area– slope thresholds. In: Kirkby, M.J. (Ed.), Process Models and Theoretical Geomorphology. Wiley, Chichester, UK, pp. 221–246. Morgan, R.P.C., 2005. Soil Erosion and Conservation, third ed. Blackwell, Malden, USA. 299 pp. Morgan, R.P.C., Mngomezulu, D., 2003. Threshold conditions for initiation of valleyside gullies in the Middle Veld of Swaziland. Catena 50, 401–414. Nordstrom, K.F., 1988. Dune grading along the Oregon coast, USA: a changing environmental policy. Applied Geography 8, 101–116. Oostwoud Wijdenes, D.J., Bryan, R.B., 1994. The significance of gully headcuts as a source of sediment on low-angle slopes at Baringo, Kenya, and initial control measures. Adv. Geoecol. 27, 205–231. Oostwoud Wijdenes, D.J., Poesen, J., Vandekerckhove, L., Ghesquiere, M., 2000. Spatial distribution of gully head activity and sediment supply along an ephemeral channel in a Mediterranean environment. Catena 39, 147–167. Patton, P.C., Schumm, S.A., 1975. Gully erosion, Northwestern Colorado: a threshold phenomenon. Geology 3, 88–90. Poesen, J., Vandaele, K., van Wesemael, B., 1998. Gully erosion: importance and model implications. In: Boardman J., Favis-Mortlock D. (Eds.), Modelling Soil Erosion by Water. NATO ASI Series, vol I 55, Spinger-Verlag, Berlin. Poesen, J., Vandekerckhove, L., Nachtergaele, J., Oostwoud Wijdenes, D., Verstraeten, G., van Wesemael, B., 2002. Gully erosion in dryland environments. In: Bull, L.J., Kirkby, M.J. (Eds.), Dryland Rivers: Hydrology and Geomorphology of Semi-Arid Channels. Wiley, Chichester, UK, pp. 229–262. Poesen, J., Nachtergale, J., Vertstraeten, G., Valentin, C., 2003. Gully erosion and environmental change: importance and research needs. Catena 50 (2–4), 91– 134. Prosser, I.P., Abernethy, B., 1996. Predicting the topographic limits to a gully network using a digital terrain model and process thresholds. Water Resources Research 32 (7), 2289–2298. Romero Diaz, A., Marin Sanleandro, P., Sanchez Soriano, A., Belmonte Serrato, F., Faulkner, H., 2007. The causes of piping in a set of abandoned agricultural terraces in southeast Spain. Catena 69, 282–293. Sidorchuk, A., 1999. Dynamic and static models of gully erosion. Catena 37, 401– 414. Valentin, C., Yong Poesen, J., Li., 2005. Gully erosion: impacts, factors and control. Catena 63, 132–153. Vandaele, K., Poesen, J., Govers, G., van Wesemael, B., 1996. Geomorphic threshold conditions for ephemeral gully incision. Geomorphology 16 (2), 161–173. Vandekerckhove, L., Poesen, J., Oostwoud Wijdenes, D., Gyssels, G., Beuselinck, L., De Luna, E., 2000a. Characteristics and controlling factors of bank gullies in two semi-arid Mediterranean environments. Geomorphology 33, 37–58. Vandekerckhove, L., Poesen, J., Oostwoud Wijdenes, D., Nachtergaele, J., Kosmas, C., Roxo, M.J., Figueiredo, T.DE., 2000b. Thresholds for gully initiation and sedimentation in Mediterranean Europe. Earth Surface Processes and Landforms 25, 1201–1220. Wasson, R.J., Caitcheon, G., Murray, A.S., McCulloch, M., Quade, J., 2002. Sourcing sediment using multiple tracers in the catchment of Lake Argyle, Nortwestern Australia. Environmental Management 29 (5), 634–646. Willgoose, G.R., Bras, R.L., Rodriguez-Iturbe, I., 1991a. A physical explanation of an observed link area–slope relationship. Water Resources Research 27 (7), 1697– 1702. Willgoose, G.R., Bras, R.L., Rodriguez-Iturbe, I., 1991b. Results from a new model of river basin evolution. Earth Surface Processes and Landforms 16, 237–254. Zachar, D., 1982. Soil Erosion. Elsevier Science Publication, North Holland. 522 pp. Zevenbergen, L.W., Thorne, C.R., 1987. Quantitative analysis of land surface topography. Earth Surface processes and Landforms 12, 47–56.