Next Article in Journal
An Efficient Indoor Wi-Fi Positioning Method Using Virtual Location of AP
Next Article in Special Issue
Application of Hybrid Prediction Methods in Spatial Assessment of Inland Excess Water Hazard
Previous Article in Journal
Selecting Prices Determinants and Including Spatial Effects in Peer-to-Peer Accommodation
Previous Article in Special Issue
Machine Learning for Gully Feature Extraction Based on a Pan-Sharpened Multispectral Image: Multiclass vs. Binary Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Cartographic, Geodetic and Photogrammetric Engineering, University of Jaén, Campus de las Lagunillas s/n, 23071 Jaén, Spain
2
Centre for Advanced Studies in Earth Sciences, Energy and Environment (CEACTEMA), University of Jaén, Campus de las Lagunillas s/n, 23071 Jaén, Spain
3
Department of Geology, University of Jaén, Campus de las Lagunillas s/n, 23071 Jaén, Spain
4
Department of Physics, University of Jaén, Campus de las Lagunillas s/n, 23071 Jaén, Spain
*
Author to whom correspondence should be addressed.
ISPRS Int. J. Geo-Inf. 2020, 9(4), 260; https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9040260
Submission received: 29 February 2020 / Revised: 5 April 2020 / Accepted: 15 April 2020 / Published: 19 April 2020

Abstract

:
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 km2 in olive groves. Several historical flights (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 differences (DoDs) allowed the identification of gullies, the calculation of the affected areas as well as the estimation of height differences and volumes between models. These analyses result in an average depletion of 0.50 m and volume loss of 85000 m3 in the gully area, with some periods (2009–2011 and 2011–2013) showing rates of 10,000–20,000 m3/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 influence of this operation on the erosion calculations, finding that it is not significant except in gully areas with a very steep shape.

Graphical Abstract

1. Introduction

Globally, soil erosion is one of the most disturbing phenomena of environmental degradation [1,2]. The degraded arable area worldwide accounts for 39% of the total land suitable for agriculture [3], a percentage that might be increased significantly in coming decades due to the acceleration of erosion caused by global warming [4]. Soil water erosion is a complex process that must be addressed at multiple spatial and temporal scales. There are studies, carried out on a plot scale, in which the physical processes of laminar and rill erosion predominate [5]. Usually, empirical measurements [6] or different variants of the RUSLE (Revised Unified Soil Loss Equation) model [7,8,9] have been applied in these case studies.
Meanwhile, some authors [10] demonstrate that concentrated flow processes such as gully erosion could explain between 50% and 90% of total erosion at basin scales. Gully erosion is one of the most important soil degradation processes, although it has not received much attention until recent years due to the difficulty of its study. Now it is a fast-growing field [11]. Gullies develop at different spatial and temporal scales. They range from ephemeral gullies with sections of less than 1 m2, lengths of around 10 m and persistence less than an annual cycle, to large permanent gullies with sections of hundreds of m2, several km of development and decades of persistence [10].
Geomatics techniques can support the precise geometric characterization of gully systems since these techniques allow the capture and acquisition of images and ground points coordinates at different spatiotemporal resolutions. Discrete points in gully data surveys can be measured with several methods and instruments such as total station, TS [12,13,14,15,16,17,18], global navigation satellite systems (GNSS) [19,20,21,22,23], or laser altimeters and laser distance meters [16,17,24], although these techniques are generally costly in time and resources. Meanwhile massive point clouds can be captured with Light Detection and Ranging (LiDAR) techniques such as terrestrial laser scanner (TLS) [15,16,21,22,25,26,27], and aerial laser scanner (ALS) [15,22,23,26,28,29,30,31,32,33]. Besides, images can be acquired by means of digital cameras and sensors, in many occasions combined with some point capture techniques. Thus, there are many studies based on images of different resolution and precision, captured with different platforms, from close range [20,34,35] and terrestrial systems [16,23,36,37,38,39,40] to airborne systems using kites [41], unmanned aerial vehicles, UAV [21,23,24,39,42,43,44], conventional aerial platforms [15,19,20,22,28,29,32,33,37,40,43,45,46,47,48,49] or satellites [21,50,51].
Images and orthoimages can be used for surface characterization, so conventional aerial orthophotographs or high-resolution satellite images can suffice for these applications [15,29,33,43,49]. This has allowed to estimate the lengths, widths and densities of the gully systems, analyzing their evolution over time and specifically the displacements or lateral growth of some gullies [15,32,40,43,49,52]. In some cases, the influence of determinant factors on the formation and development of the gullies are analyzed, which allows the application of predictive models [33,43,53,54,55,56].
In addition to planimetric studies, image-based methods permit the generation of digital elevation models (DEMs) for the estimation of gully depths and volumes as well as their variations. Thus, different approaches have been applied, either based on conventional photogrammetric techniques [19,20,24,32,37,45,46,47,48,49] or new computer vision techniques such as Structure from Motion, SfM, and Multi Video Stereo, MVS [21,23,25,38,39]. The availability of models at different epochs leads to multitemporal analysis of gully systems by comparing these models in order to obtain differential models (DEMs of Differences or DoDs). These models can be derived from conventional aerial photogrammetry, which has recently undergone a great development in this field [20,28,32,38,39,40,44]. This approach allows us to survey relatively large areas with high resolutions and also to carry out time-based studies from historical flights available since the mid-20th century [20,28,29,32,47,52]. All of these image-based methods are often combined with point capture techniques such as LiDAR [15,23,29,36], GNSS [20,22,23,45] or conventional surveying [11,15,16,49]. The use of LiDAR leads to the generation of high resolution and reliable DEMs [15,22,25,26,28,32,33,38] and enables us to obtain digital terrain models (DTMs) by classifying and filtering the high-density point clouds.
In addition, Unmanned Aerial Vehicles (UAVs) are well suited for very high resolution and precise surveys in areas of about 0.01 to 100 km2. UAVs are very suitable for intermediate scales between terrestrial techniques (GNSS, close range photogrammetry and TLS) and aerial or space surveys (conventional aerial photogrammetry, LiDAR and VHR satellite imagery) but keeping low or moderate costs and allowing high temporal resolution studies. Most of the current studies using multicopter UAVs are of centimeter resolution [21,23,25,38,39,44].
All these geomatics approaches lead to three-dimensional measurements of gullies’ width, depth and cross-sectional areas from which volumetric estimations are made [14,16,17,23,27,36,50]. They have enabled to determine that gully soil losses at least doubled the rates measured in experimental plots or RUSLE models [49]. However, the most appropriate approach is to generate digital elevation models, both DSMs and DTMs (filtered and/or edited) of high resolution and precision [15,20,21,22,23,24,25,28,29,32,34,38,39,40,41,42,45,46,47,48,51]. Based on these DEMs, morphometric measurements and volumetric estimations can also be performed [15,20,21,28,32,33,38,39,40,41,42,45,46,47,48].
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 (1980–2016), 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 km2 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.

2. Materials and Methods

2.1. Study Area

The study area, with an extension of approximately 7.5 km2, 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 sandstones with subsidiary amounts of carbonates and gypsum as metric tectonic slices. Miocene sediments consist of white marls and marly silts [60].
The area is affected by intense erosion, both laminar and gully, in addition to other superficial processes such as landslides [57,61]. Some sections of the gully area studied, which sometimes affect rural roads and paths (Figure 2c), are shown in Figure 2.
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].

2.2. 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). The image ground sample distance (GSD) of the different flights varies from of 0.27 m in the 1980 flight to 0.30 m in the 1996 and 2001 flights, and 0.45 m in the remaining flights from 2005 to 2016. The 1980, 1996, 2001 and 2005 flights are analog (film-based cameras). The 2005 flight is a Color Infrared (CIR) flight, while the 1977, 1996 and 2001 flights are panchromatic. Original film images had been digitized at different scanning resolutions between 1000 and 1700 ppi, which has implied different image pixel sizes of 0.015–0.025 mm. The remaining flights (2009, 2011, 2013 and 2016) are digital and color-near infrared (RGB-NIR). 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/m2.
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.

2.3. Methodology

The methodology summarized in Figure 3 consists of several phases partially described in previous papers [57,58,70]:
  • Orientation of the historical flights.
  • Generation of the DSMs and orthophotographs.
  • Calculation of DSMs of differences (DoDs).
  • Delimitation of gully areas and DSMs edition of detailed area.
  • Estimation of height differences and volumes between models.

2.3.1. 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/m2 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 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 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.50 and 0.80 m in the GCP and 0.40–0.55 m in the CHK. Meanwhile, Z errors are lower than 0.16 m (except for the 1980 flight where it reaches a value near to 0.3 m) in GCP, but they present a greater variability in the CHK. Thus, the values are of 0.25 m in the 2009 flight; 0.37 m in the 2001 flight; 0.46 m in the 1980 flight; 0.60–0.70 m in the 2005, 2011 and 2016 flights and 0.80–0.85 m in the 1996 and 2013 flights.
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):
Prop YEAR = (RMS LiDAR 2 + RMS PHOTO − YEAR 2) 0.5
Thus, the propagated errors in the horizontal (XY) component range from 0.80 to 0.90 m. Meanwhile, in the vertical (Z) component, they range from 0.50 to 0.65 in the 1980, 2001 and 2008 flights and from 0.80 to 0.95 m in the 1996, 2005, 2011, 2013 and 2016 flights.

2.3.2. 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.

2.3.3. 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]:
Unc. DoD YEAR1-YEAR2 = (Unc. DSM YEAR1 2 + Unc. DSM YEAR2 2) 0.5
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.

2.3.4. 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 (1980–2016).
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.

2.3.5. 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/m3. 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/m3.
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.

3. 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–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.

3.1. 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 km2 and represents 2.25% of the study area of 7.45 km2. 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–2009 and 2013–2016), except in 1996–2001, 2009–2011 and 2011–2013, when the depletion area is about twice the deposition area. 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. Average rates of depletion height differences can be observed of between 0.01 and 0.20 m/year, with the highest values in the periods 1996–2001, 2009–2011 and 2011–2013, being the value for the complete period of 0.02 m/year. Meanwhile, the average rates for deposition height differences are very low for every period (under 0.06 m/year), with a rate for the complete period close to 0. Finally, the rates for balance height differences reach negative values of 0.12 m/year in the period 2009–2011, 0.06 m/year in 2011–2013 and 0.03 m/year in 1996–2001, while the values are practically null in the remaining periods and of 0.014 m/year in the complete period.
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 m3 for depletion and about 17,000 m3 for deposition processes, which results in a negative balance of waste material of 85,000 m3. The rates are near 2800 m3/year for depletion processes and about 500 m3/year for deposition, with a balance (waste) rate of almost 2400 m3/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 m3/year, deposition rates around 8000–10,000 m3/year and waste rates of 10,000–20,000 m3/year. The remaining periods present depletion rates lower than 10,000 m3/year, deposition rates lower than 5000 m3/year with negative balance rates usually lower than 5000 m3/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.

3.2. 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 average rates for depletion height difference are below 0.10 m/year for most periods, except for 2009–2011 and 2011–2013 when the rates reach 0.50 and 0.30 m/year, respectively, 0.05 m/year being the accumulated rate for the complete period. Meanwhile, the rates for deposition show values lower than 0.10 m/year for every period, with the rate for the complete period being close to 0. Finally, the balance height differences present usually rates lower than 0.05 m/year in absolute value, except for 2009–2011 and 2011–2013 when the rates were of −0.40 and −0.30 m/year, respectively, while the balance rate for the complete period was about −0.05 m/year.
The height differences listed in Table 8b for the edited models of the periods 2009–2011 and 2011–2013 show only small differences regarding the values obtained for the models generated automatically, when the break lines are not used. Thus, the rates of depletion height differences reached an average value of 0.54 m/year and 0.33 m/year, respectively, when the resolution is of 2.5 m as well as of 1 m. Meanwhile, the rates for deposition height are below 0.05 m/year, and the balance heights are very similar to those for the depletion. When the resolution is 1 m but break lines are used, the rates for depletion height differences reach 0.75 m/year in the period 2009–2011 and 0.33 m/year in the period 2011–2013. The deposition values maintain values below 0.10 m/year in both periods, thus the balance rates are very similar to the depletion ones in every case.
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–1996) are higher than 30,000 m3 (about a third part of the volumes for the whole area), but the deposition volumes are insignificant. Thus, the volume rate is around 850 m3/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 m3/year for most of the periods, reaching values of between 6000 and 9000 m3/year in the periods 2009–2011 and 2011–2013, while the deposition rates present values equal to or lower than 1000 m3/year in all the periods. Thus, the balance (waste) rates are usually very low but they reach values of around 6000–7000 m3/year in the periods 2009–2011 and 2011–2013. The mass rates present extreme values (higher than 300 t/ha*year) for depletion processes and the balance in the periods 2009–2011 and 2011–2013, and moderate values (lower than 50 t/ha*year) in the remaining periods. The values for deposition processes are low in all periods.
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 to 10,000 m3/year and 6000 m3/year for depletion volumes in the periods 2009–2011 and 2011–2013, while they are lower than 1000 m3/year for deposition processes. Translated to mass rates, the depletion processes reach extreme values higher than 600 t/ha*year in the period 2009–2011 and higher than 300 t/ha*year in the period 2011–2013, while the deposition processes present rates much lower, about 50 t/ha*year. When the resolution increases to 1 m and break lines are considered, the rates of depletion volumes are of 13,500 m3/year in the period 2009–2011 and of 6000 m3/year in 2011–2013, with the deposition rates under 1000 m3/year. Translated to mass rates, the depletion processes present values of 800 t/ha*year in the period 2009–2011 and 350 t/ha*year in the period 2011–2013, with deposition processes being much lower, 37 t/ha*year and 41 t/ha*year, in the same periods.

4. 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.

4.1. 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 Table 3 and Table 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.

4.2. 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 km2 of a total study area of 7.45 km2. Thus, it represents 2.25% of the study area, the gully density being 3.90 km/km2. 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 km2 (Arroyo del Cortijo de la Piedra, Figure 1) are 1.35% and 3.15 km/km2, 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–1996) is about 85,000 m3 (Table 7) which must be evacuated by the drainage network to which the gullies lead. The rates are 2800 m3/year for depletion and 500 m3/year for deposition, which produces a negative balance (waste) near 2400 m3/year. Considering the whole study area (745 ha) and a bulk density of soils of 1.5 t/m3, 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 106 ha order are considered [10,43], rates around and even lower than 1 t/ha*year are found. When the areas are between 102 and 106 ha, the rates are between 1–10 t/ha*year. Finally, areas smaller than 102 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–2009 and 2013–2016, while it is about twice in the periods 1996–2001, 2009–2011 and 2011–2013, especially in 2009–2011 (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 period 2009–2011 (values higher than 0.10 m/year), followed by the period 2011–2013 (0.10–0.06 m/year) and the period 1996–2001 (values lower than 0.05 m/year).
Regarding the volumes, the analysis is similar with three periods, 1996–2001, 2009–2011 and 2011–2013, showing depletion and balance (negative) volumes higher than the other periods, while the deposition volumes are usually low (Table 7). Some periods such as the first one (1980–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 m3/year in the period 2009–2011, around 10,000 m3/year in the period 2011–2013 and 5000 m3/year in the period 1996–2001, while in the remaining periods, the rates are insignificant.
The higher erosion activity, regarding the other ones, in the periods 1996–2001, 2009–2011 and 2011–2013 can be related to rainfall (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, reaching the maximum absolute value of 310 mm in January 2010 (Figure 9c). Thus, the maximum erosion activity that occurs in the period 2009–2011 is related to the rainfalls of the autumn–winter of 2009–2010 and 2010–2011, especially the former; the activity of the period 2011–2013 is related to the rainfalls of the autumn-winter of 2012; and the activity of the period 1996–2001 is related to the rainfalls of the autumn-winter of 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 (1956–2013) 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 Table 6 and Table 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.

4.3. 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 m3, 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 m3/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 m3/year and 5000 m3/year) than the other ones (always under 800 m3/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 m3/year to 13,500 m3/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 (Figure 7e and Figure 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 (Figure 7f and Figure 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.
The analysis of the sections obtained in 2009, 2011 and 2013 allows the estimation of gully depths. Thus, considering the DTMs in this point of the gully system where the erosion processes are very intense, the depth is about 1.3 m in 2009, about 5.5 m in 2011 and almost 6.0 m in 2013.

5. 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 km2, 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 m3, reaching maximum rates of 10,000–20,000 m3/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/m3. 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–2001, 2009–2011 and 2011–2013. These coincide with the periods of highest activity of the erosion processes. The relative lower erosion activity in the rainy period 1996–2001 with respect to 2009–2011 and 2011–2103 leads to further analysis of the influence of other factors. For instance, the presence or the construction of new rural ways and roads, changes in the land use and practices and the mechanisms of gullies’ growth, in addition to the rainfall, should be analyzed when additional data are available. The analysis performed in the detailed area confirms these data and the relationships with triggering factors. The average for the balance difference heights in the complete period is 1.70 m of depletion, with maximum rates of 0.30–0.50 m/year in the periods 2009–2011 and 2011–2013. The balance volumes are around 30,000 m3 for the complete period, reaching values between 5000–8000 m3/year (275–425 t/ha*year) in the periods 2009–2011 and 2011–2013.
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.

Author Contributions

Conceptualization: Tomás Fernández, José Luis Pérez-García, Mario Sánchez-Gómez and Julio Calero; Processing of photogrammetric data: José Luis Pérez-García, José Miguel Gómez-López, Jorge Delgado and Javier Cardenal-Escarcena; Processing and analysis of rainfall data: Joaquín Tovar-Pescador; GIS analysis: Tomás Fernández and José Miguel Gómez-López; Data interpretation: Tomás Fernández and Julio Calero; Supervision: Mario Sánchez-Gómez and Jorge Delgado; Writing, edition and preparation of figures: All authors; Funding acquisition: Jorge Delgado and Julio Calero. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been financed with the following funds: Center for Advanced Studies on Earth Sciences, Energy and Environment of the University of Jaén; Convene “Risks associated to the Road Network of the Jaén Province” of the Deputy of Jaén and the University of Jaén; Group TEP213 of the Andalusian Research, Development and Innovation Plan (PAIDI).

Acknowledgments

We also acknowledge the help of the groups RNM325, RNM127 and TEP220 (Matras) of PAIDI, and the research project SUSTAINOLIVE (PRIMA projects of the UE). We especially acknowledge the work of Margarita Castillo and Francisco Moya in helping with field and lab work. We also thank to Nick Snow for his valuable review of the English manuscript. Finally, we recognize very sincerely the constructive comments of all the reviewers and the editor in order to improve the work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Boardman, J. Soil erosion science: Reflections on the limitations of current approaches. Catena 2006, 68, 73–86. [Google Scholar] [CrossRef]
  2. Borrelli, P.; Robinson, D.; Fleischer, L.R.; Lugato, E.; Balladio, C.; Alewell, C.; Meusburger, K. An assessment of the global impact of 21s century land use change on soil erosion. Nat. Commun. 2013, 8, 1–13. [Google Scholar]
  3. Lal, R. Soil Quality and Agricultural Sustainability. In Soil Quality and Agricultural Sustainability; Lal, R., Ed.; Sleeping Bear Press Inc.: Chelsea, MA, USA, 1998; pp. 3–12. [Google Scholar]
  4. Yang, D.; Kanae, S.; Oki, T.; Koike, T.; Musiake, K. Global potential soil erosion with reference to land use and climate changes. Hydrol. Process. 2003, 17, 2913–2928. [Google Scholar] [CrossRef]
  5. Merrit, W.S.; Letcher, R.A.; Jakeman, A.J. A review of erosion and sediment transport models. Environ. Model. Softw. 2003, 18, 761–799. [Google Scholar] [CrossRef]
  6. Gómez, J.A.; Guzmán, M.G.; Giráldez, J.V.; Fereres, E. The influence of cover crops and tillage on water and sediment yield, and on nutrient, and organic matter losses in an olive orchard on a sandy loam soil. Soil Tillage Res. 2009, 106, 137–144. [Google Scholar] [CrossRef] [Green Version]
  7. Gómez, J.A.; Battany, M.; Renschler, C.S.; Fereres, E. Evaluation the impact of soil management on soil loss in olive orchards. Soil Use Manag. 2003, 19, 127–134. [Google Scholar] [CrossRef]
  8. Vanwalleghem, T.; Infante Amate, J.; González de Molina, M.; Soto Fernández, D.; Gómez, J.A. Quantifying the effect of historical soil management on soil erosion rates in Mediterranean olive orchards. Agric. Ecosyst. Environ. 2011, 142, 341–351. [Google Scholar] [CrossRef]
  9. Parras-Alcántara, L.; Lozano-García, B.; Keestra, S.; Cerdà, A.; Brevik, E. Long-term effects of soil management on ecosystem services and soil loss estimation in olive grove top soils. Sci. Total Environ. 2016, 571, 498–506. [Google Scholar] [CrossRef]
  10. Poesen, J.; Nachtergaele, J.; Verstraeten, G.; Valentin, C. Gully erosion and environmental change: Importance and research needs. Catena 2003, 50, 91–133. [Google Scholar] [CrossRef]
  11. Castillo, C.; Gomez, J.A. A century of gully erosion research: Urgency, complexiy and study approaches. Earth-Sci. Rev. 2016, 160, 300–319. [Google Scholar] [CrossRef]
  12. Lane, S.N.; Richards, K.S.; Chandler, J.H. Morphological estimation of the time-integrated bed load transport rate. Water Resour. Res. 1995, 31, 761–772. [Google Scholar] [CrossRef]
  13. Wheaton, J.M.; Brasington, J.; Darby, S.E.; Sear, D.A. Accounting for uncertainty in DEMs from repeat topographic surveys: Improved sediment budgets. Earth Surf. Process. Landf. 2010, 35, 136–156. [Google Scholar] [CrossRef]
  14. Giménez, R.; Marzolff, I.; Campo, M.A.; Seeger, M.; Ries, J.B.; Casalí, J.; Álvarez-Mozos, J. Accuracy of high-resolution photogrammetric measurements of gullies with contrasting morphology. Earth Surf. Process. Landf. 2009, 34, 1915–1926. [Google Scholar] [CrossRef]
  15. Perroy, R.L.; Bookhagen, B.; Asner, G.P.; Chadwick, O.A. Comparison of gully erosion estimates using airborne and ground-based LiDAR on Santa Cruz Island, California. Geomorphology 2010, 118, 288–300. [Google Scholar] [CrossRef]
  16. Castillo, C.; Pérez, R.; James, M.R.; Quinton, J.N.; Taguas, E.V.; Gomez, J.A. Comparing the Accuracy of Several Field Methods for Measuring Gully Erosion. Soil Sci. Soc. Am. J. 2012, 76, 1319–1332. [Google Scholar] [CrossRef] [Green Version]
  17. Caraballo-Arias, N.A.; Conoscenti, C.; Di Stefano, C.; Ferro, V.; Gómez-Gutiérrez, A. Morphometric and hydraulic geometry assessment of a gully in SW Spain. Geomorphology 2016, 274, 143–151. [Google Scholar] [CrossRef] [Green Version]
  18. Ding, L.; Qin, F.C.; Fang, H.D.; Liu, H.; Zhang, B.; Shu, C.Q.; Deng, Q.C.; Liu, G.C.; Yang, Q.Q. Morphology and controlling factors of the longitudinal profile of gullies in the Yuanmou dry-hot valley. J. Mt Sci. 2017, 14, 674–693. [Google Scholar] [CrossRef]
  19. Brasington, J.; Rumsby, B.T.; McVey, R.A. Monitoring and modelling morphological change in a braided gravel-bedriver using high resolution GPS-based survey. Earth Surf. Proc. Landf. 2000, 25, 973–990. [Google Scholar] [CrossRef]
  20. Rumsby, B.T.; Brasington, J.; Langham, J.A.; McLelland, S.J.; Middleton, R.; Rollinson, G. Monitoring and modelling particle and reach-scale morphological change in gravel-bed rivers: Applications and challenges. Geomorphology 2008, 93, 40–54. [Google Scholar] [CrossRef]
  21. Wang, R.; Zhang, S.; Pu, L.; Yang, J.; Yang, C.; Chen, J.; Guan, C.; Wang, Q.; Chen, D.; Fu, B.; et al. Gully Erosion Mapping and Monitoring at Multiple Scales Based on Multi-Source Remote Sensing Data of the Sancha River Catchment, Northeast China. ISPRS Int. J. Geo-Inf. 2016, 5, 200. [Google Scholar] [CrossRef] [Green Version]
  22. Poulton, P.L.; Caccetta, P.; Wu, X.; Kinsey-Henderson, A.E. Evaluating the Utility of Photogrammetry to Identify and Map Regions at Risk from Gully Erosion; Report to Department of Agriculture and Water Resources, CSIRO Agriculture and Food: Canberra, Australia, 2018; p. 53. [Google Scholar]
  23. Koci, J.; Jarihani, B.; Leon, J.X.; Sidle, R.C.; Wilkinson, S.N.; Bartley, R. Assessment of UAV and Ground-Based Structure from Motion with Multi-View Stereo Photogrammetry in a Gullied Savanna Catchment. ISPRS Int. J. Geo-Inf. 2017, 6, 328. [Google Scholar] [CrossRef] [Green Version]
  24. Liu, K.; Ding, H.; Tang, G.; Na, J.; Huang, X.; Xue, Z.; Yang, X.; Li, F. Detection of Catchment-Scale Gully-Affected Areas Using Unmanned Aerial Vehicle (UAV) on the Chinese Loess Plateau. ISPRS Int. J. Geo-Inf. 2016, 5, 238. [Google Scholar] [CrossRef]
  25. Nadal-Romero, E.; Revuelto, J.; Errea, P.; López-Moreno, J.I. The application of terrestrial laser scanner and SfM photogrammetry in measuring erosion and deposition processes in two opposite slopes in a humid badlands area (central Spanish Pyrenees). Soil 2015, 1, 561–573. [Google Scholar] [CrossRef] [Green Version]
  26. Bartley, R.; Goodwin, N.; Henderson, A.; Hawdon, A.; Tindall, D.; Wilkinson, S.; Baker, B. A Comparison of Tools for Monitoring and Evaluating Channel Change; Queensland Department of Science, Information Technology and Innovation (DSITI), National Environmental Science Programme: Canberra, Australia, 2016; p. 29.
  27. Liu, H.; Qian, F.; Ding, W.; Gómez, J.A. Using 3D scanner to study gully evolution and its hydrological analysis in the deep weathering of southern China. Catena 2019, 183, 104–218. [Google Scholar] [CrossRef]
  28. Lane, S.N.; Westaway, R.M.; Hicks, D.M. Estimation of erosion and deposition volumes in a large, gravel-bed, braided river using synoptic remote sensing. Earth Surf. Process. Landf. 2003, 28, 249–271. [Google Scholar] [CrossRef]
  29. James, L.A.; Singer, M.B.; Ghoshal, S.; Megison, M. Historical Channel Changes in the Lower Yuba and Feather Rivers, California: Long-term Effects of Contrasting River-Management Strategies. In Management and Restoration of Fluvial Systems with Broad Historical Changes and Human Impacts; James, L.A., Rathburn, S.L., Whittecar, G.R., Eds.; Geological Society of America Special Papers: Boulder, CO, USA, 2009; Volume 451, pp. 57–82. [Google Scholar]
  30. James, L.A.; Hodgson, M.E.; Ghoshal, S.; Megison, M. Geomorphic change detection using historic maps and DEM differencing: The temporal dimension of geospatial analysis. Geomorphology 2012, 137, 181–198. [Google Scholar] [CrossRef]
  31. Evans, M.; Lindsey, J. High resolution quantification of gully erosion in upland peatlands at the landscape scale. Earth Surf. Process. Landf. 2010, 35, 876–886. [Google Scholar] [CrossRef]
  32. Ghoshal, S.; James, L.A.; Singer, M.B.; Aalto, R. Channel and Floodplain Change Analysis over a 100-Year Period: Lower Yuba River, California. Remote Sens. 2010, 2, 1797–1825. [Google Scholar] [CrossRef] [Green Version]
  33. Tindall, D.; Marchand, B.; Gilad, U.; Goodwin, N.; Byer, S.; Denham, R. Gully Mapping and Drivers in the Grazing Lands of the Burdekin Catchment; RP66G Summary Report; Remote Sensing Centre, Department of Science, Information Technology, Innovation and the Arts: Brisbane, Australia, 2014; p. 16. [Google Scholar]
  34. Moritani, S.; Yamamoto, T.; Andry, A.H.; Inoue, M.; Kaneuchi, T. Using digital photogrammetry to monitor soil erosion under conditions of simulated rainfall and wind. Aust. J. Soil Res. 2010, 48, 36–42. [Google Scholar] [CrossRef]
  35. Gesch, K.R. Validation and Application of Close-Range Photogrammetry to Quantify Ephemeral Gully Erosion. Master’s Thesis, Iowa State University, Ames, IA, USA, 2015. [Google Scholar]
  36. Sneddon, J.; Williams, G.; Savage, J.K.; Newman, C.T. Erosion of a Gully in Duplex Soils. Results of a Long-term Photogrammetric Monitoring Program. Aust. J. Soil Res. 1988, 26, 401–408. [Google Scholar] [CrossRef]
  37. Lane, S.N. The measurement of river channel morphology using digital photogrammetry. Photogramm. Record 2000, 16, 937–961. [Google Scholar] [CrossRef]
  38. Kaiser, A.; Neugirg, F.; Rock, G.; Müller, C.; Haas, F.; Ries, J.; Schmidt, J. Small-Scale Surface Reconstruction and Volume Calculation of Soil Erosion in Complex Moroccan Gully Morphology Using Structure from Motion. Remote Sens. 2014, 6, 7050–7080. [Google Scholar] [CrossRef] [Green Version]
  39. Stöcker, C.; Eltner, A.; Karrasch, P. Measuring gullies by synergetic application of UAV and close range photogrammetry—A case study from Andalusia, Spain. Catena 2015, 132, 1–11. [Google Scholar] [CrossRef]
  40. Cox, S.E.; Doncaster, D.L.; Godfrey, P.E.; Londe, M.D. Aerial and terrestrial-based monitoring of channel erosion, headcutting, and sinuosity. Environ. Monit. Assess. 2018, 190, 717. [Google Scholar] [CrossRef] [PubMed]
  41. Marzolff, I.; Poesen, J. The potential of 3D gully monitoring with GIS using high-resolution aerial photography and digital photogrammetry system. Geomorphology 2009, 111, 48–60. [Google Scholar] [CrossRef]
  42. D’Oleire-Oltmanns, S.; Marzolff, I.; Peter, K.D.; Ries, J.B. Unmanned Aerial Vehicle (UAV) for Monitoring Soil Erosion in Morocco. Remote Sens. 2012, 4, 3390–3416. [Google Scholar] [CrossRef] [Green Version]
  43. Hughes, A.O.; Prosser, I.P. Gully erosion prediction across a large region: Murray–Darling Basin, Australia. Soil Res. 2012, 50, 267–277. [Google Scholar] [CrossRef]
  44. Yang, S.; Guan, Y.; Zhao, C.; Zhang, C.; Bai, J.; Chen, K. Determining the influence of catchment area on intensity of gully erosion using high-resolution aerial imagery: A 40-year case study from the Loess Plateau, northern China. Geoderma 2019, 347, 90–102. [Google Scholar] [CrossRef]
  45. Brasington, J.; Langham, J.; Rumsby, B. Methodological sensitivity of morphometric estimates of coarse fluvial sediment transport. Geomorphology 2003, 53, 299–316. [Google Scholar] [CrossRef]
  46. Daba, S.; Rieger, W.; Strauss, P. Assessment of gully erosion in eastern Ethiopia using photogrammetric techniques. Catena 2003, 50, 273–291. [Google Scholar] [CrossRef]
  47. Martinez-Casasnovas, J.A. A spatial information technology approach for the mapping and quantification of gully erosion. Catena 2003, 50, 293–308. [Google Scholar] [CrossRef]
  48. Martinez-Casasnovas, J.A.; Ramos, M.C.; Poesen, J. Assessment of sidewall erosion in large gullies using multi-temporal DEMs and logistic regression analysis. Geomorphology 2004, 58, 305–321. [Google Scholar] [CrossRef]
  49. Hayas, A.; Vanwalleghem, T.; Laguna, A.; Peña, A.; Giráldez, J.V. Reconstructing long-term gully dynamics in Mediterranean agricultural areas. Hydrol. Earth Syst. Sci. 2017, 21, 235–249. [Google Scholar] [CrossRef] [Green Version]
  50. Fiorucci, F.; Ardizzone, F.; Rossi, M.; Torri, D. The Use of Stereoscopic Satellite Images to Map Rills and Ephemeral Gullies. Remote Sens. 2015, 7, 14151–14178. [Google Scholar] [CrossRef] [Green Version]
  51. Xu, Q.; Kou, P.; Wang, C.; Yunus, A.P.; Xu, J.; Peng, S.; He, C. Evaluation of gully head retreat and fill rates based on high-resolution satellite images in the loess region of China. Environ. Earth Sci. 2019, 78, 465–480. [Google Scholar] [CrossRef]
  52. Hao, L.; Cruse, R.; Liu, X.; Zhang, X. Effects of Topography and Land Use Change on Gully Development in Typical Mollisol Region of Northeast China. China Geogr. Sci. 2016, 26, 779–788. [Google Scholar]
  53. Conoscenti, C.; Agnesi, V.; Angileri, S.; Cappadonia, C.; Rotigliano, E.; Märker, M. A GIS-based approach for gully erosion susceptibility modelling: A test in Sicily, Italy. Environ Earth Sci. 2013, 70, 1179–1195. [Google Scholar] [CrossRef]
  54. Garosi, Y.; Sheklabadi, M.; Conoscenti, C.; Pourghasemi, H.R.; Van Oost, K. Assessing the performance of GIS-based machine learning models with different accuracy measures for determining susceptibility to gully erosion. Sci. Total Environ. 2019, 664, 1117–1132. [Google Scholar] [CrossRef]
  55. Rahmati, O.; Tahmasebipour, N.; Haghizadeh, A.; Pourghasemi, H.R.; Feizizadeh, B. Evaluating the influence of geo-environmental factors on gully erosion in a semi-arid region of Iran: An integrated framework. Sci. Total Environ. 2017, 579, 913–927. [Google Scholar] [CrossRef]
  56. Dewitte, O.; Daoudi, M.; Bosco, C.; Van Den Eeckhaut, M. Predicting the susceptibility to gully initiation in data-poor regions. Geomorphology 2015, 228, 101–115. [Google Scholar] [CrossRef]
  57. Fernández, T.; Pérez, J.L.; Cardenal, F.J.; Gómez, J.M.; Colomo, C.; Delgado, J. Analysis of landslide evolution affecting olive groves using UAV and photogrammetric techniques. Remote Sens. 2016, 8, 837. [Google Scholar] [CrossRef] [Green Version]
  58. Fernández, T.; Pérez, J.L.; Colomo, C.; Cardenal, J.; Delgado, J.; Palenzuela, J.A.; Irigaray, C.; Chacón, J. Assessment of the Evolution of a Landslide Using Digital Photogrammetry and LiDAR Techniques in the Alpujarras Region (Granada, Southeastern Spain). Geosciences 2017, 7, 32. [Google Scholar] [CrossRef] [Green Version]
  59. Pérez-Valera, F.; Sánchez-Gómez, M.; Pérez-López, A.; Pérez-Valera, L.A. An evaporite-bearing accretionary complex in the northern front of the Betic-Rif orogeny. Tectonics 2017, 6, 1006–1036. [Google Scholar] [CrossRef] [Green Version]
  60. Roldán, F.J.; Lupiani, E.; Jerez, L. Mapa Geológico de España, Escala 1:50.000, Mapa y Memoria Explicativa; Instituto Geológico Nacional: Madrid, Spain, 1988.
  61. Carpena, R.L.; Mellado, I.; Moya, F.; Colomo, C.; Bédmar, P.; Calero, J.; Pérez, A.; Fernández, T.; Sánchez-Gómez, M.; Tovas, J. Análisis de riesgos asociados a las infraestructuras viarias de la Diputación Provincial de Jaén. In Proceedings of the IX Simposio Nacional Sobre Laderas y Taludes Inestables, Santander, Spain, 27–30 June 2017; Volume 1, pp. 335–346. [Google Scholar]
  62. MAPAMA. Estadísticas 2015; Ministerio de Agricultura y Pesca, Alimentación y Medio ambiente: Madrid, Spain, 2016.
  63. FAOSTAT Data. Available online: http://www.fao.org/faostat (accessed on 31 January 2020).
  64. Calero, J.; Aranda, V.; Montejo-Ráez, A.; Martín-García, J.M. A new soil quality index based on morpho-pedological indicators as a sitespecific web service applied to olive groves in the Province of Jaen (South Spain). Comput. Electron. Agric. 2018, 146, 66–76. [Google Scholar] [CrossRef]
  65. MAPAMA. Inventario Nacional de Erosión de Suelos, Provincia de Jaén, Escala 1:250.000; Ministerio de Agricultura, Alimentación y Medio Ambiente: Madrid, Spain, 2006; p. 398.
  66. Calero, J.; Del Moral, J.D.; García-García, F. Evolución morfosedimentaria de la transformación de un embalse en un humedal (Embalse de Doña Aldonza, Alto Guadalquivir). In La Precipitación y los Procesos Erosivos; Proceedings of IV Jornadas de Ingeniería del Agua: Córdoba, Spain, 2015; pp. 145–154. [Google Scholar]
  67. Instituto Geográfico Nacional (IGN), Fototeca Digital. Available online: http://fototeca.cnig.es/ (accessed on 31 January 2020).
  68. Instituto de Estadística y Cartografía de Andalucía (IECA), Fototeca. Available online: http://www.juntadeandalucia.es/institutodeestadisticaycartografia/fototeca/ (accessed on 31 January 2020).
  69. Instituto Geográfico Nacional (IGN), Centro de Descargas CNIG. Available online: http://centrodedescargas.cnig.es/CentroDescargas/index.jsp (accessed on 31 January 2020).
  70. Cardenal, J.; Fernández, T.; Pérez-García, J.L.; Gómez-López, J.M. Measurement of Road Surface Deformation Using Images Captured from UAVs. Remote Sens. 2019, 11, 1507. [Google Scholar] [CrossRef] [Green Version]
  71. Socet Set 5.6; Bae Systems Plc.: London, UK, 2011.
  72. Cardenal, J.; Delgado, J.; Mata, E.; González, A.; Olague, I. Use of historical flight for landslide monitoring. In Proceedings of the Spatial Accuracy 2006, Lisbonne, Portugal, 5–7 July 2006; pp. 129–138. [Google Scholar]
  73. González-Díez, A.; Fernández-Maroto, G.; Doughty, M.W.; Díaz de Terán, J.R.; Bruschi, V.; Cardenal, J.; Pérez, J.L.; Mata, E.; Delgado, J. Development of a methodological approach for the accurate measurement of slope changes due to landslides, using digital photogrammetry. Landslides 2014, 11, 615–628. [Google Scholar] [CrossRef]
  74. Korsgaard, N.; Nuth, C.; Khan, S.; Kjeldsen, K.K.; Bjørk, A.A.; Schomacker, A.; Kjaer, K.H. Digital elevation model and orthophotographs of Greenland based on aerial photographs from 1978–1987. Sci. Data 2016, 3, 160032. [Google Scholar] [CrossRef] [Green Version]
  75. QGIS 3. A Free and Open Source Geographic Information System. 2020. Available online: https://www.qgis.org/en/site/ (accessed on 31 January 2020).
  76. Tarolli, P. High-resolution topography for understanding Earth surface processes: Opportunities and challenges. Geomorphology 2014, 216, 295–312. [Google Scholar] [CrossRef]
  77. Prokešová, R.; Kardoš, M.; Medved’ová, A. Landslide dynamics from high-resolution aerial photographs: A case study from W Carpathians, Slovakia. Geomorphology 2010, 115, 90–101. [Google Scholar] [CrossRef]
  78. DeRose, R.C.; Gomez, B.; Marden, M.; Trustrum, N.A. Gully erosion in Mangatu forest, New Zealand, estimated from digital elevation models. Earth Surf. Process. Landf. 1998, 23, 1045–1053. [Google Scholar] [CrossRef]
Figure 1. Geographical and geological location: (a) geographical location; (b) geological setting; (c) study area: lithological units and gully network delimited from the analysis made in this study (DoD of the complete period) and photointerpretation. Coordinates are in ETRS89-UTM30.
Figure 1. Geographical and geological location: (a) geographical location; (b) geological setting; (c) study area: lithological units and gully network delimited from the analysis made in this study (DoD of the complete period) and photointerpretation. Coordinates are in ETRS89-UTM30.
Ijgi 09 00260 g001
Figure 2. Photographs of different sectors of the gully area: (a), (b) and (c) sections; (c1), (c2) and (c3): rural way affected by erosion processes in January 2016, January 2017 and January 2018, respectively.
Figure 2. Photographs of different sectors of the gully area: (a), (b) and (c) sections; (c1), (c2) and (c3): rural way affected by erosion processes in January 2016, January 2017 and January 2018, respectively.
Ijgi 09 00260 g002
Figure 3. Flow chart of the methodology.
Figure 3. Flow chart of the methodology.
Ijgi 09 00260 g003
Figure 4. Measurement of roof corners for ground control/check points: (a) aerial image of two agricultural sheds (1,2) used for ground control/check points extraction; (b) LiDAR point cloud with the flat roofs of the sheds marked; (c) minimum bounding rectangles adjusted to the LiDAR points located on top of the flat shed roofs. LiDAR point cloud is also represented; (d) minimum bounding rectangles adjusted to the LiDAR points located on top of the flat shed roofs. DSM is represented; (e) corner of the minimum bounding box and field-measured GNSS. The distance represents the error between the corner measured in field (with GNSS) and that computed with the minimum bounding box adjusted to the LiDAR points on top of the roof; (f) GNSS field-measurement of a shed roof corner.
Figure 4. Measurement of roof corners for ground control/check points: (a) aerial image of two agricultural sheds (1,2) used for ground control/check points extraction; (b) LiDAR point cloud with the flat roofs of the sheds marked; (c) minimum bounding rectangles adjusted to the LiDAR points located on top of the flat shed roofs. LiDAR point cloud is also represented; (d) minimum bounding rectangles adjusted to the LiDAR points located on top of the flat shed roofs. DSM is represented; (e) corner of the minimum bounding box and field-measured GNSS. The distance represents the error between the corner measured in field (with GNSS) and that computed with the minimum bounding box adjusted to the LiDAR points on top of the roof; (f) GNSS field-measurement of a shed roof corner.
Ijgi 09 00260 g004
Figure 5. DSMs and orthophotographs: (a) DSM derived from LiDAR data (2014); (b) orthophotograph of 1980; (c) orthophotograph of 2016. Coordinates are in ETRS89-UTM30.
Figure 5. DSMs and orthophotographs: (a) DSM derived from LiDAR data (2014); (b) orthophotograph of 1980; (c) orthophotograph of 2016. Coordinates are in ETRS89-UTM30.
Ijgi 09 00260 g005
Figure 6. DSMs of differences (DoDs) of the whole area: (a) period 1980–1996; (b) period 1996–2001; (c) period 2001–2005; (d) period 2005–2009; (e) period 2009–2011. Coordinates are in ETRS89/UTM30; DSMs of differences (DoDs) of the whole area: (f) period 2011–2013; (g) period 2013–2016; (h) complete period (1980–2016). Coordinates are in ETRS89/UTM30.
Figure 6. DSMs of differences (DoDs) of the whole area: (a) period 1980–1996; (b) period 1996–2001; (c) period 2001–2005; (d) period 2005–2009; (e) period 2009–2011. Coordinates are in ETRS89/UTM30; DSMs of differences (DoDs) of the whole area: (f) period 2011–2013; (g) period 2013–2016; (h) complete period (1980–2016). Coordinates are in ETRS89/UTM30.
Ijgi 09 00260 g006aIjgi 09 00260 g006b
Figure 7. DSMs of differences (DoDs) of the detailed area: (a) period 1980–1996; (b) period 1996–2001; (c) period 2001–2005; (d) period 2005–2009; (e) period 2009–2011; (f) period 2011–2013; (g) period 2013–2016; (h) complete period (1980–2016).
Figure 7. DSMs of differences (DoDs) of the detailed area: (a) period 1980–1996; (b) period 1996–2001; (c) period 2001–2005; (d) period 2005–2009; (e) period 2009–2011; (f) period 2011–2013; (g) period 2013–2016; (h) complete period (1980–2016).
Ijgi 09 00260 g007
Figure 8. DTMs of differences of the detailed area: (a) period 2009–2011; (b) period 2011–2013.
Figure 8. DTMs of differences of the detailed area: (a) period 2009–2011; (b) period 2011–2013.
Ijgi 09 00260 g008
Figure 9. Relationships between erosion processes and rainfall: (a) precipitations in 24 hours (1 day); (b) precipitations in 7 days (1 week); (c) precipitations in 1 month.
Figure 9. Relationships between erosion processes and rainfall: (a) precipitations in 24 hours (1 day); (b) precipitations in 7 days (1 week); (c) precipitations in 1 month.
Ijgi 09 00260 g009
Figure 10. Orthophotographs and DSMs/DTMs sections: (a) orthophotograph of 2009; (b) orthophotograph of 2011; (c) orthophotograph of 2013; (d) sections of 2009; (e) sections of 2011; (f) sections of 2013.
Figure 10. Orthophotographs and DSMs/DTMs sections: (a) orthophotograph of 2009; (b) orthophotograph of 2011; (c) orthophotograph of 2013; (d) sections of 2009; (e) sections of 2011; (f) sections of 2013.
Ijgi 09 00260 g010
Table 1. Properties of the images and LiDAR data sets.
Table 1. Properties of the images and LiDAR data sets.
Camera and image properties
DateBandsFormatScaleCameraDigitalization resolution (mm)GSD (m)
1980PanchromaticFilm1:18,000Wild RC100.0250.27
1996PanchromaticFilm1:20,000Wild RC100.0200.30
2001PanchromaticFilm1:20,000Leica RC300.0150.30
2005CIRFilm1:30,000Leica RC300.0150.45
2009RGB-NIRDigital1:30,000Z/I DMC120---0.45
2011RGB-NIRDigital1:30,000Z/I DMC120---0.45
2013RGB-NIRDigital1:30,000Vexcel UCXp---0.45
2016RGB-NIRDigital1:30,000Vexcel UCXp---0.45
LiDAR data set
DateSystemFlying Height m (above sea level, asl)points/m2
2014Leica ALS6027000.5–1
Table 2. Error of LiDAR point cloud regarding GNSS-based check points 1.
Table 2. Error of LiDAR point cloud regarding GNSS-based check points 1.
StatisticXYXYZ
Minimum−0.655−0.8080.256−0.067
Maximum0.6861.4421.5211.915
Mean0.0900.1790.6490.236
Standard deviation0.4050.5910.3370.403
RMS0.4050.6040.7270.459
1 Errors are expressed in m.
Table 3. Errors of orientation process at control (GCP) and check (CHK) points 1.
Table 3. Errors of orientation process at control (GCP) and check (CHK) points 1.
DateImage NumberGCP/CHK NumberTie Points NumberRMS PixelRMS GCPRMS CHKRMS Prop.
XYZXYZXYZ
198037/11970.1870.7290.2710.3970.4600.8280.650
199638/131050.2660.5070.0240.5450.8090.9080.930
2001310/15970.2560.7680.1230.3960.3670.8280.588
2005312/131450.4220.6910.1650.5080.7060.8870.842
2009314/161430.1480.5780.1590.4730.2470.8680.521
2011313/121470.1470.5750.1670.5330.6640.9010.807
2013316/181270.1290.6750.0920.3980.8320.8290.950
2016318/251230.1660.5450.1820.3610.6320.8120.781
1 RMS errors are expressed in m.
Table 4. Z errors and vertical uncertainties of the models (DSMs and DoDs) 1.
Table 4. Z errors and vertical uncertainties of the models (DSMs and DoDs) 1.
DSM DateDSM UncertaintyIntervalDoD Uncertainty
19800.650
19960.9301980–19961.135
20010.5881996–20011.100
20050.8422001–20051.027
20090.5212005–20090.990
20110.8072009–20110.961
20130.9502011–20131.247
20160.7812013–20161.230
1 Errors and uncertainties are expressed in m.
Table 5. Areas affected by the gullies: depletion and deposition 1.
Table 5. Areas affected by the gullies: depletion and deposition 1.
PeriodWith the Uncertainty AreaWithout the Uncertainty Area
Depletion DepositionUncertaintyDepletionDeposition
Area%Area%Area%Area%Area%
1980–19961879411.222205013.1612669475.628343849.808410650.20
1996–20012980017.79106256.3412711975.8710601963.286152536.72
2001–200587195.2086565.1715016989.636633139.5910121360.41
2005–2009107316.41120817.2114473186.387966947.558787552.45
2009–20113465020.68127317.6012016371.7210453162.396301337.61
2011–20131943811.6075884.5314051983.879084454.227670045.78
2013–201657063.4191755.4815266391.128781952.427972547.58
1980–20165148130.73116506.9510441362.3210616963.376137536.63
1 All the areas are m2.
Table 6. Height differences and rates in the gullies of the whole study area 1,2.
Table 6. Height differences and rates in the gullies of the whole study area 1,2.
PeriodDepletionDepositionBalance
MinAverageRateMaxAverageRateAverageRate
1980–1996−4.182−0.165−0.0103.3400.1740.0110.0090.001
1996–2001−4.236−0.250−0.0504.0220.1020.020−0.148−0.030
2001–2005−4.880−0.075−0.0194.5700.0970.0240.0210.005
2005–2009−4.897−0.105−0.0265.1250.1140.0290.0090.002
2009–2011−5.910−0.356−0.17817.9840.1230.061−0.233−0.117
2011–2013−14.883−0.206−0.1033.8720.0800.040−0.125−0.063
2013–2016−3.615−0.073−0.0244.3670.0930.0310.0200.007
1980–2016−6.824−0.611−0.0172.9900.1000.003−0.510−0.014
1 Average, minimum and maximum height differences are in m. 2 Rates are in m/year.
Table 7. Volumes and rates in the gullies of the whole study area 1,2,3.
Table 7. Volumes and rates in the gullies of the whole study area 1,2,3.
PeriodDepletion DepositionBalance
Vol.Vol. RMass RVol.Vol. RMass RVol.Vol. RMass R
1980–1996−27,644−1728−3.4829,09318183.661449910.18
1996–2001−41,874−8375−16.8617,11034226.89−24,765−4953−9.97
2001–2005−12,641−3160−6.3616,17040438.1435298821.78
2005–2009−17,593−4398−8.8619,11547799.6215223800.77
2009–2011−59,614−29807−60.0120,53310,26720.67−39,081−19,540−39.34
2011–2013−34,462−17231−34.6913,464673213.55−20,997−10,499−21.14
2013–2016−12,164−4055−8.1615,567518910.45340411352.28
1980–2016−102,333−2843−5.7216,8034670.94−85,530−2376−4.78
1 Volumes are in m3; 2 Volume rates are in m3/year; 3 Mass rates are t/ha*year.
Table 8. (a) Height differences and rates in the gullies of the detailed area for DSMs 1,2. (b) Height differences and rates in the gullies of the detailed area for edited DTMs 1,2.
(a)
(a)
PeriodDepletionDepositionBalance
MinAverageRateMaxAverageRateAverageRate
1980–1996−2.2610.0070.0002.3210.1650.0100.1720.011
1996–2001−3.281−0.233−0.0471.7250.0370.007−0.196−0.039
2001–2005−3.497−0.162−0.0403.0390.0490.012−0.113−0.028
2005–2009−3.868−0.098−0.0253.9740.0700.018−0.028−0.007
2009–2011−5.910−0.967−0.4843.9610.1410.070−0.827−0.413
2011–2013−5.720−0.641−0.3203.0270.0200.010−0.620−0.310
2013–2016−3.355−0.116−0.0393.8570.2120.0710.0960.032
1980–2016−6.462−1.691−0.0472.0530.0160.000−1.675−0.047
1 Average, minimum and maximum height differences are in m. 2 Rates are in m/year.
(b)
(b)
PeriodDepletionDepositionBalance
MinAverageRateMaxAverageRateAverageRate
2009–2011 3−5.910−0.967−0.4843.9610.1410.070−0.827−0.413
2011–2013−5.720−0.641−0.3203.0270.0200.010−0.620−0.310
2009–2011 4−6.442−1.075−0.5385.4810.0980.049−0.978−0.489
2011–2013−5.377−0.659−0.3303.3880.0750.038−0.584−0.292
2009–2011 5−6.365−1.077−0.5395.4810.0940.047−0.983−0.492
2011–2013−5.856−0.665−0.3323.3410.0760.038−0.589−0.295
2009–2011 6−6.612−1.495−0.7473.4270.0690.034−1.426−0.713
2011–2013−6.099−0.662−0.3313.5080.0760.038−0.586−0.293
1 Average, minimum and maximum height differences are in m. 2 Rates are in m/year. 3 DSM with resolution of 2.5 m. 4 Edited DTM without break lines and resolution of 2.5 m. 5 Edited DTM without break lines and resolution of 1 m. 6 Edited DTM with break lines and resolution of 1 m.
Table 9. (a) Depletion and deposition volumes and rates in detailed area for DSMs 1,2,3. (b) Depletion and deposition volumes and rates in detailed area for edited DTMs 1,2,3.
(a)
(a)
PeriodDepletion DepositionBalance
Vol.Vol. RMass RVol.Vol. RMass RVol.Vol. RMass R
1980–199611870.44299818711.24311619511.69
1996–2001−4223−845−50.686711348.06−3552−710−42.62
2001–2005−2938−735−44.0889622413.43−2043−511−30.64
2005–2009−1785−446−26.78127131819.07−514−128−7.71
2009–2011−17,557−8778−526.702552127676.55−15,005−7503−450.15
2011–2013−11,626−5813−348.7836718311.00−11,259−5630−337.78
2013–2016−2102−701−42.053841128076.81173857934.76
1980–2016−30,685−852−51.1428280.47−30,403−845−50.67
1 Volumes are in m3. 2 Volume rates are in m3/year. 3 Mass rates are t/ha*year.
(b)
(b)
PeriodDepletion DepositionBalance
Vol.Vol. RMass RVol.Vol. RMass RVol.Vol. RMass R
2009–2011 4−17,557−8778−526.702552127676.55−15,005−7503−450.15
2011–2013−11,626−5813−348.7836718311.00−11,259−5630−337.78
2009–2011 5−19,528−9764−585.85177488753.22−17,754−8877−532.63
2011–2013−11,972−5986−359.17136768341.00−10,606−5303−318.17
2009–2011 6−19,454−9727−583.61169985050.98−17,754−8877−532.63
2011–2013−12,007−6003−360.20137068541.09−10,637−5319−319.11
2009–2011 7−26,987−13,493−809.61124462237.32−25,743−12872−772.29
2011–2013−11,952−5976−358.56138069041.39−10,573−5286−317.18
1 Volumes are in m3. 2 Volume rates are in m3/year. 3 Mass rates are t/ha*year. 4 DSM with resolution of 2.5 m. 5 Edited DTM without break lines and resolution of 2.5 m. 6 Edited DTM without break lines and resolution of 1 m. 7 Edited DTM with break lines and resolution of 1 m.

Share and Cite

MDPI and ACS Style

Fernández, T.; Pérez-García, J.L.; Gómez-López, J.M.; Cardenal, J.; Calero, J.; Sánchez-Gómez, M.; Delgado, J.; Tovar-Pescador, J. Multitemporal Analysis of Gully Erosion in Olive Groves by Means of Digital Elevation Models Obtained with Aerial Photogrammetric and LiDAR Data. ISPRS Int. J. Geo-Inf. 2020, 9, 260. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9040260

AMA Style

Fernández T, Pérez-García JL, Gómez-López JM, Cardenal J, Calero J, Sánchez-Gómez M, Delgado J, Tovar-Pescador J. Multitemporal Analysis of Gully Erosion in Olive Groves by Means of Digital Elevation Models Obtained with Aerial Photogrammetric and LiDAR Data. ISPRS International Journal of Geo-Information. 2020; 9(4):260. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9040260

Chicago/Turabian Style

Fernández, Tomás, José Luis Pérez-García, José Miguel Gómez-López, Javier Cardenal, Julio Calero, Mario Sánchez-Gómez, Jorge Delgado, and Joaquín Tovar-Pescador. 2020. "Multitemporal Analysis of Gully Erosion in Olive Groves by Means of Digital Elevation Models Obtained with Aerial Photogrammetric and LiDAR Data" ISPRS International Journal of Geo-Information 9, no. 4: 260. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9040260

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop