Multitemporal Analysis of Gully Erosion in Olive Groves by Means of Digital Elevation Models Obtained with Aerial Photogrammetric and LiDAR Data

: Gully erosion is one of the main processes of soil degradation, representing 50%–90% of total erosion at basin scales. Thus, its precise characterization has received growing attention in recent years. Geomatics techniques, mainly photogrammetry and LiDAR, can support the quantitative analysis of gully development. This paper deals with the application of these techniques using aerial photographs and airborne LiDAR data available from public database servers to identify and quantify gully erosion through a long period (1980–2016) in an area of 7.5 km 2 in olive groves. Several historical ﬂights (1980, 1996, 2001, 2005, 2009, 2011, 2013 and 2016) were aligned in a common coordinate reference system with the LiDAR point cloud, and then, digital surface models (DSMs) and orthophotographs were obtained. Next, the analysis of the DSM of di ﬀ erences (DoDs) allowed the identiﬁcation of gullies, the calculation of the a ﬀ ected areas as well as the estimation of height di ﬀ erences and volumes between models. These analyses result in an average depletion of 0.50 m and volume loss of 85000 m 3 in the gully area, with some periods (2009–2011 and 2011–2013) showing rates of 10,000–20,000 m 3 / year (20–40 t / ha*year). The manual edition of DSMs in order to obtain digital elevation models (DTMs) in a detailed sector has facilitated an analysis of the inﬂuence of this operation on the erosion calculations, ﬁnding that it is not signiﬁcant except in gully areas with a very steep shape. data are available from public download services, the approach used in this study can be applied without the need of using field-measured GCP/CHK points, and then, historical flights can be oriented using only the public LiDAR dataset. Moreover, the use of these second-order GCPs [57,58], extracted from the LiDAR data, has guaranteed that all flights were processed under the same common reference system (ETRS89-UTM 30N). The GCPs, the flight and camera parameters (camera GNSS positions, inertial IMU data, camera calibration parameters, etc.) and a number of tie points were used in the orientation or alignment process of the historical photogrammetric flights, using a digital photogrammetric workstation (DPW) and the software Socet Set 5.6 [71]. Due to the lack of calibration data in the old analog film-based cameras, approaches for the use of historical flights must be considered [72,73].

A crucial aspect in these studies is the accuracy or quality of the data collected with the different techniques, since the validity of the results is strongly linked to it. Accuracy analysis can be made from qualitative visual observations but also from different quantitative methods [16,32,41]. In some cases the accuracy can be estimated from errors (usually the root mean square or RMS) of the image orientation or the laser point cloud registration processes, at ground control and/or check points [14,21,23,32,[39][40][41][42][43]47,48,57,58]. Other accuracy analyses are those based on the comparison of the height of profile points or DEMs, with respect to points measured with a more accurate method such as the TS/GNSS [15,20,23,37,39,42,45] or TLS [16,25,38]. Finally, methods based on the comparison of repeated measurements of the same surface can be found, generally by calculating the standard deviation of these measures in a sample of points [13,19,26,35,57]. Moreover, sometimes the error analyses extend from the DEMs to the DoDs derived from them [13,19,30,57,58].
This paper aims to describe a methodology for the identification and quantification of gully erosion through a long period , based on aerial photogrammetry and LiDAR techniques. For this objective, only those data available on public spatial data infrastructures (SDI) or download servers have been used. Ground control points (GCPs) from the LiDAR point clouds were extracted, allowing the alignment of historical photogrammetric flights in the same common coordinate reference system without the need for field surveyed GCPs, although some check points were measured with GNSS in the field in order to monitor the whole process. From these flights, DSMs, DTMs and differential models (DoDs) were generated in order to analyze the effects of gully erosion on olive groves in a selected study area of 7.5 km 2 in the province of Jaén (Southern Spain). This methodology has been validated by means of the estimation of errors and uncertainties from the models obtained. The analysis has been addressed calculating the height differences and volumes between models and their corresponding rates along the different periods considered. Furthermore, a detailed area where the erosion processes are more intense was considered to monitor the gully evolution and to study the relationships with the rainfalls as triggering factor. Finally, to complete the validation analysis and limitations of the methodology employed, DSMs were edited to obtain DTMs, comparing the results obtained in both approaches.

Study Area
The study area, with an extension of approximately 7.5 km 2 , is located in the western part of the province of Jaén (Andalusia, Spain), at a distance of about 25 km from the province capital (Figure 1a,b). It has an altitude of between 355 m and 565 m and an average slope of 5.37 • . It is located within the natural region of the Eastern Guadalquivir river basin, which includes the headwater areas from Cazorla and Segura Ranges. Discarding the Betic Mountains to the south, the basin varies in altitude from about 200 to 1000 m, with slopes generally below 10 • .
From the geological point of view, this basin is made up of the Guadalquivir Units [59], a set of materials of diverse lithology tectonically intercalated with autochthonous loamy-clay marine sediments from the Miocene age ( Figure 1b). The predominant allochthonous lithologies are Triassic lutites, evaporites and carbonates, as well as Cretaceous-Paleogene marls and clays. Specifically, in the study area, Triassic and Miocene materials outcrop (Figure 1c). Triassic units comprise lutites and The province of Jaén is the area with the largest extension of olive groves in the world. It occupies 48% of the arable area of the province [62] and represents 59% of the surface of this crop in Spain, 30% in Europe and 19% in the world [63]. The olive grove connects the economic, cultural and social life of the province, being also its most emblematic landscape. However, its sustainability is exposed to several threats, mainly due to the intensification of the crop, which has reduced the quality of the the agricultural production along with environmental protection and human health and well-being. Currently water erosion is the main cause of soil degradation in the olive grove, as happens in the province of Jaén whose affected area has become the largest in Andalusia and Spain. Thus, 80% of the Jaén olive grove has erosion rates above 10 t/ha*year, with an average value of 32 t/ha*year [65],

Materials
The methodology focuses on the use of photogrammetry techniques with data captured from conventional aerial platforms with sub-metric resolution, accomplished by LiDAR data. A set of historical flights from 1980 [67], 1996   The province of Jaén is the area with the largest extension of olive groves in the world. It occupies 48% of the arable area of the province [62] and represents 59% of the surface of this crop in Spain, 30% in Europe and 19% in the world [63]. The olive grove connects the economic, cultural and social life of the province, being also its most emblematic landscape. However, its sustainability is exposed to several threats, mainly due to the intensification of the crop, which has reduced the quality of the soils in recent decades [64]. Thus, the traditional alternation of Mediterranean crops between cereal and olive is being replaced by a monoculture of the olive grove, affecting the soil´s ability to support the agricultural production along with environmental protection and human health and well-being. Currently water erosion is the main cause of soil degradation in the olive grove, as happens in the province of Jaén whose affected area has become the largest in Andalusia and Spain. Thus, 80% of the Jaén olive grove has erosion rates above 10 t/ha*year, with an average value of 32 t/ha*year [65], which means losses far exceeding the capacity for natural regeneration of the soil [7]. Erosive processes also represent an economic loss, especially related to the reduction of soil productivity but also with other effects such as the filling of reservoirs and damage to infrastructures [66].

Materials
The methodology focuses on the use of photogrammetry techniques with data captured from conventional aerial platforms with sub-metric resolution, accomplished by LiDAR data. A set of historical flights from 1980 [67], 1996 and 2001 [68], 2005, 2009 2011, 2013 and 2016 [67] obtained from different cartographic agencies have been processed (Table 1) In addition, the LiDAR data correspond to the first national LiDAR cover, which in Andalusia was captured in 2014 with a resolution of 1 point/m 2 . Both photogrammetric and LiDAR data are available from several public Spatial Data Infrastructures and download services. Thus, photographs of the national flights (Interministerial and PNOA flights) are available with different image formats (TIFF/JPG/ECW) in the photo-library of the National Geographic Institute of Spain (IGN) [67], and the photographs of the regional flights in the photo-library of Andalusia [68]. Finally, the LIDAR data, stored as tiles (2 x 2 km) in ASPRS LAS format, can be accessed via the download service of IGN [69]. These data consisted of a previously classified point cloud. Raw LiDAR data were not available in the download service.

Methodology
The methodology summarized in Figure 3 consists of several phases partially described in previous papers [57,58,70]: Orientation of the historical flights.

2.
Generation of the DSMs and orthophotographs.

4.
Delimitation of gully areas and DSMs edition of detailed area.

5.
Estimation of height differences and volumes between models.  Thus, the corners were obtained by adjusting a flat surface to the LiDAR points located on the roofs and computing the minimum bounding rectangle (MBR). Figure 4 shows this process. Then these corners were compared with the field-measured ones. Table 2 shows the results of this comparison. The overall mean X and Y errors are lower than 0.20 m while the root mean square (RMS) and standard deviation (SD) are about 0.50 m; combined, mean and RMS horizontal (XY) errors are of about 0.7 m. Meanwhile, vertical (Z) errors are the following: mean, 0.24 m; SD, 0.40 m and RMS, 0.46 m. These results validate the methodology developed to obtain ground control points in order to align and georeferenced the aerial images. Thus, from this point cloud, a set of additional photogrammetric ground control points (GCP) and check points (CHK) was extracted according to the following criteria: • The points should be well distributed throughout the study area. • They should be unequivocally identifiable features (artifacts and constructions). • They should remain stable and visible in most of the flights considered. Although corners and identifiable details could have been measured accurately in the images and the field to be used as conventional GCPs, these points were extracted from the LiDAR data, as indicated in the preceding paragraphs. Since all data are available from public download services, the approach used in this study can be applied without the need of using field-measured GCP/CHK points, and then, historical flights can be oriented using only the public LiDAR dataset. Moreover, the use of these second-order GCPs [57,58], extracted from the LiDAR data, has guaranteed that all flights were processed under the same common reference system (ETRS89-UTM 30N).
The GCPs, the flight and camera parameters (camera GNSS positions, inertial IMU data, camera calibration parameters, etc.) and a number of tie points were used in the orientation or alignment process of the historical photogrammetric flights, using a digital photogrammetric workstation (DPW) and the software Socet Set 5.6 [71]. Due to the lack of calibration data in the old analog filmbased cameras, approaches for the use of historical flights must be considered [72,73].

Orientation of the Historical Flights
First, the LiDAR point cloud was obtained by merging the corresponding 2 x 2 km tiles, downloaded from the above-mentioned service. Data consisted of a classified point cloud of 1 point/m 2 density. Thus, both artefacts and noise had previously been corrected, which was confirmed by means of a visual inspection. Nevertheless, in order to control the quality of the LiDAR data and to develop an approach to extracting reliable control points from the LiDAR point cloud, a set of 21 check points well defined in the point cloud was measured in the field with differential GNSS techniques. These points had to be unequivocally identifiable features (for example constructions and building roof corners). Thus, some agricultural buildings and sheds with flat roofs were selected, and the roof corners were measured with GNSS techniques. Since these corners could not be accurately identified in the LiDAR point cloud due to the low data density, the coordinates of the roof corners had to be computed from the LiDAR data themselves.
Thus, the corners were obtained by adjusting a flat surface to the LiDAR points located on the roofs and computing the minimum bounding rectangle (MBR). Figure 4 shows this process. Then these corners were compared with the field-measured ones. Table 2   These results validate the methodology developed to obtain ground control points in order to align and georeferenced the aerial images. Thus, from this point cloud, a set of additional photogrammetric ground control points (GCP) and check points (CHK) was extracted according to the following criteria: • The points should be well distributed throughout the study area.

•
They should be unequivocally identifiable features (artifacts and constructions).

•
They should remain stable and visible in most of the flights considered.
Although corners and identifiable details could have been measured accurately in the images and the field to be used as conventional GCPs, these points were extracted from the LiDAR data, ISPRS Int. J. Geo-Inf. 2020, 9, 260 8 of 30 as indicated in the preceding paragraphs. Since all data are available from public download services, the approach used in this study can be applied without the need of using field-measured GCP/CHK points, and then, historical flights can be oriented using only the public LiDAR dataset. Moreover, the use of these second-order GCPs [57,58], extracted from the LiDAR data, has guaranteed that all flights were processed under the same common reference system (ETRS89-UTM 30N).
The GCPs, the flight and camera parameters (camera GNSS positions, inertial IMU data, camera calibration parameters, etc.) and a number of tie points were used in the orientation or alignment process of the historical photogrammetric flights, using a digital photogrammetric workstation (DPW) and the software Socet Set 5.6 [71]. Due to the lack of calibration data in the old analog film-based cameras, approaches for the use of historical flights must be considered [72,73].   Table 3 shows the main parameters of the photogrammetric process and the errors involved, both in the GCP and in the CHK points. The XY errors range between 0. Taking into account that the LiDAR control points used in the orientation process are affected by the errors shown in Table 2   Taking into account that the LiDAR control points used in the orientation process are affected by the errors shown in Table 2, the propagation error can be computed in the check points by means of Equation (1):

Generation of DSMs and Orthophotographs
Digital surface models (DSM) and orthophotographs for each flight were generated using automatic correlation techniques (dense matching) implemented with the NGATE (Next Generation Automated Terrain Extraction) module of Socet Set 5.6 [71,74]. A DEM resolution of 5 times (2.5 m) the value of the original imagery GSD (0.5 m) has been considered, according the conventional norm in photogrammetry of using a multiple value of GSD for DEM generation. At present, this rule is gradually being reduced due to the improvements on new algorithms for dense matching and the use of patterns of photogrammetric flights with higher overlaps, as those planned with unmanned aerial systems (UAS). From the DSMs generated, grids files with point spacing of 2.5 m have been exported in text format. These grids were opened in a GIS software, QGIS 3 [75], and transformed to raster in TIFF format with the same resolution (2.5 m), by means of the rasterization tool of QGIS. Thus, the grid values resulting of the DSM generation were not modified by other raster interpolation methods. Meanwhile, the orthophotographs were exported to a raster format (TIFF), with a GSD equal to the original imagery GSD (0.5 m). The DSM for the LiDAR data and the orthophotographs for the 1980 and 2016 flights are shown in Figure 5.
In this paper, the propagated Z errors calculated at the check points are assumed to be the DSMs uncertainties (Table 4), as will be discussed later.

Calculation of DSMs of Differences (DoDs)
Differential models or DSM of differences (DoDs) have been calculated from the DSMs, which has objectively allowed the detection of the areas that undergo vertical displacements of the ground surface between successive dates. Displacements may be negative or positive, depending on whether each model compared with a reference model lies below or above it, allowing the identification of areas of ground descent (mass erosion or depletion) or ascent (mass accumulation or deposition), respectively.
In this paper as in previous ones [57,58,70], DSMs have been used preferably to DTMs. This is because all the models used for computing DoDs are of photogrammetric origin, and thus, the processes of automatic classification and filtering of point clouds do not ensure good results in areas where there is a certain vegetation cover or other elements. Moreover, editing models in order to obtain DTMs in a wide area would be a very time-consuming process and would not ensure an improvement. Nevertheless, a small area has been edited in order to compare and calibrate the results, as described below. The work with DSMs can produce a distortion in the detection and quantification of the gullies, because some vertical changes are due to vegetation (growing, clearing, etc.), constructions or other elements.

Calculation of DSMs of Differences (DoDs)
Differential models or DSM of differences (DoDs) have been calculated from the DSMs, which has objectively allowed the detection of the areas that undergo vertical displacements of the ground surface between successive dates. Displacements may be negative or positive, depending on whether each model compared with a reference model lies below or above it, allowing the identification of areas of ground descent (mass erosion or depletion) or ascent (mass accumulation or deposition), respectively.
In this paper as in previous ones [57,58,70], DSMs have been used preferably to DTMs. This is because all the models used for computing DoDs are of photogrammetric origin, and thus, the processes of automatic classification and filtering of point clouds do not ensure good results in areas where there is a certain vegetation cover or other elements. Moreover, editing models in order to obtain DTMs in a wide area would be a very time-consuming process and would not ensure an improvement. Nevertheless, a small area has been edited in order to compare and calibrate the results, as described below. The work with DSMs can produce a distortion in the detection and quantification of the gullies, because some vertical changes are due to vegetation (growing, clearing, etc.), constructions or other elements.
Regarding the vertical uncertainties of the DoDs, they are estimated as follows [19,58,76]: The calculated uncertainties for DoDs are around 1 m (Table 4). Based on this uncertainty, firstly a threshold filter of ±1 m was applied in the GIS software for visual identification of the areas affected by erosion processes, discarding those areas with vertical displacements lower than 1 m in absolute value. As will be discussed later, according to [13,19,28,45] this filter eliminates 68% of the uncertainty.

Delimitation of Gully Areas and DSM Edition of Detailed Area
The identification and delimitation of gully areas have been based on a semiautomatic procedure. First the application of the filter defined before, and especially the areas with negative values higher than 1 m (in absolute terms), allows the identification of potential gullies, which will be validated through interpretation of orthophotographs. These areas have been delimited by polygon-drawn tools in the GIS software. The gully network shown in Figure 1c results of the application of this procedure to the complete period analyzed .
Moreover, a detailed sector with an intense progress of the gully in recent years has been edited manually in order to obtain DTM by means of Socet Set, using the stereo tools of the DPW. Specifically, the edition has been applied to the DSMs for which the DoDs were more relevant (2009-2011 and 2011-2013), and it implied vegetation clearing and artifacts deletion. Several tests have been carried out with different resolutions and using or not using break lines in the DTM generation. The results were also exported to grid text files that were opened in the GIS software and transformed to raster in TIFF format by means of the rasterization tool of QGIS. In the same way as before, DoDs are obtained from DTMs.

Estimation of Height Differences and Volumes between Models
The GIS analysis of DoDs allows the estimation of the vertical depletion (erosion) or deposition (accumulation) of soil material. First, the areas affected by gully erosion can be calculated as an estimation of the importance of the process in the whole study area. Moreover, the area of depletion and deposition can be also calculated, considering or not the threshold of uncertainty of ±1 m. These calculations allow a first approach for evaluating the predominance of the depletion or deposition processes.
Meanwhile, the calculation of the mean or average values from the DoD (height differences) in the gully areas allows us to analyze in a deeper way the balance between depletion and deposition processes, estimating the general losses or gains of soil material. If the average balance is negative, the depletion processes predominate, and if it is positive, the deposition processes do. In this case, an approach based on the probability distribution of the errors has been addressed. Thus, according to [13,19,28,45], the average balance of height differences can be estimated by weighting each height difference according to its probability of being free of uncertainty. For instance, height differences greater than ±1 m (±RMS Z) have a probability of 68%, height differences of more than ±2 m (±2RMS Z) have a probability of 95% and so on. The average values of negative (depletion) and positive (deposition) height differences are also calculated, considering in each case the corresponding areas. The rates of depletion, deposition and balance height differences can be estimated by dividing the corresponding values by the time interval between models.
Finally, the volumetric calculations are carried out by means of the GIS software (QGIS), estimating the volume balance of soil material (losses if the balance is negative or gains if the balance is positive) in a given area. In this case, the depletion, deposition and balance volumes are estimated following the same approach based on the probabilities described in the height differences. In the same way, the volume rates can also be calculated by dividing those values by the time interval. Finally, the volume rates have been 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 density of 1.5 t/m 3 . This is a rounded value taken from the bulk density estimated from 4 undisturbed cylindrical soil cores (5 cm in diameter, 5 cm long) in each one of sections a and c (Figure 2), yielding an average density of 1.512 (±0.201) t/m 3 .
The height differences and volumes between models were also analyzed in the detailed area for the DSMs as well as the DTMs obtained by means of a manual edition of the DSMs. The analysis with the DTMs has been made in order to validate the results obtained from DSMs. Thus, some tests have been applied to evaluate the influence of model resolution and the use of break lines in the interpolation.

Results
In this section, the results of the analysis of DoDs obtained from DSMs, both in the whole area and in the detailed area, are presented. For the whole area, three analysis are made: first we analyze the areas, second the height differences and third the volumes for the complete period (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) as well as for the partial periods. For the detailed area, the analysis is focused in the height differences and the volumes, also for the complete and the partial periods. Additionally, for the detailed area, an analysis of DODs from edited DTMs has been made, calculating also the height differences and the volumes.

Areas, Height Differences and Volumes for the Whole Area
The DoDs of the different periods for the whole area are shown in Figure 6. The area affected by gully erosion has an extension of 0.17 km 2 and represents 2.25% of the study area of 7.45 km 2 . Within it, the distribution of depletion and deposition areas for the different periods is listed in Table 5. As shown in this table, the uncertainty area, where the height differences are in the range of ±1 m, is much larger than the depletion and deposition areas in every period. Meanwhile the depletion and deposition areas have a similar extent in most periods (1980-1996, 2001-2005, 2005- The maximum percentage of depletion area regarding the total gully area considered was reached in 2009-2011. Regarding the whole period analyzed, the percentage of depletion area becomes much larger (31%) than the percentage of deposition area (7%). Discarding the uncertainty area (62% of the total gully area), the percentages reach 63% and 37%, respectively.
The results of the height differences between models (depletion and deposition heights) for the whole area are presented in Table 6    The results of volume calculations are shown in Table 7, where it is observed that the volumes involved in the complete period are about 100,000 m 3 for depletion and about 17,000 m 3 for deposition processes, which results in a negative balance of waste material of 85,000 m 3 . The rates are near 2800 m 3 /year for depletion processes and about 500 m 3 /year for deposition, with a balance (waste) rate of almost 2400 m 3 /year in absolute value. Translated to mass rates, the values estimated for the whole period are of almost 8 t/ha*year for depletion processes, 2 t/ha*year for deposition and 5.5 t/ha*year for the balance. However, the volumes and rates vary greatly by periods, and thus, there are periods with larger amounts of material motion such as 2009-2011 and 2011-2013 that present depletion rates of 20,000-30,000 m 3 /year, deposition rates around 8000-10,000 m 3 /year and waste rates of 10000-20,000 m 3 /year. The remaining periods present depletion rates lower than 10,000 m 3 /year, deposition rates lower than 5000 m 3 /year with negative balance rates usually lower than 5000 m 3 /year. The mass rates present in some periods such as 2009-2011 and 2011-2013 values higher than 30-60 t/ha*year for depletion processes, around 13-20 t/ha*year for deposition processes and between 20 and 40 t/ha*year for the negative balance. The remaining periods present rates lower than 10 t/ha*year for  Table 7, where it is observed that the volumes involved in the complete period are about 100,000 m 3 for depletion and about 17,000 m 3 for deposition processes, which results in a negative balance of waste material of 85,000 m 3 . The rates are near 2800 m 3 /year for depletion processes and about 500 m 3 /year for deposition, with a balance (waste) rate of almost 2400 m 3 /year in absolute value. Translated to mass rates, the values estimated for the whole period are of almost 8 t/ha*year for depletion processes, 2 t/ha*year for deposition and 5.5 t/ha*year for the balance. However, the volumes and rates vary greatly by periods, and thus, there are periods with larger amounts of material motion such as 2009-2011 and 2011-2013 that present depletion rates of 20,000-30,000 m 3 /year, deposition rates around 8000-10,000 m 3 /year and waste rates of 10,000-20,000 m 3 /year. The remaining periods present depletion rates lower than 10,000 m 3 /year, deposition rates lower than 5000 m 3 /year with negative balance rates usually lower than 5000 m 3 /year. The mass rates present in some periods such as 2009-2011 and 2011-2013 values higher than 30-60 t/ha*year for depletion processes, around 13-20 t/ha*year for deposition processes and between 20 and 40 t/ha*year for the negative balance. The remaining periods present rates lower than 10 t/ha*year for depletion and deposition processes and usually lower than 2 t/ha*year for the balance, except for the period 1996-2001, which almost reaches 10 t/ha*year.

Height Differences and Volumes for the Detailed Area
The DoDs for the detailed area are shown in Figure 7 for the DSMs and in Figure 8 for the edited DTMs (in the latter, for the periods of higher activity, 2009-2011 and 2011-2013). Meanwhile, the calculations of height differences in the detailed area are shown in Table 8a for the DSMs and in Table 8b for the edited DTMs.
Regarding DSMs (Table 8a)      The results obtained for volume analysis of the detailed area are shown in Table 9a for the DSMs and in Table 9b for the edited DTMs. From Table 9a, it is observed that the depletion and balance (waste) volumes for the whole period (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) are higher than 30,000 m 3 (about a third part of the volumes for the whole area), but the deposition volumes are insignificant. Thus, the volume rate is  The results obtained for volume analysis of the detailed area are shown in Table 9a for the DSMs and in Table 9b for the edited DTMs. From Table 9a, it is observed that the depletion and balance (waste) volumes for the whole period (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) are higher than 30,000 m 3 (about a third part of the volumes for the whole area), but the deposition volumes are insignificant. Thus, the volume rate is around 850 m 3 /year for the depletion processes and the balance that gives a mass rate by surface and time near 50 t/ha*year, being these rates insignificant in deposition processes. By periods, the depletion volumes are larger than the deposition volumes for most periods, which results usually in negative balance (waste volumes). The rates for depletion processes are lower than 1000 m 3   The volumes computed for the edited models (Table 9b) are also very similar to those calculated from the automatic models, again when break lines are not taken into account. Thus, the rates are near

Discussion
In this section we first present a discussion about the LiDAR and the photogrammetric flights' uncertainties. After that, we discuss the results obtained from the height differences and volumes in the whole area and in the detailed area. The relationships of erosion rates by periods to the rainfalls are also analyzed. Finally, for the detailed area, the comparison between the results from DSMs and DTMs with different interpolation conditions is discussed.

Accuracies and Uncertainties
First, a quality control of the LiDAR points cloud was carried out with some field GNSS surveyed points. This control showed that horizontal (XY) mean and RMS errors were around 0.7 m, while the vertical (Z) mean error was of 0.24 m in absolute value, and the RMS and SD errors were about 0.40-0.45 m. From these errors, the uncertainties, estimated as the values of RMS or SD as in previous studies [57,58,70], were of 0.7 m for the horizontal component and 0.5 m for the vertical one. Then the horizontal uncertainty for X and Y would be about 0.5 m for each one, which coincides approximately with half of the LiDAR data resolution and is equal to the orthophotographs resolution. Meanwhile, the mean vertical (Z) error of 0.24 m informs us about a good general fit of the aerial LiDAR point cloud to the ground surface, as the SD or RMS values under 0.5 m do about the good fit of particular points or sectors. Moreover, it is in accordance with some previous studies that have used LiDAR data to study Earth surface processes [15,22,26,28,30,31,33,58,76]. In these papers the accuracy, expressed as an instrumental specification or estimated from comparison of LiDAR data with more precise methods (GNSS, TS, TLS), ranges mostly from 0.1 m to 0.5 m. Lower values can only be obtained by means of TLS surveys, in which the accuracy ranges between 0.02 m and 0.1 m [15,22,25,26].
Regarding the uncertainty of the photogrammetric products (DSMs and orthophotographs), in some previous studies [57,58,77], it has been observed that the vertical uncertainties in DEMs obtained by means of matching from photogrammetric flights are about two or three times the orientation residual errors in Z, calculated in the control points (GCPs). However, in this paper, we have considered a propagation error from LiDAR data to the photogrammetric orientation process using external check points. In general, since an external check has been used, errors have been higher than the estimations of the previously mentioned papers. As can be observed in Tables 3 and 4, the Z error in check points is higher than two or even three times the Z error in control points. Thus, the RMS errors propagated of the vertical component (Z) show values between 0.5 and 0.9 m, which are assumed as the uncertainties of the DSMs. From these, the uncertainties of DoDs are calculated, resulting in values around 1 m, in the same order of magnitude to those estimated in reference studies with similar methodologies and resolutions [20,28,47,48,58]. Higher accuracies are only obtained when images of higher resolution are processed: 10 −3 m order and even lower in close-range photogrammetry [20,34,36]; 10 −2 to 10 −1 m order in terrestrial photogrammetry [23,25,35,[38][39][40] or UAV photogrammetry [21,23,24,[39][40][41][42]. These last accuracies are comparable with those obtained by means of GNSS or TS surveys [13,19,20].
Following the methodology developed in the studies of [13,19,28,45], a minimum level of detection (minLoD) threshold could be defined from the uncertainties in order to perform the height differences or volumes calculations from DEMs in a reliable way. Therefore, a simple approach could be discarding the values lower to the minLOD defined as the SD or RMS errors. In a normal distribution, these values represent 68% of the confidence interval, so the calculations are free of uncertainty in this percentage. Usually, it is considered twice these errors, which represents 95% of the confidence interval, and so is the accuracy of the calculations. However, height differences and volume calculations are very sensitive to this minLOD [45], and a large amount of information can also be lost with this procedure [13]. Thus, there are other approaches based on the idea that each height difference presents a confidence interval and so a probability of being free of uncertainty can be estimated. From the DoD of each period analyzed, a T-score can be computed and then the probability of being free of a given uncertainty [13]. Finally, the calculations of height differences and volume can be done by weighting each height difference by its probability.

Areas, Vertical Heights and Volumes in the Whole Area
From the results section, firstly a study area can be observed affected by a gully erosion of variable intensity along the different periods. Qualitatively, Figure 6 shows that depletion or erosion processes (negative height differences) predominate over deposition in the gully area. This is particularly observed in the map of the complete period (Figure 6h) where the gully system extends and branches throughout the whole area.
The area affected by gully erosion is 0.17 km 2 of a total study area of 7.45 km 2 . Thus, it represents 2.25% of the study area, the gully density being 3.90 km/km 2 . Both are partial data calculated from a gully stretch and an area that does not correspond to a catchment. Meanwhile, the percentage and density data estimated for the entire catchment area of about 20 km 2 (Arroyo del Cortijo de la Piedra, Figure 1) are 1.35% and 3.15 km/km 2 , respectively. These values are in the same order of magnitude as those found in previous studies for catchment areas of similar extension [49] or even larger [21,24,33,43,53,55]. Nevertheless, some extreme cases have been documented where values between 16% and 60% are reached [15,25,47,56], although the areas considered are so different that the results are hardly comparable.
Within the gully area, in turn, the percentage of depleted area for the whole period of analysis (1980-2016) becomes much larger (31%) than the percentage of depositional area (7%), although the uncertainty area, considered here as the interval between ±1 m (±RMS Z), reaches more than 60% (Table 5). This uncertainty area is not an area where erosion and deposition processes do not occur, as it has been discussed before, and a correction based on the error probability has been applied to the calculation of height differences and volumes. However, discarding the uncertainty area, depletion or erosion predominate over the deposition or accumulation of soil material (63% and 37%, respectively). These results coincide with those observed in some studies [47], unlike other ones relative to fluvial systems where both processes are balanced [13,32] or even the accumulation processes predominate [28].
In the same sense, the balance between the average of depletion (0.61 m) and deposition height differences (0.10 m) in the gully area became 0.51 m of depletion ( Table 6). The rates are of 0.017 m/year for depletion, 0.003 m/year for deposition and 0.014 m/year of negative balance, which confirms the predominance of depletion or erosion processes with respect to the deposition or accumulation processes. These heights of decimeter-meter order are comparable to some reference values found in approaches based on aerial photogrammetry and LiDAR [15,28,30,32,45,47,48,78].
In volume, it can be observed that the depletion volumes are generally larger than the deposition volumes, in absolute values, so the negative balance volume for the complete period (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) is about 85,000 m 3 (Table 7) which must be evacuated by the drainage network to which the gullies lead. The rates are 2800 m 3 /year for depletion and 500 m 3 /year for deposition, which produces a negative balance (waste) near 2400 m 3 /year. Considering the whole study area (745 ha) and a bulk density of soils of 1.5 t/m 3 , the rates are of 5.72 t/ha*year, 0.94 t/ha*year and 4.78 t/ha*year, respectively. These rates are in accordance with the values estimated for the province of Jaén [65], taking into account that these rates are for all erosion processes including laminar, rill and gully erosion. Moreover, they are within the range found in the references [10,11,49], although in their lower values. Nevertheless, this range is very wide depending on the area extension and environment, if the catchment or gully area is considered, or how the measurement is performed (for instance, depletion volume or sediment production). Usually, when very large areas of 10 6 ha order are considered [10,43], rates around and even lower than 1 t/ha*year are found. When the areas are between 10 2 and 10 6 ha, the rates are between 1-10 t/ha*year. Finally, areas smaller than 10 2 ha present rates higher than 10 t/ha*year [10,11,78]. However, most study areas with similar approaches and extension as those considered in this study present rates much higher, even when the whole catchment is considered: 39.7 t/ha*year [49], 331 t/ha*year [47], between 160 and 430 t/ha*year [48], around 700 t/ha*year [28], 871 t/ha*year [40] and extreme values of 1500-2500 t/ha*year [78]. If only the gully area (16.75 ha) is considered, the rates estimated in our analysis reach values of 255 t/ha*year for depletion, 42 t/ha*year for deposition and 213 t/ha*year for the (negative) balance, more similar to those found in the references.
By periods, the analysis is much richer, as the different periods present great variations regarding the average values for the whole area. Thus, in the maps of Figure 6 intense depletion or erosion processes can be observed in the periods 1996-2001, 2009-2011 and 2011-2013, especially in the last two. Quantitatively, the higher proportion of uncertainty areas corresponds to periods in which the processes are less intense. Discarding the uncertainty areas, the proportion of depletion areas to deposition areas is balanced in the periods 1980-1996, 2001-2005, 2005- (Table 6).
This coincides with the higher values of average height differences of depletion (decimeter order), which occurs in the same periods. It is also in accordance with the higher (negative) values of average balance, given that the average deposition values present in general low values in every period (centimeter order). The rates even increase the differences between periods, and therefore the maximum average depletion and balance height differences in absolute terms occur in the  (Table 7). Some periods such as the first one (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) present large absolute volumes, both of depletion and deposition. This is due more to the uncertainties in the models obtained from photographs of a worse quality (digitized original film images) than to real differences between DSMs over a long period of time. The rates show more clearly the different erosion activity in the periods with values for balance rates near to 20,000 m 3 (Figure 9). The weekly rainfalls show maximum values (higher than 100 mm) in the autumn-winter of the years 1996-1997, 1997-1998, 1999-2000, 2009-2010, 2010-2011 and 2012-2013. In fact, the absolute maximum value of 150 mm is reached in November 2012 ( Figure 9b). Meanwhile, the monthly rainfalls also present maximum values (around 250 mm) in the winter of 1995-1996, 1996-1997, 1997-1998, 2009-2010, 2010-2011 and 2012-2013 1996-1997 and 1997-1998. The relationship between rainfall and erosion has been well established in different environments all over the world, rainfall being the most important triggering factor of these erosion processes [49]. In fact, in the Andalusian region the analysis of the evolution of drainage networks from historical orthophotography  shows similar patterns to the results of this study [49]. Thus, the maximum values of precipitation in 24 h and the mean annual rainfalls which occurred in 1997 and 2010 generated an increase of the drainage network. These patterns and relationships have been observed in other processes such as landslides in the region [57,70] under an irregular rainfall regime in which dry periods of several years alternate with wet periods concentrated throughout the autumn-winter seasons of two or three consecutive years. Nevertheless, as can be observed in Figure 9 and Tables 6 and 7, the rainfalls related to the period 1996-2001 are approximately 80% of the rainfalls of the period 2009-2011, but their consequences in the erosion processes are about 50%. Thus, some questions about the influence of other factors, in addition to the rainfalls, should be analyzed. Among these the influence of rural trails, paths and roads (as shown in Figure 2c) and changes in land use and practices but also the mechanism of gullies' growth once the process has started. Although the influence of precipitations and changes in the land use have been analyzed in areas near to the study area of this paper [49], additional data is necessary in order to extract conclusions about the regional influence of all these factors.

Areas, Vertical Heights and Volumes in the Detailed Area
In general, all the observations made in the previous section for the whole area are valid in the detailed area, but the heights, volumes and rates are different. Thus, the depletion processes are clearly observed in the gully area (Figure 7). The catchment area covers a surface of 30.83 ha and the gully stretch 1.81 ha, which represents a gully intensity of 5.87%, approximately twice the percentage estimated for the whole area.
In the same way, in this sector, the average height difference for depletion in the complete period 1980-2016 is almost three times (1.70 m) the value computed in the whole area (0.61 m). Meanwhile, the average for deposition is practically insignificant and so the average for the balance is practically equal to the depletion (Table 8a). Then the average rates for the depletion and balance height differences reach values of around 0.05 m/year. This means that in this active sector of the gully, the erosion processes predominate absolutely over the deposition processes. This is confirmed by the values shown by the volume analyses (Table 9a). The average depletion and balance volumes are around 30,000 m 3 , approximately a third of the volume of the gully stretch of the whole study area, while the deposition volumes are practically negligible. The volume rates of depletion and balance reach values near 900 m 3 /year, which translated to mass rates surpasses 40 t/ha*year. These values are significantly higher than the values calculated for the study area and in the middle range of those found in the province and the references [10,11,49]. If only the gully areas are considered, the rates become near to 700 t/ha*year, comparable with some extreme values of the references [28,40,47,78].
As for the whole area, there are periods in the detailed sector with significantly higher values, such as the periods 2009-2011 and 2011-2013. This can be observed in the maps of Figure 7 (especially in Figure 7e,f), where the negative values of the height differences (depletion) reached in the gully areas in these periods predominate over the null or positive values (deposition). Moreover, the resolution of these figures allows us to observe the distribution of the depletion areas that in the period 2009-2011 were concentrated in the center or axis of the gully, while in the period 2011-2013 they were displaced towards the lateral walls of the gully. This led to an evolution first by means of vertical incision (2009-2011) and later by means of gully wall retraction. Quantitatively, these periods present notable rates of depletion and balance height differences (0.50 m/year and 0.30 m/year) regarding the remaining periods in which the rates are insignificant (Table 8a). Even in the period 1996-2001, when rates in the whole area reach a remarkable level, in this particular area, they are practically negligible. The volume analysis is even clearer (Figure 9a), since the depletion and balance rates of these periods are higher in an order of magnitude (8000 m 3 /year and 5000 m 3 /year) than the other ones (always under 800 m 3 /year). The translation of these volumes to mass rates allows extreme values to be reached, up to 275-425 t/ha*year for balance, as those found in the references [28,40,47].
The relationship of the erosional activity in this area to rainfalls leads to similar conclusions as those for the whole area. Thus, the maximum values both in the weekly and monthly rainfalls in the periods 2009-2011 and 2011-2013 are related to the higher erosion and loss of soils in these periods ( Figure 9). However, the notable activity of the period 1996-2001 in the whole area, related to the rainy years 1996-1998, is not found here. This leads us to consider the influence of local triggering factors such as the roads and rural ways and land use and practices that accelerate or inhibit the erosion processes. Regional studies of determinant and triggering factors as well as hazard analysis, which exceed the objective of this study, could lead to a greater knowledge about the causes and evolution of erosion processes in the region.
Regarding the influence of DEMs resolution and data processing, some tests were performed. First, only manual edition to obtain DTMs from DSMs did not produce a significant improvement in the models and calculations, and the increase of resolution (until 1 m) did not either. Meanwhile, the combination of the edition procedure with the introduction of break lines at a resolution of 1 m (Figure 8) allowed the significant improvement of the DTMs definition. Thus, when processing DTMs of the period 2009-2011 (edited and interpolated with break lines), the rates of depletion and balance height differences increase by 50% (from 0.50 m/year to 0.75 m/year, Table 8b). However, these rates do not change in the period 2011-2013 (0.33 m/year). Meanwhile, the rates for depletion and balance volume also increase, in this case by 40% (from 10,000 m 3 /year to 13,500 m 3 /year) for the period 2009-2011, but they do not change either in the period 2011-2013 (Table 9b).
From these data, first a subestimation of the intensity of the erosion processes is observed when only the DSMs are considered and linear interpolation is applied to obtain the models. Thus, the depletion heights and volumes can become 40%-50% higher if the DSMs are edited and break lines are considered in the interpolation, as in the period 2009-2011. This sub-estimation should be taken into account in the quantification of gully erosion in addition to the uncertainties derived from the orientation processes. However, this behavior is not observed in all cases, such as in the period 2011-2013 in which the heights differences, volumes and rates are similar in both edited and non-edited models. These differences in the behavior of the DTMs after edition could be related to the shape of the gully section and thus to the erosion processes involved. Thus, from the observation of profiles or sections of DSMs and DTMs, the DoDs and the orthophotographs (Figure 10), it is observed that the gully growth in depth until the period 2009-2011 is mainly in the vertical component (gully incision). In this case, the resulting shape of the 2011 section is very steep and narrow. This is in accordance with the observation of the height difference map of this period (Figures 7e and 8a), which is higher in the gully central axe. Therefore, in these steep section, shapes linear interpolation can smooth the model and sub-estimate the depth and thus the height differences and the volumes. However, the growth of the gully in the period 2011-2013 is mainly in the horizontal component by lateral walls retraction more than the incision. This leads to a wider and less steep section, where the interpolation does not smooth so much the models, and thus, the differences between the DSMs without break lines and the edited DTMs with break lines are practically insignificant. It also confirms the evolution observed in the height difference maps (Figures 7f and 8b), where the height differences are displaced from the gully central axe to the lateral walls.
Despite the observed sub-estimation of height differences and volumes in some given gully morphologies, the technique applied with automatically generated DSMs is very useful for a first identification and even quantification of gully processes. If more accurate measurements are needed, a model edition is required in some cases. It means the application of time-consuming procedures, which are not always justified for large areas. In any case, the availability of more quality models such as LiDAR datasets and those coming from aerial or UAV surveys of very high resolution and quality would contribute to reducing the errors and uncertainties of models. smooth the model and sub-estimate the depth and thus the height differences and the volumes. However, the growth of the gully in the period 2011-2013 is mainly in the horizontal component by lateral walls retraction more than the incision. This leads to a wider and less steep section, where the interpolation does not smooth so much the models, and thus, the differences between the DSMs without break lines and the edited DTMs with break lines are practically insignificant. It also confirms Despite the observed sub-estimation of height differences and volumes in some given gully morphologies, the technique applied with automatically generated DSMs is very useful for a first identification and even quantification of gully processes. If more accurate measurements are needed, a model edition is required in some cases. It means the application of time-consuming procedures, which are not always justified for large areas. In any case, the availability of more quality models such as LiDAR datasets and those coming from aerial or UAV surveys of very high resolution and quality would contribute to reducing the errors and uncertainties of models. The

Conclusions
The usefulness of geomatics techniques, mainly photogrammetry and LiDAR, for the study of superficial processes like gully erosion and landslides is proved again with success. The accuracy and consistency of these techniques enables the detection of change in detailed terrain features over long periods and large areas.
The methodology developed in this study is based on the extraction of GCPs from LiDAR data that allows the orientation of historical flights in a common coordinate reference system, avoiding or reducing the field GNSS measurements of the GCPs. Then, DMS, DoDs and orthophotographs are obtained and analyzed. In this study all data, both the different images and the LiDAR point cloud, are public and available from different spatial data infrastructures and database servers. This approach allows these studies to be addressed without the necessity for data capture, which is a costly process in time and money.
Regarding the accuracy of the methodology, the propagated RMS error estimated for the LiDAR data and the image orientation process is about 0.80-0.90 m for the horizontal (XY) component and 0.50-0.95 m for the vertical (Z) component. These propagated errors for the vertical (Z) component lead to consider the uncertainties of DSM of Differences (DoDs) of about ±1 m. Such uncertainties are enough to identify and map gully erosion processes and to quantify them through the calculation of height differences and volumes between successive digital elevation models.
The gully area, delimitated from the interpretation of DoDs and orthophotographs, is of 0.17 km 2 , which represents 2.25% of the whole study area. The height differences present an average balance of 0.51 m for the gully area in the complete period, with maximum average rates of depletion around 0.10 m/year in some periods such as 2009-2011 and 2011-2013. Meanwhile, the waste volume in the complete period is about 85,000 m 3 , reaching maximum rates of 10,000-20,000 m 3 /year. It produces mass rates of 20-40 t/ha*year in the periods of higher activity, considering the total study area (745 ha) and a soil density of 1.5 t/m 3 . This different activity of the erosion processes over time can be related to the rainfall regime. Thus, the weekly and monthly rainfalls reach maximum values higher than 100 and 250 mm, respectively, in some episodes included in the periods 1996 The analysis of DTMs, obtained by the manual edition of some DSMs in the detailed area, allows us to extract some ideas about the approach implemented, mainly related to the models' resolution and the shape of the gully. Thus, given that there are no significant changes when only manual edition is considered (even increasing the resolution from 2.5 m to 1 m), this leads to the necessity of adding break lines to the interpolation. When these features are included the depletion heights and volumes increase by 40%-50%. Hence an underestimation of erosion processes can occur if only automatic model generation is applied. However, this behavior is observed mainly in models corresponding to gullies with steep and narrow sections (2011), but not in those gullies with smoother sections (2013).
Future advances in the methodology can lead to a greater automation of the point (GCPs) extraction and measurement from LiDAR data, as well as the creation of routines to generate DEMs and DODs, to calculate height differences and volumes and to make gully maps in a GIS environment. Further approaches to the automatic edition of models in order to obtain true DTMs, including the using of new LiDAR datasets, will be also addressed. Indeed, these procedures will contribute to apply this methodology to larger extensions in a more efficient way. Moreover, the analysis of a larger extension of gullies with different shapes and properties and, in addition, factor analysis based on classical statistical or modern machine learning techniques will permit a deeper knowledge of the dynamics of erosion processes in the study area.