ANALYSIS OF THE EVOLUTION OF GULLY EROSION IN OLIVE GROVES USING PHOTOGRAMMETRY TECHNIQUES. RELATIONSHIPS WITH RAINFALL REGIME

Gully erosion is one of the main processes of soil degradation, reaching 50%–90% of total erosion in basins. As erosion processes are related to rainfall regime, the depletion and deposition rates can be increasing in a climate change scenario. This paper deals with the quantification of erosion processes in an active gully affecting olive groves of the province of Jaén (southern Spain), using geomatics techniques (photogrammetry and LiDAR). Eight historical aerial flights from 1980 to 2016, a LiDAR dataset (2014) and 2 recent UAS surveys (2019-2020) were used and processed in a common reference system with the support of field GNSS ground control points. Then, DSMs and orthophotographs were obtained and DSMs of differences (DoDs) calculated, from which we can identify gullies, calculate the depletion and deposition areas, and estimate height differences and volumes involved. These analyses result in an average depletion of about -1.6 m (incision) and a waste volume up to 30000 m3 (soil losses), which lead to a rates of 0.05 m/year and -44 t/ha*year, respectively. These rates are very different along the considered periods, reaching the maximum values (near -300-450 t/ha*year) in 2009-2011 and 2011–2013, coinciding with the periods of higher rainfall in the last fifty years, that probably have underwent an increase of 10-30% in the last decades. Thus, the evolution of the gully area for 40 years has been analysed in relation to the rainfall regime that has been established from the daily rainfall data. * Corresponding author


INTRODUCTION
Soil erosion is one of the main phenomena of environmental degradation (Borrelli et al., 2013) that could increase significantly in the coming decades in relation with processes associated to climate change (Yang et al., 2003;Li and Fang, 2016). Many studies have been developed at different spatial and temporal scales. In the plot scale, empirical measurements or different approaches based on the USLE model (Gómez et al., 2003) have been used. Meanwhile, some authors (Poesen et al., 2003) have estimated that gully erosion processes could explain between 50% and 90% of total erosion at basin scales.
Photointerpretation and image analysis allows the estimation of lengths, widths and densities of the gully systems. In addition, photogrammetric methods allow the generation of digital elevation models (DEMs) for the estimation of gully depths and volumes, and orthoimages. The availability of data at different epochs leads to multitemporal analysis of gully areas, either in 2D (Hayas et al., 2017) or 3D approaches, usually by means of the calculation of DEMs of Differences (DoDs). These methods are often combined with point capture techniques such as LiDAR, GNSS or TS (Brasington et al, 2003;Castillo et al., 2012;Koci et al., 2017;Fernández et al., 2020a). UAS in their different modalities are well suited for very high resolution and accurate surveys in areas between 0.01-100 km 2 . They are very used in intermediate scales between terrestrial and aerial techniques, keeping low costs and allowing repetitive studies with high temporal frequency. Most of the current studies using multicopter UAS are of centimeter resolution (Eltner et al., 2013;Fernández et al., 2020b). After the image alignment, high resolution and precision DEMs are obtained with Structure from Motion (SfM) and Multi Video Stereo (MVS) techniques, from which morphometric measurements and volumetric estimations can be addressed.
The accuracy of data and models can be estimated from several methods: RMS errors of the image orientation at ground control and/or check points; the comparison of the heights with point samples measured with more precise techniques such as TS/ GNSS/TLS; and with repeated measurements of points on the same surface (Brasington et al., 2000(Brasington et al., , 2003Lane et al., 2003;Martínez-Casasnovas et al., 2004;Kaiser et al., 2014;Wang et al., 2016;Fernández et al., 2016;2020a-b;Koci et al., 2017).
Meanwhile, the influence of climate, especially of the rainfall regime, in erosion processes is well known (Poesen et al., 2003;Li and Fang, 2016). Thus, the close relationship between high intensity of precipitations and water erosion is due to: the high erosivity of raindrops in strong storms that produces laminar erosion (Mohamadi and Kavian, 2015); or the runoff of water that cannot infiltrate and then accumulates in the gully channels, removing the soil in them (Poessen et al., 2003). Then, in a global climate change scenario in which intense rainfall events are expected in coming decades (IPCC, 2013;Monjo et al., 2016), the increase of erosion rates can lead to important environmental problems (Li and Fang, 2016). The increase of erosion rates associated to climate changes have been demonstrated along the Late Cenozoic (Mutz et al., 2018).
This paper aims to the development of a methodology for the identification, mapping and quantification of gully erosion for a 40-years period  in an olive grove area. It is based on aerial (both conventional and UAS) photogrammetry and LiDAR data. For each survey, the DSMs and orthophotographs were obtained and the DSMs of differences (DoDs) were also computed, in order to estimate the height differences and volumes between models, as well as their corresponding rates for the different periods considered. Finally, the relationships with rainfall regime along the time have been analysed.

Study area
The study area corresponds to an active gully stretch with an extension of 1.81 ha in a catchment area in olive groves. It is located in the province of Jaén (southern Spain) at a distance of about 25 km from the province capital ( Figure 1). The altitude of the surrounding basin ranges between 430 and 465 m with an average slope of 9°. The basin belongs to the natural region of the Eastern Guadalquivir river basin. Geologically, it is made up of the Guadalquivir Units (Pérez-Varela et al., 2017), a set of materials, intercalated by means of tectonic structures in sedimentary rocks from the Miocene age ( Figure 1). Lutites, evaporites and carbonates of Triassic age are the predominant lithologies as well as marls and clays of Cretaceous-Paleogene. In the study area, Triassic lutites and sandstones outcrop, partially covered by alluvial and colluvial Quaternary deposits.
The area is affected by intense erosion processes, both laminar and gully, in addition to other surface processes such as landslides (Fernández et al., 2016). Some sections of the gully area, which sometimes affect rural roads and paths, are shown in Figure 2. Erosive processes are very intense in the Jaén province, causing a reduction of soil productivity and other effects such as filling of reservoirs and damage to infrastructures (Calero et al., 2018). Figure 2 shows the area and the main and a tributary gully with the characteristic V-shape as well as soil slabs falls at the steep sidewalls. Moreover, upstream erosion is affecting to a rural path, which can cause its collapse in coming years as shown in the photographs from 2016 to 2019 in Figure 2c.

Materials
The methodology is based on photogrammetric techniques, accomplished by LiDAR data processing. A set of eight historical aerial flights (1980, 1996, 2001, 2005, 2009 2011, 2013 and 2016) acquired from different cartographic agencies have been used. The data source and characteristics of these flights are shown in Table 1. The image ground sample distance (GSD) varies from of 0.27 to 0.45 m. The first flights are analog and the last ones are digital and color or color-near infrared (RGB-NIR). Meanwhile, the LiDAR data correspond to the first national LiDAR coverage, which was captured in 2014 in the Andalusia region with a resolution of 1 point/m 2 .
Images captured with two UAS flights (2019 and 2020) have been also considered. The UAS employed is a DJI Phantom 4 with a RTK module integrated that provides positioning of centimeter accuracy. It has allowed the dramatic reduction of the number of the field-surveyed ground control points for the photogrammetric processing, the improvement of the flight security with a flight range up to 20 minutes and the data collection for post processing kinematic (PPK). The camera is a DJI FC6310R (20 mpx and 0.0024 mm pixel size) with a 8.8 mm lens. Although the UAS carries a GNSS antenna for vehicle positioning/navigation, an additional GNSS equipment (LEICA SYSTEM 1200+) was employed for ground control and check points (GCP/CHK) field measurements.

Methodology
The methodology is based on aerial and UAS photogrammetry techniques. It has been developed in previous works of the research group for landslides (Fernández et al., 2016(Fernández et al., , 2017 but it has been adapted to the erosion studies (Fernández et al., 2020 a-b). It can be summarized in the following steps: 1. Image acquisition and field work. 2. Image processing and orientation. 3. Generation of DSMs and orthophotographs, and calculation of DSMs of differences (DoDs). 4. Estimation of height differences and volumes between models. 5. Analysis of rainfall time series and relationships between erosion rates and rainfall regime.

Image acquisition and field work.
Both aerial photogrammetric and LiDAR data are available from several public Spatial Data Infrastructures (SDI) and download services (Fernández et al., 2020a). Thus, photographs of the national flights were downloaded from the photo-library of the National Geographic Institute of Spain (IGN) and photographs of the regional flights were downloaded from the photo-library of Andalusia (Institute of Statistics and Cartography of Andalusia, IECA). The 2005 and subsequent flights are from the National Plan of Orthophotography (PNOA). The LIDAR data can be accessed via the download service of IGN. These data consisted of a previously classified point cloud of 1 point/m 2 density, organised in tiles of 2 x 2 Km.
Additionally, two UAS flights of very high resolution were made on 03-April-2019 and 13-February-2020 (Fernández et al., 2020b) (Table 1). The flight planning was conducted with the DJI desktop software. A GSD of 2-2.5 cm was selected for the images, captured from a flying height above terrain between 90-120 m, within the limits of Spanish regulations for RPAS. Each of the photogrammetric projects consisted in a block of vertical images organized in three strips adjusted to the gully area. For the first flight of 2019 a set of 30 ground control (GCPs) and check points (CHK) was measured in the field by means of differential GNSS.

Image processing and orientation.
First, the LiDAR point cloud was obtained by merging the 2 x 2 km tiles. In order to control its quality, a set of 21 check points were selected in the point cloud and measured in the field with differential GNSS. These points were measured in constructions and building roof corners, unequivocally identifiable in the LiDAR data. They were fixed by adjusting a flat surface to the LiDAR points located on the roofs and computing the minimum bounding rectangle (Fernández et al., 2020a). The mean X and Y errors calculated were lower than 0.20 m being the RMS and SD of about 0.5 m. Meanwhile, the mean vertical (Z) error was of 0.24 m, the SD was of 0.40 m and the RMS was of 0.46 m.
In addition, a point data set was extracted from the LiDAR point cloud to be used as second-order GCPs and CHK. Errors estimated in the LiDAR data allowed to validate this methodology. Since all data were available from public download services, the approach could be applied without GCP measurement in the field. These GCPs, the flight and camera parameters, and a number of tie points were used in the orientation of the historical flights, with the software Socet Set 5.6, in a digital photogrammetric workstation. Meanwhile, the UAS images were processed with Agisoft Metashape. The first flight of April 2019 was the reference flight. This flight was aligned with RTK camera positions (RTK c) and GCPs measured in the field with DGNSS. Therefore, a network of four GCPs located at the ends of the strips was selected. The remaining 26 points were considered as check points. Table 2 shows the orientation errors of the reference 2019 flight, estimated at these check points (RMS CHK). It also shows the errors in camera positions (RMS RTK c) after the block adjustment. After validating the alignment of this reference flight, some additional points were measured as second order GCPs in the images. They were well-defined and stable points that in that way could be transferred to the 2020 flight. Therefore, a common reference system (CRS) was kept for both flights without any necessity of additional field surveyed GCP for the 2020 flight. The results of the orientation of this flight are also included in the  Table 2. Orientation processing data and errors (in m).

Generation of DSMs and orthophotographs.
Digital surface models (DSMs) were generated for the historical flights, with a resolution of 2.5 m using automatic correlation techniques (dense matching) implemented with the NGATE module of Socet Set 5.6. From these DSMs, grid files were exported to a GIS software (QGIS) and then transformed to raster in TIFF format. After Fernández et al.(2016Fernández et al.( , 2020a, DSMs were used due to the inconveniences to generate DTMs. The DSM uncertainties were the Z errors calculated at CHK computed in the orientation process (Fernandez et al., 2020a). Meanwhile, the orthophotographs were exported to a raster format (TIFF), with a resolution equal to the original imagery GSD (0.5 m).
For UAS flights, DSMs were also generated from the dense point clouds in Agisoft Metashape. DSM resolution was of 0.10 m with an uncertainty also of 0.10 m (Fernández et al., 2010b). However, DSMs were resampled to a pixel size of 2.5 m to make them comparable with historical DSMs and the uncertainty of these was assumed. Next, orthoimages were generated with a resolution of 0.05 m, but they were also resampled to the resolution of historical flights (0.5 m).
Then, DSMs of Differences (DoDs) were calculated from DSMs, which allowed the identification of the areas in which vertical changes of the ground surface occur in the periods considered. DSMs differences may be negative or positive, that allows to distinguish areas of ground descent (soil depletion or erosion) or ascent (soil deposition or accumulation), respectively. Moreover, the vertical uncertainties of the DoDs are estimated as follows (Brasington et al., 2000): Unc. DoD YEAR1-YEAR2 = (Unc. DSM YEAR1 2 + Unc. DSM YEAR2 2 ) 0.5 The uncertainties estimated for DoDs were around 1 m. Then, a threshold filter of ± 1 m was applied for identification and delimitation of the areas affected by erosion processes, discarding the areas with DoDs values lower than 1 m.

Estimation of height differences and volumes between models.
The GIS analysis of DoDs allows the estimation of the soil depletion (erosion) or deposition (accumulation). First, in the GIS we delimitated manually the gully area based on the filter applied to DoDs. Thus we considered as gully those areas with overall depletion higher than 1 m, which have to be validated through interpretation of orthophotographs. Then, the depletion and deposition areas within the gully area delimitated could be calculated in order to estimate the relative importance of erosion and accumulation processes in these zones.
Meanwhile, the calculation of the mean or average values from the DoDs or height differences in the gully area allowed us to determine more accurately the balance between depletion and deposition processes. The depletion processes predominate when the balance is negative, while the deposition processes do when it is positive. An approach based on the probability distribution of the errors were addressed (Brasington et al., 2000;Lane et al., 2003;Fernández et al., 2020a). Moreover, the average values of negative (depletion) and positive (deposition) height differences were also calculated. From these, the rates of depletion, deposition and balance of the height differences could be estimated by dividing the corresponding values by the time period between models.
Finally, the estimation of the depletion, deposition and balance volumes of soil material were made following the same approach based on probabilities mentioned before. Soils losses occur if the balance is negative or gains if positive. In the same way, the volume rates could also be calculated. Finally, the volume rates were transformed to mass rates by surface and time (t/ha*year), that is, the units in which the erosion data are usually expressed, dividing by the area in ha and considering an average soil density of 1.5 t/m 3 .

Analysis of rainfall time series and relationships between erosion rates and rainfall regime.
First, climatic data were obtained from the ERA-Interim re-analysis database of the European Center (ECMWF) to which downscaling techniques were applied to obtain grid resolutions of 10 km from 1970. Moreover, the data available in different meteorological networks have been used to make a more precise reconstruction of the precipitation values for each point, reaching resolutions of 1 km through geostatistical techniques. In particular, data from the Meteorological State Agency (AEMET), the Environmental Department of the Andalusian Government (phytosanitary network) and the Hydrological Information System of the Guadalquivir Basin have been used.
From daily precipitation data of the nearest grid point to the study area, several statistics have been calculated in order to analyse the rainfall regime in the last fifty years. Thus, besides the mean (Me), standard deviation (Std), maximum and minimum of the temporal data, the number of days exceeding some extreme rainfalls were computed. The same analysis was performed for the accumulated rainfalls of 7 days (week), 1 month and 3 months. Finally, the relationships with erosion and deposition rates were analysed in a qualitative way.

RESULTS
The DoDs for the different periods are shown in Figure 3. The calculations obtained from them are described in next sections.

Areas affected by erosion processes
The areas affected by the erosion processes are shown in Table  3. Thus, in this particular gully stretch, depletion processes predominated clearly (57%) over deposition processes (7%), being the uncertainty area of 36% for the whole period. These percentages were not constant along the different periods analysed. In the periods 2009-2011 and 2011-2013, depletion area reached values around 40%, being the uncertainty area around 50-60%, while in the remaining periods the uncertainty area reached values higher than 80%. Anyway, most periods presented a predominance of erosion over deposition except in 1980-1996, 2013-2016 and 2016-2019.

Height differences
The height differences are shown in Table 4 (Fernández et al, 2020a).

Rainfalls
The rainfall data are shown in Tables 6-8 and Figure 4. From Table 6, it can be observed that the mean annual rainfall was 470 mm with a standard deviation of 158, ranging from 230 and 903 mm. In general a slightly increase (over 10%) of the main rainfalls were observed from values under 500 mm in the first 5-years periods to values over 500 mm in the last periods (except for the very last). The maximum values accumulated in 1 day (24 h), 7 days (a week), 1 month and 3 months also showed some increase in the last 5-years periods regarding the first ones. This was near 30% in 1-day rainfalls and between 5-10% in the other accumulated rainfalls. Regarding monthly rainfalls, the average was of 39.8 mm with a high Std of 41.5, ranging from 3.2 mm in July to 65.3 mm in December.  Table 6. Rainfall data (mm). AnnMe (annual means), maximum in 1 day, 7 days, 1 month and 3 months (accumulated).  : a: 1980-1996; b: 1996-2001; c: 2001-2005; d: 2005-2009; e: 2009-2011. In the same way, the number of days in which a rainfall amount has been exceeded was higher in the last 5-years periods than in the first ones. Thus, the number of days with precipitations higher than 50 mm in 1 day reached a maximum value in 1996-2000 and 2011-2015 (Table 7). The total number was 12 days with a return period of about 4 years. Meanwhile, the number of days with accumulated precipitations higher than 100 mm in a week, 200 mm in a month and 400 mm in 3 months was also maximum in 1996-200, 2006-2010 and 2011-2015. The total number of days were of 32, 164 and 184, respectively. Taking into account the periods between photogrammetric surveys, the periods with maximum number of days with a high amount of rainfalls were 1996-2001, 2009-2011 and 2011-2013, while for a week, a month and 3 months, the rainy days were maximum especially in 2009-2011 (Table 8).  Table 7. No. days in which a rainfall (mm) has been exceeded.

Periods of 5 years
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume VI-3/W1-2020, 2020 Gi4DM 2020 -13th GeoInformation for Disaster Management conference, 30 November-4 December 2020, Sydney, Australia (online)  Table 8. Number of days in which a rainfall (mm) has been exceeded for the periods considered (photogrammetric surveys).

DISCUSSION
Regarding the uncertainties, first, the RMS errors of the Z component propagated from LiDAR data to photogrammetric orientation of aerial flights (ranging from 0.5 to 0.9 m) are assumed as the uncertainties of the corresponding DSMs. This leads to DoDs uncertainties of around +1 m, in the same order of magnitude of those estimated in previous works (Lane et al., 2003;Fernández et al., 2017;2020a). Regarding the uncertainties of DSMs from UAS flights, these were derived from the RMS errors obtained in the alignment process. They are lower than 0.05 m in Z, similar to the estimated in previous works (Fernández et al., 2016(Fernández et al., , 2020b. However, as the UAS products were resampled to the lower resolution of the other ones, their uncertainties (+1 m for DoDs) have been adopted.
The gully stretch presents an intense activity, especially of erosion or depletion processes respect to the deposition ones (the depletion area was about 7-8 times the deposition area). Thus, the balance was about -1.6 m for the average of height differences (incision) and -30000 m 3 for the volumes (soil losses), which lead to rates of -0.05 m/year and -750 m 3 /year (-44 t/ha*year), respectively. The stretch belong to a gully system that, in a wider area of 7.45 Km 2 , has rates of -0.015 m/year for average height differences and -2400 m 3 /year for volumes (-5 t/ha*year) (Fernández et al., 2020a). The values for the gully stretch considered are within the range found in the region and in the references (Poesen et al., 2003;Hayas et al., 2017), while for the wide area they are in the lower part of this range.  The volume analysis shows even clearer results, since the depletion and balance mass rates of the more active periods are higher in an order of magnitude (330-450 t/ha*year in absolute terms) than the other ones (below 40 t/ha*year). These values are higher than those usually found in the references, near the most extreme cases (Lane et al., 2003;Martinez-Casasnovas et al., 2004). The importance of erosion processes can be observed in the maps of Figure 3 (e-f). Moreover, the depletion processes that until the 2009-2011 period were concentrated in the gully center, from 2011-2013 displaced towards its lateral walls. Thus, the evolution was first by means of vertical incision and later by means of gully walls retraction which was confirmed by detailed studies with UAS (Fernández et al., 2020 b).
The relationship between rainfall and erosion has been well established in different environments all over the world (Poesen et al., 2003) as well as in Andalusia (Hayas et al., 2017). In this study, these relationships were found in a qualitative way. Thus, the higher activity in 2009-2011 and 2011-2013 seems to be related to the increase of rainfalls in these periods, for the precipitations in 24 h as well as for the precipitations in 7 days, 1 month and 3 months. In the autumn-winter of 2009-2010, 2010-2011 and 2012-2013, several days had rainfalls higher than 30 mm, and there were events of weekly rainfalls higher than 100 mm, monthly rainfalls higher than 200 mm and 3months rainfalls higher than 400 mm. However, the appreciable rainfalls of the period 1996-2001 was not reflected in the erosional activity in the gully stretch (Fernández et al., 2020a).
More difficult is to conclude about the influence of climate change on the erosional activity, mainly due the lack of enough evidences of change in the rainfall regime. It can be observed an increase of the precipitations, especially of the 1 day rainfall (about 30%), but also in 1 week, 1 month and 3 months (5-10%), but probably it is necessary to have longer datasets (about 100 years) and especially analyses of the rainfall regime in more points scattered in a wider area. Moreover, some observations such as the existence of wet periods with low erosional activity (1996)(1997)(1998)(1999)(2000)(2001) lead to take into account the possible influence of other factors, such as the construction of rural trails, paths and roads and the changes in land use and practices (Fernández et al., 2020a). Then, we need additional data in order to extract clearer conclusions about the regional influence of all these factors, and especially the role of climate change in the erosion processes.

CONCLUSIONS
We have developed an approach to study the erosion processes in an active gully stretch of about 1.81 ha using aerial and UAS photogrammetry techniques combined with LIDAR data. Data were acquired from public data servers and also captured with ad hoc UAS flights. The image processing has been carried out by means of conventional photogrammetric techniques and SfM-MVS techniques (UAS). Digital Surface Models (DSMs) and orthophotographs were obtained with resolutions of 2.5 m and 0.5 m, respectively. DSMs of differences (DoDs) were also calculated, with an accumulated uncertainty of about ±1 m.
The height differences and volumes for depletion, deposition and balance were estimated by GIS procedures. The results were an average depletion of -1.6 m (incision) and a waste volume of -30000 m 3 (soil losses), which lead to rates of -0.04 m/year and -44 t/ha*year, respectively. These rates were very different along the considered periods, reaching the maximum values (near -300-450 t/ha*year) in 2009-2011 and 2011-2013. These rates are related to rainfalls, being the periods with a more intense erosion processes those in which the precipitations are higher. Rainfall regime has undergone an increase in the last years (a 30% for 1 day rainfalls), but it is necessary more data and analysis to conclude about the influence of climatic change in erosion processes. Thus, future works will be lead to improve the accuracy of photogrammetric and LiDAR techniques and their automation to be applied to larger extensions. Then, more advanced techniques of climatic data analysis will be used.