Reconstructing long-term gully dynamics in Mediterranean agricultural areas

Gully erosion is an important erosive process, especially in Mediterranean basins. However, the long-term dynamics of gully networks and the variation of sediment production in gullies is not well known. Available studies are often done over a few years only, while many gully networks form, grow, and change in response to environmental and land use or management changes over a long period. In order to clarify the effect of these changes, it is important to analyze the evolution of the gully network with a high temporal resolution. This study aims at analyzing gully morphodynamics over a long time scale 5 (1956-2013) in a large Mediterranean area in order to quantify gully erosion processes and its contribution to overall sediment dynamics. A gully network of 20 km2 located in SW Spain, has been analyzed using a sequence of 10 aerial photographs in the period 1956-2013. The extension of the gully network both increased and decreased in the study period. Gully drainage density varied between 1.93 km km-2 in 1956, with a minimum of 1.37 km km-2 in 1980 and a maximum of 5.40 km km-2 in 2013. The main 10 controlling factor of gully activity appeared to be rainfall, while land use changes were found to have only an indirect effect. A new Monte Carlo-based approach was proposed to reconstruct gully erosion rates from orthophotos. Gully erosion rates were found to be relatively stable between 1956-2009, with a mean value of 11.2 ton ha-1yr-1, while in the period 2009-2011, characterized by extreme winter rainfalls, this value increased significantly, to 591 ton ha-1yr-1. These results show that gully erosion rates are highly variable and that a simple interpolation between the start and end date would highly underestimate 15 gully contribution during certain years, such as for example between 2009-2011. This illustrates the importance of the applied methodology using a high temporal resolution of orthophotos.

Gentle hills prevail in the study area except from the south and the centre east where steeper ones exist (up to 32%). Altitudes range from 233 to 558 m high and mean slopes are 13%. The soils in the area are dominated by Vertisols, formed mainly in marls and calcareous sandstones deposited during the Miopliocene.
Currently the dominating land uses are olive orchards and herbaceous crops covering almost the whole area, except some 5% of the surface area occupied by grassland. Mean annual precipitation varies between 500 and 600 mm (Córdoba Airport station 5 and Baena RIA station). The distribution of the precipitation shows a marked dry season between June and September while the main wet period occurs from October to May.

Rainfall characterization
Characterization of the rainfall regime was performed from daily rainfall collected in the period 1956-2013 at Castro del Río meteorological station (37.69 • N, 4.47 • W), belonging to the Spanish National Meteorological Agency (AEMET). Iso-10 lated data gaps between 1970 and 1971 were completed from the data recorded at Cañete de las Torres meteorological station (37.83 • N, 4.36 • W, Phytosanitary Warnings Network of Andalusia, RAIF) and Córdoba Airport meteorological station (37.84 • N 4.84 • W, AEMET). Anomalies in annual rainfall were evaluated by means of normalization, through average and standard deviation of annual rainfall for a 57 years period , following Martínez-Casasnovas et al. (2003). Values falling outside the interval R mean (average rainfall) ± sd (standard deviation), which correspond to the normalized values >1 15 and <-1, were considered anomalies.
The frequency distribution of daily rainfall above a threshold value of 13 mm was analysed, considering this as the minimum rainfall that produces erosive effects as proposed by Wischmeier and Smith (1978) and Renard et al. (1997). In addition, frequency distribution of records above the average daily rainfall event plus the standard deviation were also analysed, assuming that these events represent the extreme rainfall events within the study period.  Table 1. The working scale in the photointerpretation processes was established to 1:5000 for the whole dataset.

Land use 25
Land use in the study area for 2001, 2005, 2009, 2011, and 2013 was derived from the respective orthophotos while for the rest of the years (1956, 1980, 1984, 1999, 2003, and 2007) existing Maps of the Land Use and Vegetation Cover of Andalusia (Red de Información Ambiental de Andalucía, REDIAM) were used. Different land uses present in the area were simplified to three classes as shown in Table 2.

Gully network length
Gully length was obtained by digitizing the extension of the gully network for each available year, distinguishing between newly incised and infilling stretches. Gully network was decomposed in m y segments, where subscript y indicates the year.
Each segment comprises the length between consecutive junctions (Fig. 2). Due to changes in the drainage network during the study period, the number of segments ranged between 108 in 1980 and 940 in 2013. The total length of the drainage network 5 for a given year, L y , was calculated as the sum of the lengths of individual segments, l y,i with m y equal to the total number of individual segments of the gully network for each digitalized year.

Gully network width
In order to measure gully width in a representative way, 35 stretches were selected from the earliest digitalized gully network 10 of 1956, covering a wide range of widths. Gully width was measured at the same locations on later orthophotos, allowing the evaluation of the widening process during the complete study period.

Field campaign
During 2013 and 2014 several field campaigns were conducted to measured current gully widths and depths with measuring tape and a clinometer (Suunto PM-5/360 PC). Gully top width and depth were measured at 27 representative sections distributed 15 randomly over the gullies catchments. These representative sections covered the entire range of width and depth variability, including different landscape positions, from upstream close to the divide to the junction with the stream network, and both in gullies on herbaceous crops and under olive orchards gullies.

Monte Carlo-based simulations
Although gully length for the different years between 1956 and 2013 could be determined directly from observations using the 20 available air photographs, determination of the gully volume was not straightforward. As we used freely available orthophotos, it was only possible to measure the size of the gullies in two dimensions and no measure of depth was readily available.
Also observations of gully width for each year were limited to the representative sections measured on the orthophotos of that particular year and therefore included a term of uncertainty as the real population mean remained unknown.
Estimation of overall gully network volume for each year,V y , was therefore tackled by conducting a Monte Carlo simulation 25 in which a volume and an associated uncertainty were calculated for every single gully segment, l y,i , described in paragraph 2.3.2 (Fig. 2). For each year, y, a set of n = 1000 estimated cross area sections, S y,i = {s y,i,j , j = 1, ..., n} for every single segment, l y,i , were generated as show in Figure 3, which required the generation of sets of width and depth values for each year. Each generated section is calculated as where k is a shape factor, and w y,i,j , and d y,i,j , the simulated gully width and depth respectively. Field observations suggested 5 that a triangular section is a reasonable approximation of most gully sections, so a shape factor k = 0.5 was adopted in order to compute the simulated sections.
To generate a representative measure of gully width, first of all, the gully width distribution measured for each year by photointerpretation at the representative sections was fitted to different probability distribution functions (normal or Gaussian, gamma, lognormal, exponential and Weibull) using the maximum likelihood method. Next, goodness of fit was evaluated for these 10 different distributions by means of the Kolmogorov-Smirnov statistics. Finally, the best overall fitting theoretical probability distribution was selected to obtain the necessary parameters (µ y , σ y ) to generate n random simulations of representative gully widths for any particular year.
The estimation of gully depth for each year was based on the field data gathered in 2013-14. In order to estimate depth for previous years, firstly a width-depth relationship was estimated by linear regression analysis from the collected field data. Such 15 a relationship could only be established for the present-day situation. Uncertainty on this linear width-depth relation was then taken into account by computing the estimated intercept, slope and their respective standard deviations (a, b, s a , s b ). Assuming a normal distribution, a set of one thousand slope and intercept pairs were simulated. Depths for unique segments (D y,i ) were then derivate from simulated widths and slope-intercept pairs.
Finally, a set of n simulated volumes V y,i = {v y,i,j , j = 1, ..., n} was calculated for each year and segment multiplying indi-20 vidual measured lengths by the simulated sections ( Fig. 3) v y,i,j = s y,i,j l y,i A set of n different simulated volumes of the complete gully network for a particular year V y was eventually calculated as the sum of volumes of single segments v y,i,j Finally average volume of the total gully network for a given year,V y , was computed as Erosion rates were then obtained from the different between pairs of simulated volumes in consecutive dates divided by the duration of the period.

Rainfall characteristics during the study period
The annual rainfall depths in the analysed period ranged between 180 mm In 1984-1999 eight out of fifteen records were greater than the 0.75 percentile, and 6 of them were considered anomalies since 10 they were greater than the average annual rainfall plus the standard deviation. In the period 2009-2011, both years recorded annual rainfall amounts higher than the standard deviation and can thus be considered an anomalous extreme rainy period.  Land use experienced a progressive conversion from herbaceous crops to olive orchards as shown in Figure 6. In the study period, olive orchards grew from 13% to 63% of occupation of the land use in the study area. At the same time herbaceous crops decreased from 85% to 35% of the occupied land. The main land use change occurred between 1984 and 1999, when the olive orchards passed from occupying 25% to 48% of the total area. The highest rates of change however were observed in the period 2005-2007 with more than 4% rate of annual land use change from herbaceous crop to olive orchards.  analysis on the length and area ratio showed that the drainage density has grown from 17.2 m ha -1 to 53.3 m ha -1 . In most of the analyzed period variations on drainage density occurs are small. However, there are two significant periods where the increase is very high and that account for the main increases in the overall value. From 1984 to 1999 and 2009 to 2011 there was an increment of 14.6 m ha -1 and 23.6 m ha -1 respectively, which account for 84% of the total drainage density growth. When comparing these gully length dynamics to controlling factors of land use and rainfall, it can be seen in table 3 that these rapid 5 growth could be related to extreme rainfall events that occurred in 1997 and anomalous rainy periods in 2009-2011. In contrast, in some periods as for instance in 1956-1980, 1999-2001, 2001-2005 and 2007-2009 the gully network experimented several decreases in the drainage density, although in no case this decrease was more than 4 m ha -1 , and can therefore be considered modest. These decreases can be directly related to farming operations where farmers fill in the upstream gully stretches that are limited in depth and can be considered ephemeral gullies.
10 Figure 8 shows the frequency distribution of headcut growth and infilling of individual gullies for the different periods between 1956-2013. Some of the observation periods show a balance between infilling and growing reaches, which leads to a very minor overall change of the total gully network length. During a few distinct intervals however, 1984-1999 and 2009-2011, this balance shifts drastically and results in a fast increase of the gully network's total length, as can be seen in Figure 9. This can be in part explained because in these two periods infillings are almost negligible ( Fig. 8 and Fig. 9). However, in Figure   15 9 growth of the gully at the end of those periods (1999 and 2011) is much higher (31 km and 49 km) than those from the other end periods (13 km as the highest value), which clearly shows that gully growth is the dominant process controlling gully dynamics in those periods. during some other periods, as for instance in 1956-1980, 1999-2001, 2001-2005 and 2007-2009, the balance between infilling and growing stretches resulted in a net reduction of the total gully network length. Infilling gully stretches identified during photointerpretation, may be classified in two different types: those made while regular tilling operations at the end of the summer, usually in the order of several tens of meters and those resulting from land levelling during phases of land use change, which may reach some hundreds of meters.

25
Extraordinary annual rainfalls as well as individual extreme precipitation events seem to be the main factors that can be linked to gully retreat (Table 3). Land use does not seem to control these observed peaks in gully length increase directly. However, land use change could have contributed to the rainfall extremes inducing high peak discharges, because since 1956 a shift from cereal crops to olive orchards occurred in half of the study area, and which was especially intensive from 1984 forward. Young olive trees with limited root systems and small canopies leave an important bare soil surface and little protection to overland 30 flow or gully headcut advance. However, further analysis should be done in order to confirm this hypothesis.

Gully network width dynamics
Top width at the representative cross sections, as derived from the orthophotos dataset, experienced a continued widening over time (Fig. 10). While at the beginning of the study period (1956)

Width and Depth relationship
In order to compute the volume of the gully network, depths at the different stretches were derived from the Monte Carlo simulated widths using a width-depth relation derived from field work, shown in Figure 11. A coefficient of determination R 2 = 0.83 was obtained from a logarithm-based fitting, with slope, intercept and their standard deviation respectively 1.73 ± 0.16 and 0.55 ± 0.32. Normal deviates based on those coefficients were used to generate 1000 width and depth pairs.

Gully erosion rate dynamics
Dynamics of gully erosion rate are shown in Figure 13. Maximum erosion rate was reached in the period 2009-2011 when 591 ton ha -1 yr -1 were lost according to the simulation process. Minimum erosion rate (-5.21 ton ha -1 yr -1 ) was registered in 5 the period [2005][2006][2007]. Negative values here reflect the decrease of the gully network volume, and therefore should not be considered an erosion rate but an infilling rate. Average erosion rate for the whole study period was 39.7 ton ha -1 yr -1 .

Discussion
The average gully erosion rate of 39.7 ton ha -1 yr -1 obtained in this study, by means of photo-interpretation techniques combined with stochastic methods, are of the same order of magnitude with those found in literature in Mediterranean basins. Oostwoud Most importantly, the results show a wide variability in gully erosion rates, ranging between -5.21 and 591 ton ha -1 yr -1 . This includes periods dominated by infilling and rapid growth, underlining the importance of measuring erosion rates at the finest 20 temporal resolution possible in order to overcome under-and/or overestimations in sediment production. Such variability is in part explained by the inherent irregularity of the local rainfall regime which appears to be the main controlling factor for gully erosion at this site. However, land use change has played an important role intensifying in some cases and masking in other cases the rates of gully erosion. For instance, in the initial period between 1956 and 1980, erosion rate shows a negative value.
However, given the length of this period and since there are some particular years (i.e. 1961-1962) with extreme rainfalls it is 25 likely that positive gully growth occurred during this period that was later masked by infilling. Moreover, the data presented here clearly show that in Mediterranean areas the gully growth dynamics are different than in temperate areas. A review of different studies on gully growth over time by Poesen et al. (2006) indicated a rapid initial growth, followed by a stable phase with slow growth for "mature" gullies. Data for this study was from the temperate loess belt or from lab experiments under constant discharge conditions. In our case, with a high variability in natural rainfall, even after several decades, intense growth phases were observed. As stated before, these could mainly be attributed to an increase of the gully's cross sections, and less to a gully headcut retreat. Therefore, models such as CHILD or REGEM, which have been applied with success to gully modelling, but focus mainly on headcut activities, would probably not yield good results in this case.
From a wider geomorphological perspective, other phenomena such as lowering of the base level and incision of the river bed could be suggested as a cause of the progressive increase on the erosion rate. Nevertheless, there was no field evidence of this.

5
In addition, the Guadalquivir river basin is a highly regulated river, with many dams, which could be expected to limit any such effects.
Gully erosion rates computed between the start and the end of the study period would incur in gross underestimation. Erosion  (2007) found that on average the three largest daily events per year accounted for more than 50% of the total sediment exported from the basin. The so-called time compression of Mediterranean climate with respect to soil erosion is therefore very high, as is demonstrated by the data from this study.

15
The Monte Carlo stochastic modelling performed also allow to identify that while gully length dynamics ( Fig. 9) could explain some of the rapid increases in the volume and erosion rate computed, widening processes ( Fig. 10) determine the shape of volume curve (Fig. 12) pointing out the importance of that parameter in the computed volume as opposite, in this particular case, as suggested by other authors who found for other areas and climates that the leading controlling parameter is gully length (Nachtergaele and Poesen, 1999). This observation will lead future field work and modelling efforts, which should not only 20 consider gully headcut advance, but also on the mechanisms of gully sidewall collapse and erosion. Possibly a very important factor here, also in order to control the gully growth, is the effect that roots may have on stabilizing the gully walls (De Baets et al., 2008).
The main advantage of the new method described here, is that by means of Monte Carlo simulation, an estimation of the uncertainty associate along with the measure of gully erosion volume is generated. This is especially relevant when suitable 25 knowledge of erosion dynamics is required and management systems need to be evaluated or compared.

Conclusions
A new method was presented to evaluate gully growth over decadal time scales, combining airphotos interpretation with a stochastic approach through Monte Carlo modelling for the channel section parameters. This resulted to be a reliable procedure to determine gully network dynamics over time. Uncertainty ranges obtained in the simulation provide an unprecedented rates were in accordance with previous studies in Mediterranean basins. The fluctuations in erosion rates were mainly attributed to the variability in rainfall regime variations, likely exacerbated by land use changes, although further research -using physical modelling of runoff, gully headcut retreat rates and sidewall dynamics-should be made at this last point.
Simple interpolation between the start and end date would highly underestimate gully contribution during certain years, as it could be verified when comparing average erosion rate (39.7 ton ha -1 yr -1 ) with punctual erosion rates at the end of the study    1956 1980 1984 1999 2001 2005 2007 2009 Table 3. Land use, rainfall indicators and gully growth. f h and fo: fractions of surface dedicated to herbaceous and olive crops, in the first year of each period. nle: number of 24 hours rainfall events per year higher than 13 mm, nleo: number of 24 hours rainfall events per year over the average 24 hours rainfall plus the standard deviation, Rmax: highest daily rain depth registered within the period, M AR: Mean annual rainfall in the period. ∆L: total, and ∆L/∆t, partial increase in gully length, and GH: gully headcut growth, averaged over the area.