New Method for Automated Quantification of Vertical Spatio-Temporal Changes within Gully Cross-Sections Based on Very-High-Resolution Models

Abstract

Gully erosion is one of the most prominent natural denudation processes of the Mediterranean. It causes significant soil degradation and sediment yield. Most traditional field methods for measurement of erosion-induced spatio-temporal changes are time and labor consuming, while their accuracy and precision are highly influenced by various factors. The main research question of this study was how the measurement approach of traditional field sampling methods can be automated and upgraded, while satisfying the required measurement accuracy. The VERTICAL method was developed as a fully automated raster-based method for detection and quantification of vertical spatio-temporal changes within a large number of gully cross-sections (GCs). The developed method was tested on the example of gully Santiš, located at Pag Island, Croatia. Repeat unmanned aerial vehicle (UAV) photogrammetry was used, as a cost-effective and practical method for the creation of very-high-resolution (VHR) digital surface models (DSMs) of the chosen gully site. A repeat aerophotogrammetric system (RAPS) was successfully assembled and integrated into one functional operating system. RAPS was successfully applied for derivation of interval (the two-year research period) DSMs (1.9 cm/pix) of gully Santiš with the accuracy of ±5 cm. VERTICAL generated and measured 2379 GCs, along the 110 m long thalweg of gully Santiš, within which 749 052 height points were sampled in total. VERTICAL proved to be a fast and reliable method for automated detection and calculation of spatio-temporal changes in a large number of GCs, which solved some significant shortcomings of traditional field methods. The versatility and adaptability of VERTICAL allow its application for other, similar scientific purposes, where multitemporal accurate measurement of spatio-temporal changes in GCs is required (e.g., river material dynamics, ice mass dynamics, tufa sedimentation and erosion).

1. Introduction

Gully erosion is one of the most prominent natural denudation processes of the Mediterranean [1,2,3,4,5,6,7,8,9], which causes significant soil degradation and sediment yield [10,11,12,13,14,15]. Active gullies tend to grow as long as predisposing factors, such as lithology, vegetation cover, land use, terrain attributes, and climatic factors, sustain erosion processes and soil removal [16,17,18]. Research about gully erosion is mostly aimed at quantification of different aspects of spatio-temporal changes in gully geometry, for instance, gully headcut retreat rate [4,15,19,20,21] or changes within gully cross-sections (GCs) [22,23,24,25]. Within this research, we concentrated on the accurate quantification of spatio-temporal changes in GCs as indicators of overall gully evolution. Spatio-temporal changes within GCs include all changes in cross-section geometry caused by erosion and/or accumulation processes that have occurred at the chosen study area, within a certain time period (e.g., day, month, year).

The choice of the most appropriate measurement method for the detection and quantification of spatio-temporal changes depends primarily on the aim and spatial scope of the research, as well as on the desirable and achievable accuracy of the measurements [26,27]. A variety of traditional sampling methods for field study of changes in gully cross-sections exist (e.g., pole [27], tape and ruler [26,28,29], total station [27,30], profilometer [22,24,25,26,27,31,32,33,34]). These methods differ in terms of precision and accuracy, time and cost-effectiveness, spatial coverage, complexity, and required expertise of researchers [26,27]. Characteristics, as well as detailed accuracy assessment of various field sampling methods, were given by [26]. While field sampling methods can have satisfactory accuracy, their main lack is that they are, due to the high time ineffectiveness and difficulty of field sampling [26,35], spatially limited to measurement of changes in a small number of GCs per surveyed gully (

Table 1. Review of the recent studies focused on the measurement of GCs.

ID	Case Study (Authors)	Measurement Method (Method Type)	${\mathit{N}}_{\mathit{o}}\mathbf{of}\mathbf{Sampled}\mathbf{Gullies}$	${\mathit{N}}_{\mathit{o}}$ of Sampled GCs Per Gully	Measured GCs Application
1	Sapphire Mountains, Montana, USA [77]	Tape (direct method)	6	5	Volume; SF; D; W; W/D ratio
2	Umbulo catchment, Ethiopia [28]	Tape (direct method)	15	1	Volume; W/D ratio; D; W
3	Bardenas Reales, Navarre, Spain [22]	Laser profiliometer (indirect method); Aerial photogrammetry (indirect method)	5	4–6	D; TW; BW; CSA; W/D
4	Avon-Richardson Catchment, Victoria, Australia [36]	Aerial Photo Interpretation (indirect method)	89	1	W; D; CSA
5	Ethiopia [71]	Tape (direct method)	811	1	W; D; TW; BW; CSA; W/D; SF; volume
6	Pravara River, Western India [78]	Profilometer; Erosion pins (direct methods)	5	1	CSA; volume
7	Belgium; Ethiopia [66]	Ground photogrammetry (indirect method); Tape (direct method)	4	1–2	W; D; TW; BW; CSA
8	Yuanmou Dry-Hot Valley, Yunnan Province, China [24]	Laser distance meter (indirect method)	152	3	26 different morphological GCs parameters
9	Extremadura, SW Spain [30]	Laser total station (indirect method)	1	28	W; D; CSA; volume
10	Loess Plateau, China [67]	Terestric laser scanning (indirect method)	44	2–3	D; TW; BW; CSA; W/D
11	Yuanmou Dry-Hot Valley, Yunnan Province, China [25]	Laser distance meter (indirect method)	152	3	26 different morphological GCs parameters
12	Cordoba, Spain [68]	FreeXSapp (indirect method); tape (direct method)	1	10	W; D; CSA; volume
13	Loess Plateau, China [69]	Terestric laser scanning (indirect method)	31	6	D; W/D
14	New Brunswick, Canada [29]	Ground photogrammetry (indirect method); Tape (direct method)	1	10	CSA;

). An additional disadvantage is that most traditional sampling methods are direct measurement techniques, that require direct physical contact between measurement probe and soil surface. Such direct contact can cause the occurrence of soil surface compaction, resulting in measurement overestimation [29,36,37].

In recent years, different remote sensing methods have emerged as cost and time-effective indirect techniques for accurate data collection [20,27,38,39,40,41] that allow for the creation of very-high-resolution (VHR) digital surface models (DSMs). The use of Structure-from-Motion (SfM) algorithms and unmanned aerial vehicles (UAV) photogrammetry has been recently particularly popular for the creation of VHR DSMs and as such it has wide application in various geomorphological researches [42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63]. In particular, the application of UAV photogrammetry and SfM algorithms have found wide application in soil erosion related research [51,52,53,54,55,56,57,58,59,60,64]. VHR DSMs provide the perfect basis for a pixel-orientated analysis of morphological changes and quantification of spatio-temporal changes at a sub-decimeter scale [40,65]. Although interval VHR DSMs are mostly used for detection and quantification of areal and/or volumetric spatio-temporal changes, lately the interest in the use of DSMs for accurate measurement and evaluation of GCs characteristics has been growing [29,66,67,68,69,70]. From a scientific perspective, research of spatio-temporal changes in GCs is still relevant and important, as cross-sections reveal valuable geomorphological information that cannot be observed from volumetric measurements derived from continuous surfaces [24,68]. While volumetric spatio-temporal changes indicate overall erosion or accumulation of soil material, measurement of spatio-temporal changes within GCs can provide insight into: cross-sectional metrics (e.g., cross-section width (W); gully top width (TW); bottom width (BW); gully depth (D); gully cross-sectional area (CSA); width–depth ratio (W/D ration); shape factor (SF)) [24,70,71], erosion rate [24], gully evolutionary stage [72] (72 in 68), soil resistance [73,74] (73,74 in 68), dominance of certain types of erosive processes [75,76] (75,76 in 68), etc. A detailed description of more than 20 different geometric and morphological parameters that can be extracted for every measured GCs is given in [24].

A review of various case studies in which direct and/or indirect measurement methods were applied for sampling and measurement of GCs is given in Table 1. It is evident that measurement methods and $N_o$ of sampled GCs per gully vary considerably.

The main research question of this research was how to automate and improve the measurement process of gully cross-sections (GCs), while achieving high measurement accuracy. Therefore, the main research objective (1) was to develop the pixel-based (interval VHR DSMs) methodological approach that will allow automated measurement of vertical spatio-temporal changes in a large number of GCs. A conceptual basis for the development of the VERTICAL method was based on the principle of profilometer, as one of the most commonly used devices for manual field measurement of vertical spatio-temporal changes within chosen cross-sections [22,24,25,26,27,31,32,33,34,79]. A profilometer is a simple device that allows interval field measurement of depth at point samples distributed within identical distance intervals along the chosen cross-section [26,80,81]. Older versions of profilometer devices have aluminum pins (Figure 1) that allow mechanical measurement of depth by lowering the pin until it reaches the surface [79,81,82]. Newer versions are equipped with optical devices, e.g., laser distance gauge [22,34].

Accurate VHR DSMs of the chosen gully are essential for accurate measurement of spatio-temporal changes in GCs by the developed VERTICAL method. Therefore, repeat UAV photogrammetry was used, as a cost-effective and practical method [83,84,85] for the derivation of interval (within a two-year research period) DSMs. Although numerous ready-to-fly UAVs are most often used for different geomorphic surveys, they have limited flight capabilities and camera characteristics [85,86,87,88]. Therefore, a repeat aerophotogrammetric system (RAPS) was assembled and functionally integrated into a high-end UAV system that was applied for two interval aerial surveys.

The second research objective (2) of our study was to determine and interpret the intensity of spatio-temporal changes at GCs of the chosen gully within the 2-year research period. For that purpose, the VERTICAL method was developed and applied on interval VHR DSMs, of the chosen gully at Pag Island, Croatia. With several hundred recorded active gullies [89,90] Pag Island can be considered as a very suitable location for gully erosion research. Overall, studies considering soil erosion in Croatia are scarce [89,90,91,92,93], especially in regard to the application of modern geospatial technologies and advanced research methods for temporal monitoring of spatio-temporal changes. Therefore, thorough longitudinal research of gully erosion is of crucial importance for a better understanding of the intensity of overall soil erosion dynamics.

2. Study Area

Pag Island (284 km²) is the largest island in the Northern Dalmatia archipelago (Croatia) (Figure 2A). Structurally Pag is characterized by alteration of several folds of Dinaric direction (NW-SE) [94]. Prevailing parent material is composed of Upper Cretaceous and Eocene limestones, occasionally covered with Dalmatian Flysch and scarce Kalkocambisol sediments [95]. Deposition of diluvial sediments occurred during the Pleistocene, while most recent Holocene layers within Pag Island are represented by alluvium and organogenic-swamp sediments [95] (Figure 2B).

The island has a mild temperate climate with hot summers (Cfa) [96]. In the period between 1981 and 2011 mean annual amount of precipitation was 977.5 mm, with pronounced seasonal distribution during the autumn and early winter. The rainiest months are November (124.1 mm) and September (115.6 mm) [97]. The average annual temperature is 15.5 °C [97]. The absence of protective vegetation cover leads to exposure of surface soil materials to the influence of various exogenous processes and anthropogenic influences. Among the exogenous processes, climatic characteristics have a great significance for relief formation [89]. A highly developed karstic landscape is characterized by shallow and scarce soil cover. There is no permanent surface runoff, although periodic runoff frequently occurs due to intense precipitation. Such intense periodic surface runoff is important for the formation of gullies.

The chosen study site is gully Santiš (1163 m²) (Figure 2D), located within the south-eastern part of Pag Island (Figure 2B). Gully Santiš is located at the very end of a larger (0.18 km²) drainage basin (Figure 2C), where it is formed in relatively deep accumulated soil sediments of Kalkocambisol, a soil type formed by the dissolution of limestone and/or dolomite [98]. Steep, almost vertical headcut, about 15 m wide, forms the initial part of the gully Santiš, where most intensive and complex gully erosion processes had been observed during the field research. After the headcut, the gully narrows towards the direction of the main channel, where parental material is less homogenous and where fewer traces of active erosion can be found. Around 80 m from the headcut, at the contact zone with the Adriatic Sea, a gully forms a small pebble beach. Gully Santiš was chosen for this research as a simple, unbranched, and relatively short gully, with very scarce vegetation cover (e.g., short grass). Recent intensive active gully erosion traces are visible and there are no obvious anthropogenic influences that could disturb natural erosion processes.

3. Materials and Methods

3.1. Assembling and Functional Integration of Repeat Aerophotogrametric System (RAPS)

Repeat aerophotogrammetric system (RAPS) successfully integrates different high-grade components (Figure 3). DJI Matrice 600 PRO (Figure 3A) is a professional UAV, that represents a basis of developed RAPS, chosen for its advanced flight capabilities, the ability to carry a payload up to 5.9 kg, and compatibility with various gimbals and add-ons (e.g., professional-grade DSLR cameras, multispectral and thermal cameras, aeroLiDAR solutions) [99]. Chosen UAV was upgraded with Gremsy T3 gimbal (Figure 3B), which was chosen due to its advanced camera stabilization, 360° horizontal and 90° vertical rotation, and compatibility with various professional cameras and different lenses [100]. During the functional integration of RAPS, it was not possible to connect the DSLR camera to the UAV’s GPS receiver, and thus collected images lacked the information needed for georeferencing. Therefore, Reach M+ GNSS module for UAV mapping (Figure 3C), which allows accurate determination of xyz coordinates of every collected aerial image [101] was added into the RAPS setup. Professional Sony Alpha A7RII (42 MP) DSLR camera equipped with 20 mm lens (Figure 3D), which allowed the acquisition of VHR aerial imagery, was mounted on Gremsy T3 gimbal. Although all components of RAPS are important in the derivation of VHR DSMs, the professional-grade Sony Alpha A7RII (42 MP) DSLR camera had the highest impact on the quality of created models. Detailed characteristics of RAPS are given in

Table 2. Component (A), parameters (B) and characteristics (C) of repeat aerophotogrammetric system (RAPS).

#	Component (A)	Parameter (B)	Value (C)
1.	DJI Matrice 600 PRO	Flight time (min)	16–32 min
2.		Max takeoff weight (kg)	15.5
3.		Max wind resistance (m/s)	8
4.		Max height above sea level (m)	2500
5.		Max transmission distance (m)	5000
6.	Sony Alpha A7RII	Sensor size	861.6 mm² (35.90mm × 24.00mm)
7.		Camera weight (kg)	0.64
8.		Aperture	f/3.5–f/22
9.		Sensor (px)	7952 × 5304
10.		ISO	100–25600
11.		Shutter Speed	1/8000–30 sec
12.		Focal Length (mm)	28–70

.

3.2. Field Data Acquisition

Repeat UAV photogrammetry was used as a cost-effective and practical method for interval acquisition of aerial imagery [83,84,85]. Two aerial surveys were conducted with developed RAPS within the interval of nearly 2 years. The initial survey was carried out on 26 May 2017 (Survey A), while the final survey (Survey B) was conducted on 11 March 2019. UAVs flight missions for both surveys were planned and conducted through the Universal Ground Control Software (UgCS) commercial UAV flight planning application. It allows user to define various flight mission parameters (e.g., flight profiles, flight height, percentage of side and forward overlap) [102] (Figure 4). Identical UAVs flight mission was used for both aerial surveys. Double-grid flight profiles were chosen and 85% forward and side overlap was set. In total 233 aerial images were collected per each survey.

In order to reduce systematic error of DSMs created from acquired data, collected aerial images were georeferenced accordingly to the seven fixed ground control points (GCPs) scattered throughout the study area (Figure 2D). The number of GCPs was determined following examples of good practice [45,103]. Along with seven GCPs, additional seven fixed checkpoints (CPs) were marked in the field as a quality measure [51]. The number of GCPs and CPs is also conditioned by land cover type in the wider area of the case study. Namely, shallow-brown soil predominates in this area. On such, extremely dynamical and erosion susceptible surface, it was not possible to mark and construct a lot of permanent, fixed geodetic points, i.e., to measure the GCPs and CPs using the proposed methodology. Therefore, the GCPs were limited to predominantly rocky areas within the study area. All GCPs and CPs were identical for both surveys, thus ensuring consistency of created models. GCP and CP locations were marked before the initial and final aerial survey by 15 × 15 cm X signs (Figure S1). Additionally, a Bosch hammer was used to drill a small hole for the rod of the RTK-GPS rover at the center of the marked cross. Red weather-resistant spray paint was used for marking X signs, which stands out well from surrounding light carbonate rocks (Figure 2D). Spray paint applied on the carbonate rock base remained visible after two years, but it was enhanced with new spray before the final aerial survey. Precise coordinates (XYZ) of every selected GCP and CP were collected before every survey with Stonex S10 RTK GPS, with 0,8 cm horizontal and 1,5 cm vertical precision [104], where every point was observed for a minimum of 5 minutes (300 epochs). These data were used for the validation of RTK GPS precision. Validation of RTK GPS coordinates has shown that the precision of collected coordinates slightly deviates from the stated factory precision (horizontal precision = 1.2 cm; vertical precision = 1.8 cm). Such deviation is probably caused by the remoteness of the study location, in relation to the three closest base stations of the Croatian positioning system (CROPOS) used for the correction of the collected XYZ coordinates. Identical GCPs and CPs were used in both aerial surveys. Acquired aerial photos were later used for photogrammetric processing and the derivation of VHR DSMs.

Along with aerial surveys, Hobo Onset RG3-M data logger, intended to record the amount and intensity of precipitation was installed next to the study area. Unfortunately, due to the weather-related mechanical failure, collected precipitation data are reliable only for the sixth month period, between 30 April 2017, and 1 November 2017.

3.3. Aerial Data Processing and VHR DSM Creation

Imagery acquired through aerial surveys were processed in Agisoft Metashape 1.5.1. software, which allows photogrammetric processing of aerial images and creation of VHR DSMs [105]. Agisoft Metashape is image-based 3D modelling software, with an implemented SfM algorithm and multi-view 3D reconstruction technology that enables the creation of precise point clouds and 3D structures from 2D image sets [106,107]. Processing steps and parameters, as well as user-defined options/values used for the derivation of two-interval VHR DSMs, are given in

Table 3. Processing steps (A), parameters (B) and user-defined options/values (C) used for the creation of very-high-resolution (VHR) digital surface models (DSMs) in Agisoft Metashape from collected aerial imagery.

#	Processing Step (A)	Parameter (B)	User-defined Option/Value (C)
1	Selection of aerial images	Image quality (IQ) check	Images with IQ < 0.8 removed
2	Align photos	Accuracy	High
		Pair selection	Reference
		Key point limit	40.000
		Tie point limit	10.000
3	Sparse point cloud filtration (gradual selection)	Reprojection error	< 0.27
		Projection accuracy	< 6
		Reconstruction uncertainty	< 23
4	Point cloud optimization	Optimization parameters	All parameters
5	Build dense cloud	Quality	Low
		Depth filtering	Aggressive
6	Build mesh	Surface type	Arbitrary
		Face count	High
		Interpolation	Enabled
7	Adding GCPs and CPs	7 GCPs and 7 CPs added
8	Point cloud optimization	Optimization parameters	All parameters
9	Sparse point cloud filtration (gradual selection)	Reprojection error	< 0.27
		Projection accuracy	< 6
		Reconstruction uncertainty	< 23
10	Point cloud optimization	Optimization parameters	All parameters
11	Build dense cloud	Quality	High
		Depth filtering	Aggressive
12	Build mesh	Surface type	Arbitrary
		Face count	High
		Interpolation	Enabled
13	Build texture	Mapping mode	Generic
		Blending mode	Mosaic
		Texture size	8096
		Color correction	Enabled
14	Build DEM	Coordinate system	HTRS96
		Source data	Dense cloud
		Interpolation	Enabled
		Point classes	All
15	Build orthomosaic	Surface mode	DEM
		Blending mode	Mosaic

. The image workflow process was performed according to the guidelines of the SfM photogrammetry in geomorphic research [108].

Validation of Model Accuracy and Uncertainty

Detection and interpretation of real spatio-temporal changes are directly affected by the accuracy of created VHR models, as various errors and artifacts can lead to misinterpretations [109,110]. In order to detect and extract real morphological changes, DSM uncertainty was accounted for through thorough accuracy assessment.

Within this research, accuracy assessment was based on seven CPs which were used for the determination of model uncertainty, through the calculation of mean absolute error (MAE) and root-mean-square error (RMSE) [37,44]. Calculated MAE and RMSE values were used for the determination of a minimal level of detection (${\mathrm{LoD}}_{\min }$), which was applied as a threshold for the separation of real morphological spatio-temporal changes from changes induced by systematic and/or nonsystematic errors in created models. ${\mathrm{LoD}}_{\min }$ was determined for both created DSMs, as the root sum square of errors, based on calculated MAE and RMSE values [111]:

$${\mathrm{LoD}}_{\min \left(RMSE\right)}=\sqrt{{\left(\delta xy{z}_{SurveyB}\right)}^2+{\left(\delta xy{z}_{SurveyA}\right)}^2}$$

(1)

$${\mathrm{LoD}}_{\min \left(MAE\right)}=\sqrt{{\left(\delta xy{z}_{SurveyB}\right)}^2+{\left(\delta xy{z}_{SurveyA}\right)}^2}$$

(2)

where:

• ${\mathrm{LoD}}_{\min \left(RMSE\right)}$ = minimal level of detection calculated as root sum square of errors based on RMSE values
• ${\mathrm{LoD}}_{\min \left(MAE\right)}$= minimal level of detection calculated as root sum square of errors based on MAE values
• $\delta xy{z}_{SurveyA}$ = error (RMSE or MAE) for survey A calculated in used CPs
• $\delta xy{z}_{SurveyB}$ = error (RMSE or MAE) for survey B calculated in used CPs

Final absolute detection threshold ${\mathrm{meanLoD}}_{\min }$ which represents the minimal intensity of spatio-temporal changes that VERTICAL can detect and quantify from created interval DSMs, was determined as the mean of two calculated minimal levels of detection ${\mathrm{LoD}}_{\min \left(RMSE\right)}$ and ${\mathrm{LoD}}_{\min \left(MAE\right)}$), and applied to the erosion rates determined by VERTICAL.

As CPs only demonstrate only point-based accuracy of created DSMs accuracy of created models was further evaluated through the validation of point clouds with M2C3 plugin of CloudCompare software [87]. Using the M2C3 plugin, two-interval point clouds were compared within three test sites (2 × 2 m) chosen over limestone bedrock, where no spatio-temporal changes did not occur within the study period. On average around 12,000 points were used for the comparison of point clouds within every test site.

3.4. Development of the VERTICAL Method

VERTICAL was developed in ModelBuilder application within ArcGIS 10.1 software [112], as a toolbox that combines various existing tools from 3D Analyst and Spatial Analyst ArcGIS extensions. It requires the following input parameters for the automated creation of GCs and quantification of spatio-temporal changes: main sampling line (MSL), gully (study) area, and two interval DSMs (DSM 1 and DSM 2). MSL is a very important parameter since it represents a user-defined line, whose vertices are used as starting points for the derivation of perpendicular GCs. It can be either a straight line that divides the whole study area into two identical sections, or it can be a specific, user-defined line, such as gully thalweg (Figure 5). In this research 110 m long thalweg of gully Santiš, generated from the flow accumulation model of the gully area is chosen as an MSL. The user-defined gully area parameter restricts the processing extent for VERTICAL. The DSM 1 and DSM 2 are the interval DSMs used for comparison of heights within created GCs and calculation of vertical spatio-temporal changes.

The entire process of VERTICAL method application can be divided into four main phases (Figure 5). In the first phase, the VERTICAL method determines the mean linear direction (LDM) of the input MSL which is set as default orientation (°) of every vertex within MSL (Figure 5A). Before the actual calculation of LDM input MSL is automatically split at vertices, thus allowing calculation of the linear direction for all line segments, and not only for the whole MSL.

In the second step, MSL vertices are extracted, along with their XY coordinates (Figure 5B). Precise coordinates are needed for the later creation of GCs, perpendicular to the determined LDM value. If gully thalweg is used as MSL (as in this study) then the interval between created GCs is defined by extracted vertices that represent curves defined by the flow direction. In that case interval between created GCs is not manually and subjectively determined by the user, but it is rather determined by the morphological and hydrological characteristics of the studied gully. Otherwise, if a straight line is used as an MSL, the user has to define manually the interval at which GCs will be created.

Within the third step, for every extracted MSL vertex VERTICAL automatically constructs linear GCs perpendicular to the determined LDM (Figure 5C). Cross-section lines are constructed from the MSL vertex as a starting point in two opposite perpendicular directions (left and right), according to the following formulas:

A₁ = LDM (°) + 90°

(3)

A₂ = LDM (°) – 90°

(4)

where:

• A₁ = GCs line constructed left from every MSL vertex
• A₂ = GCs line constructed right from every MSL vertex
• LDM = mean linear direction of MSL

Left and right parts of created GCs lines are then merged into a single GCs line and a unique identification number (ID) is assigned to every created GCs. Unique IDs are necessary for later analysis and extraction of separate statistics for every created GCs. The width of the derived GCs is automatically limited by the input gully area.

Then, within the fourth step, VERTICAL uses the created GCs for a sampling of height data (h) from two-interval DSMs (Figure 5D), where the sampling interval within created GCs has to be defined manually by the user. Overall sampling density directly influences the thoroughness of GCs morphology representation and quantification. We propose that the sampling interval should be defined in regard to the research goals, spatial resolution, and accuracy of used DSMs. In that way, the user can adjust the density of height point samples within created GCs according to the research needs and goals, and the quality of available interval digital surface models.

The height difference between two DSMs is calculated for every sampled point according to the following formula:

$$\Delta h=h_{DS{M}_{\mathrm{A}}}-h_{DS{M}_{\mathrm{B}}}$$

(5)

where:

• $\Delta h$ = height difference at the sampled point
• $h_{DS{M}_{\mathrm{A}}}$ = height of sampled point in initial DSM
• $h_{DS{M}_{\mathrm{B}}}$ = height of sampled point in final DSM

The final step of VERTICAL application is an automated calculation of width/depth (W/D) ratio. It is a dimensionless ratio that gives an indication of the GCs shape [24,72,113]. As one of the most frequently used cross-sectional metrics (Table 1), W/D ratio is used for the evaluation of proportion between lateral erosion and incision rates within the GCs, where higher values indicate the predominance of lateral erosion rate over incision rate [24]. Furthermore, the W/D ratio is commonly used for basic, shape-based GCs classification (e.g., V-shaped cross-sections; U-shaped cross-sections). Values of W/D ratio usually vary between 2.0 and 18.0, where lower ratio values are indicating narrower gully channels (V-shaped) and higher values indicate wider channels (U-shaped) [114](114 in 24). VERTICAL automatically calculates a unique W/D ratio for every sampled GCs. Calculation of the W/D ratio (Figure 6) is performed by the VERTICAL accordingly to the following formula:

$$\mathrm{W}/\mathrm{D} ratio=\frac{W_{GCs}}{D_{gcs}}$$

(6)

where:

• $W_{GCs}$ = width of the sampled GCs
• $D_{gcs}$ =maximal depth of the sampled GCs

4. Results

4.1. Spatial Resolution and Accuracy of Created Interval VHR Models

Aerial imagery collected during the conducted aerial surveys was used for the derivation of two-interval VHR DSMs with 1.9 centimeters spatial resolution and digital orthomosaics with 0.5 centimeters spatial resolution. Point density was 2980 points/m² for Survey A and 3180 points/m² for Survey B. Created VHR models can be seen in Figure 7.

Results of accuracy assessment are demonstrating that model uncertainty is very similar for both created VHR DSMs. Total RMSE values calculated for seven fixed checkpoints (CPs) amount to 3.76 cm for Survey A and 3.51 cm for Survey B, resulting in a mean RMSE of 3.63 cm. The total mean absolute error (MAE) was 3.37 cm for Survey A and 3.25 cm for Survey B (mean MAE = 3.31 cm). Reprojection error (RE) for Survey A is 0.250 pix and 0.263 pix for Survey B. Minimal level of detection calculated based on RMSE and MAE values was 5.143 cm (${\mathrm{LoD}}_{\min \left(RMSE\right)}$) and 4.682 cm (${\mathrm{LoD}}_{\min \left(MAE\right)}$, resulting in an absolute detection threshold (${\mathrm{meanLoD}}_{\min }$ of 4.913 cm. Point cloud accuracy validation performed by the M2C3 algorithm resulted in a 0.59 cm mean calculated distance between two point clouds, with a 2.86 cm standard deviation. If obtained values (MAE; RMSE; ${\mathrm{LoD}}_{\min })$ are compared to the values from similar published geomorphic researches [45,87], it can be concluded that the accuracy of created VHR models is sufficient for conducted spatio-temporal changes analysis.

Based on carried accuracy assessment it can be concluded that uncertainty of created interval VHR DSMs is within ±5 cm. As a result, a minimal level of detection threshold for spatio-temporal changes detected by VERTICAL was set to ±5 cm accordingly. Application of ±5 cm threshold allowed separation of real morphological changes (erosion: values above −5 cm; accumulation: values above 5 cm), from changes induced by errors in created models.

4.2. Interpretation of Determined Spatio-Temporal Change

The VERTICAL method was applied for the comparison of DSMs representing the 2-year period between surveys A and B. In total VERTICAL method created and sampled 2379 GCs, along the 110 m long thalweg of gully Santiš (Figure 8). The total number of height points sampled within all created GCs is 749,052, while on average 314.86 points were sampled per every created GCs. Since the spatial resolution of interval DSMs in our research was 1.9 cm and calculated accuracy metrics (RMSE; MAE; $mean{\mathrm{LoD}}_{min}$) were under 5 cm, we defined a 5 cm sampling interval. As a result, VERTICAL sampled the height data from input DSMs along the created GCs lines at 5 cm intervals.

Furthermore, all sampled GCs under 5 cm were excluded from the determination of spatio-temporal changes, as such these values cannot be considered accurate and reliable enough for the analysis, due to the DSM uncertainty. Therefore, from 2379 sampled GCs, only 922 GCs that had derived mean spatio-temporal changes values above the 5 cm threshold were used for the analysis of spatio-temporal changes intensity.

Calculated mean spatio-temporal changes measured within 922 GCs vary significantly from the gully headcut to the gully terminus (Figure 8). This correlates with the overall heterogeneity of the gully erosion process [92]. Erosion prevails in 356 GCs (38.61%), with the mean erosion rate per sampled GCs between 5.13 ± 5 and 51.65 ± 5 cm in the study period⁻¹. On the other hand, accumulation prevails in 566 GCs (61.38%), with the mean accumulation rate per sampled GCs between 5.04 ± 5 and 13.43 ± 5 cm in the study period⁻¹. Variations in extracted spatio-temporal change values reflect the different processes and trends that have occurred between the two surveys. Although gully Santiš represents a relatively simple and short gully, results derived by VERTICAL are showing that spatial distribution and intensity of the processes that influence its formation are more complex than expected. For example, within the first 51 analyzed GCs (Figure 8A) maximal mean erosion values were recorded, ranging between 9.21 ± 5 and 51.64 ± 5 cm in the study period⁻¹. The maximal erosion value measured within point samples of that 51 GCs is 1.43 ± 5 m in the study period⁻¹. Such pronounced erosion is a result of mass movements related to the collapse and uphill progression of the gully headcut (Figure 8 and Figure 9). Pronounced uphill progression and intensive erosion rates are typical for most gully headcuts [15], and values of measured spatio-temporal changes within the first 51 GCs are strongly indicating that this is true for the headcut of gully Santiš.

Plotted point samples from ${\mathrm{GCs}}_{50}$ confirm that between two surveys, intensive headcut collapse occurred within these GCs (Figure 9A,D,F), as well as accumulation of eroded material at the headcut base (Figure 9B,C,E). Uphill progression of the gully headcut is also evident from the horizontal shift of the line representing headcut position in Survey A and in Survey B (Figure 9). Due to the erosion-induced material collapse mean linear gully headcut retreat (GHR) between Survey A and B was 10.85 ± 5 cm study period⁻¹, while maximum values amounted up to the 49.71 ± 5 cm study period⁻¹. Since these values represent the two-year period, the maximum annual GHR amounts to 24.85 ± 5 cm year⁻¹. Determined annual GHR corresponds to the rates recorded in other researches, with similar conditions [15,115,116].

After the headcut and initial prevalence of erosion, accumulation prevails within the next several hundred GCs, until the ${\mathrm{GCs}}_{886}$ (Figure 8B). In this section, the gully gradually narrows, forming one main channel. Steep sidewalls on both sides of the main channel are being slowly eroded (Figure 10A), similarly as at the gully headcut, while eroded material is accumulated within the channel (Figure 10B–E). The prevalence of accumulation within this area is related to the accumulation of the soil material dispatched from the nearby sidewalls and main headcut. Due to the lack of stronger surface flow eroded material cannot be transported much further away from the base of sidewalls and headcut [21], thus being accumulated within the first 20 meters of the gully. Experimental research conducted by [21] have demonstrated the importance of intense rainfall-induced torrential surface flows on transportation of discharged sediment. They concluded that lower-intensity rainfall can be sufficient for the occurrence of mass movements at gully headcut, but at the same time insufficient for further transportation of material from the headcut. This causes a gradual accumulation of the eroded material. Our onsite precipitation measurements between May and November 2017 registered 626 mm of rainfall with highly heterogenic distribution, which is in very good concordance with the average values and seasonal distribution of the rainfall recorded by [97]. Summer months (June-August) record very low and relatively uniform amounts of precipitation, while autumn months recorded uneven distribution and significantly increased precipitation values. Considering that no torrential rain (I > 60 mm/h) was recorded, it can be presumed that also no torrential surface runoff occurred within the studied period. Therefore, it is possible that surface runoff intensity was insufficient for further transportation of the material. Consequently, such rainfall intensity corroborates data on soil accumulation within the first 20 meters after the headcut (Figure 8B). However, as noted by the field observations, the rain events recorded in September 2017 (96.8 mm and 140 mm) had a significant influence on the changes in gully morphology. Such changes proved the importance of the occurrence of intense rainfall-induced surface runoff for reshaping the gully morphology [21,91].

Sporadic stronger erosion values were recorded within the middle part of the gully (Figure 8C), which is related to the intensive collapse of steep walls and incision of smaller sub-channels, formed by concentrated periodic surface runoff within this part of the main gully channel. Within the lower part of the gully (Figure 8D) low-intensity erosion prevails, with erosion values in GCs up to the 11.61 ± 5 cm study period⁻¹. The mean erosion value for this section (${\mathrm{GCs}}_{1234}-{\mathrm{GCs}}_{2057}$) of gully Santiš was 6.17 ± 5 cm study period⁻¹. The absence of soil accumulation and less homogenous parental material have led to the occurrence of selective erosion within this part of the gully.

The last section of gully Santiš (Figure 8E) is influenced by the fluctuations and dynamics of the Adriatic Sea, primarily by the influence of wave activity. Generally, the eastern Adriatic Sea has microtidal characteristics, with a tidal amplitude between 0.22 and 1.2 m [117] (117 in 119). According to [118], mean annual significant wave heights within the Pag island quadrant are between 0.3 and 0.4 m, while the mean annual wind speed (Bf) ranges from 1.5 m in S, SE, and SW part, up to 2.5 m in the rest of the Island. In spite of the rather low values, the influence of waves on coastal areas is present. During the field surveys accumulation of pebbles and soil sediment was observed at the beach, while these sediments were later disrupted and dislocated by the Sea activity. Our observation corresponds well with the time of percentage with no wind in the Pag island quadrant, which is between 10 and 20% [118]. Moreover, even though the coastline is mostly formed by carbonate rocks, which are more prone to karst processes with negligible effects of coastal erosion [119], the existing pocket beach consisting of gravel deposits (Figure 11) proves coastal dynamics and erosion of soil deposits.

Results also indicate that smaller pebbles accumulated at the beach (Figure 11C,D) were dislocated by Sea activity, while larger, heavier boulders remained at the same location (Figure 11E). Another confirmation of the accuracy of used DSM models is the anthropogenic material (e.g., timber, plastic pipes), accumulated at the beach by Sea activity between two surveys. This anthropogenic material was detected and measured by VERTICAL as the accumulation of material (Figure 11A,B).

We assume that the strongest wave activity and related processes occur under the influence of Bora wind [89] which in this part of the Adriatic raises waves up to 7.2 meters [120]. Namely, Santiš is oriented in the NW direction and the average number of days with wind from the NE and NNE quadrant is 85.4 and 144.6 respectively, with a mean annual wind speed of 5.5 Bf and maximum of 30.4 Bf [97]. Field observations made a few days after the significant rain event in September 2017 (140 mm) proved the connection between rain-induced surface and subsurface runoff formation and reactivation of gully hydrological function. During the field survey carried out a few days after the significant rain event (140 mm), reactivation of spring at the pebble beach occurred, within the final part of the gully, along with the reactivation of nearby submarine springs (Figure 12B).

4.3. Interpretation of Derived W/D Ratio

VERTICAL calculated automatically W/D ratio for all 2379 sampled GCs. Values of calculated W/D ratio range between 5.23 and 7.91 for survey A, and between 5.13 and 8.15 for survey B.

The highest W/D ratio values for both survey A and B were present at the initial part of gully Santiš, from gully headcut towards the middle part. From the middle part, the gully narrows, and W/D ratio values gradually decline towards the pebble beach at the gully terminus. Such transition of W/D ratio values (from U-shaped GCs towards the V-shaped GCs) is opposite to typical general gradual increase in W/D ratio values from headcut towards the gully terminus, which was reported in some other similar case studies [e.g., 77] (77 in 24). This can be further explained by the overall heterogeneity of soil sediments within the gully Santiš. While the initial part of the gully is formed in homogenous, erosion-prone soil sediments, the middle and final parts of the gully are formed in sediments with pronounced erosion resistance. Thus, as calculated W/D ratio values indicated, strong lateral erosion that is present in the initial part of the gully (Figure 12A) is gradually subsided by more pronounced channel incision and occurrence of selective erosion towards the middle and final parts of the gully (Figure 12B).

5. Discussion

5.1. Advantages of VERTICAL over the Profilometer

Regardless of its design, the profilometer has few significant shortcomings which affect the practicality of its application and accuracy of its measurement. First, an important shortcoming is a length constraint (1) of the GCs that the profilometer can measure. Due to its design profilometer is suitable only for the measurement of smaller erosion features, e.g., rill and ephemeral gullies. For gullies with wider cross-sections, such as gully Santiš, a profilometer should be more than 10 m long, which is very impractical. A mechanical profilometer is affected by depth constraint (2), as the maximal depth that the profilometer’s pins can measure is restricted by the length (scale) of the pins. If the intensity of spatio-temporal changes overcomes the maximal scale of measurement pins, a profilometer cannot be used for the measurement of these changes. Furthermore, the whole measurement process with a profilometer is time-consuming and labor-intensive (3), since measurements are restricted to one GCs at a time. This significantly increases the duration of field measurements, especially for longer gullies and researches that require sampling of a larger number of GCs. The next shortcoming of a profilometer regards consistency, accuracy, and precision (4) of its measurements. While optical profilometers tend to underestimate depth values [27], mechanical profilometers overestimate them [79]. Errors caused by optical profilometer are sensor related [27], while errors caused by mechanical profilometer can be related to the direct measurement of depth values (e.g., compaction problem [29,121,122]. Furthermore, in order for interval measurements to be accurate and precise, the profilometer has to be placed at the exact XYZ location, which is often difficult to accomplish in the field. Even the slightest deviation along any axis (X, Y, or Z) in device placement can lead to deviation of measurement values and occurrence of significant measurement errors.

The developed VERTICAL method represents a fast and accurate measurement method that allows automated measurement of vertical spatio-temporal changes within a very large number of GCs. The main novelty and advantage of the VERTICAL method is the automation of the whole measurement process (in comparison to both traditional field techniques and most of existing raster-based methods), which significantly shortens and simplifies the detection of spatio-temporal changes, while at the same time improving the accuracy and repeatability of GCs measurement.

Furthermore, the application of VERTICAL demonstrated that it successfully resolves most of the stated limitations of profilometer and similar traditional field sampling methods. The main advantages of VERTICAL are:

Indirect measurement (1)—VERTICAL allows measurement of spatio-temporal changes in GCs without direct physical contact with the measured surface. Measurements are non-destructive, meaning that the measurement process will not alter the state of a measured surface (e.g., compaction occurrence), as some other traditional methods would.

High sampling density (2)—this method measures far more GCs and height samples that would be possible to achieve with traditional field methods. A total number of sampled GCs and height points are defined by the user-defined MSL and sampling interval. In theory, VERTICAL can measure an unlimited number of GCs and height points. On the other hand, the measurement of such a large number of GCs and height samples would be hardly achievable with any traditional field measurement technique. Studies that have used profilometer measured significantly less GCs per one gully (e.g., 1 [78], 3 [24,25], 5 [27], 28 [30]) and less height samples within particular GCs (e.g., 46 [78], 50 [26], 100 [27,34]) (Table 1.).

Absence of length and depth constraint (3)—minimal and maximal lengths of GCs that can be sampled by VERTICAL are restricted only by the user-defined MSL and gully area. The profilometer’s depth constraint is also resolved, as VERTICAL can measure any height difference from two defined DSMs. This means that no matter how wide GCs are, or how intensive certain spatio-temporal change at a given location is, VERTICAL will be able to measure it. The length of 2379 GCs sampled within gully Santiš varies from just a few meters (minimum GCs length = 1.91 m), up to tens of meters (maximal GCs length = 30.34 m). Measured height difference at some points sampled within created GCs amounted up to over a meter (maximal height difference = 1.43 m). Measurement of GCs with that variability in length and depth would be hardly achievable with a profilometer (e.g., max GCs length = 1 m [26]; 2 m [27]; 5 m [34]).

Measurement speed and efficiency (4)—in comparison to traditional field techniques, VERTICAL significantly simplifies and shortens the overall time-consuming and labor-intensive process of GCs measurement, where the whole measurement process is fully automated and the required processing time is restricted solely by available processing-power.

Measurement consistency (5)—as VERTICAL uses identical GCs and points with precise XYZ coordinates for measurement of height samples from both interval DSMs, all user-influenced deviations in measurement accuracy and precision are eliminated. This was hard to achieve in the field, as the profilometer had to be positioned et exactly the same location for every measurement.

Adaptability and flexibility (6)—VERTICAL can be applied for different research purposes and areas, not exclusively for measurement of gully erosion induced spatio-temporal changes. This method can be applied on DSMs with different spatial resolution and extend, from sub-millimeter resolution models (e.g., tufa dynamics monitoring, freeze cracking) covering a few cm², sub-meter resolution models (e.g., gully erosion, landslides, river dynamics) covering hundreds of m², up to the medium, or even low-resolution models (e.g., ice mass dynamics, tectonics, bathymetry) covering tens or hundreds of km². Furthermore, VERTICAL can be applied on models created by various different geospatial technologies (e.g., LiDAR, close-range photogrammetry, optical satellite stereo-imagery). An example of successful VERTICAL application for monitoring of tufa formation dynamics is given in Section 5.2.

Freely available (7)—this method is compatible with ArcToolbox and available as an open-source toolset from the official web page of Geospatial Analysis Laboratory (gal.unizd.hr).

5.2. Applicability of VERTICAL in the Study of Tufa Formation Dynamics (TFD)

The VERTICAL method can have prominent applicability in the monitoring of recent tufa and travertine growth and erosion dynamics. This analysis can be carried out at sub micro-level of research (<0.1 mm). Until recently, the tufa formation dynamics (TFD) was predominantly analyzed using a direct measurement method of modified micro-erosion meter (MEM) [123,124,125]. It gathers tufa elevation data in the form of individual points (n = 50) and has various drawbacks [37]. In these researches, seasonal growth dynamics from laminated tufa layers is often analyzed from specific cross-sections interpolating the acquired MEM data. These “visible” sections are created by cutting and removing a part of the formed tufa or travertine from a specific test plate (mostly limestone) which was mounted in tufa forming watercourses. This procedure interrupts the tufa forming process on the studied test surface, which in some cases lasted several years (Figure 13).

From the derived cross-section seasonal differences [125] in tufa formation dynamics (changes in cross-section thickness), differences in structural and textural characteristics of formed tufa, and specific events (macroinvertebrates accumulation and incrustation, erosion due to extreme flood events, dry periods) in the tufa formation process can be detected. However, [37] were the first to use SfM and interval high-resolution digital tufa models (DTHRM) in TFD-a analysis. Given that this new methodological approach studies TFD in the fixed local coordinate system, it is now possible to continuously monitor the dynamics of tufa formation from quantified cross-sections derived from DTHRM. Therefore, VERTICAL, in this case, enables the creation and quantification of cross-sections from generated 2.5D and 3D tufa models thus eliminating the need to destroy the recent tufa and to interrupts the tufa forming process.

Figure 14 and Figure 15 show the elevation values within specific cross-section which were derived from DTHRM using a VERTICAL tool. Two test plates were placed in different fluvial environments and analyzed. The vertical spatio-temporal changes (increment +, erosion −, no change) were tracked in each cross-section with 82 points. On the first test plate, large variability in tufa growth between the initial and sequential DTHRMs can be observed due to a wide range of factors. In the first period, the cross-sections have the highest recorded mean spatio-temporal change (MSTC) because a large number of aquatic insect larvae from the family Chironomidae (order: Diptera) accumulated on the surface. However, since the accumulated organisms and plant material were not well-bounded to the substrate, significant erosion occurred which lowered down the MSTC values measured in the following intervals.

On the second test plate, the formed tufa was cut and partially removed following the classical examples of studying the cross-sections. However, before this, the interval SfM measurement of the formed tufa was made. Then, after the cutting, the remaining tufa was measured and derived cross-sections were added to its 3D model (Figure 15). This example shows how the identification of specific laminated layers (thickness), that have pronounced seasonal characteristics (warm–cold period), can be improved. For example, Figure 13 shows that it is not possible to accurately detect the thickness of the tufa formed in a specific period (e.g., April–September). If the test plate of formed tufa is measured regularly with the SfM approach this drawback is easily addressed.

5.3. Limitations of VERTICAL Method

Despite the stated advantages, the developed VERTICAL method has certain limitations, which should be solved through future upgrades.

The first limitation is its dependence on DSM quality (1), as the accuracy and precision of used VHR DSMs directly affect the measurements performed by VERTICAL. Thus, various systematic and non-systematic DSM errors can cause significant deviation in GCs measurements. As VHR DSMs produced from SfM photogrammetry are influenced by various different factors (e.g., camera calibration, surface texture, lighting conditions, GCP characteristics, vegetation, complex terrain morphology) that can cause the occurrence of errors [44], users should take care of possible influence of such errors on GCs measurement. For example, if vegetation is not excluded from created interval DSMs, VERTICAL can measure vegetation growth, which could then be misinterpreted as erosion-induced STC. However, the application of VERTICAL for monitoring TFD has demonstrated that highly accurate VHR DSMs, which are based on local coordinate systems, significantly reduce the uncertainty of measurements.

The second important limitation of the developed VERTICAL method is its restricted applicability to 2.5D models (2), and this segment will be addressed within future upgrades. Currently VERTICAL is not suitable for the measurement of complex 3D gully morphology, like undercuttings and overhangs, as these are not represented in 2D or 2.5D models.

The third limitation is the required user expertise (3) needed for the derivation of high-quality VHR DSMs that VERTICAL uses for measurement of GCs. Although the application of VERTICAL is fast and user-friendly, the process of creation of VHR DSMs with required accuracy is not simple and straightforward. However, if VERTICAL is applied to already created, available open-source DSMs, the amount of required user expertise is significantly lower.

6. Conclusions

Several important methodological conclusions can be highlighted on the basis of the developed VERTICAL method:

(1)

VERTICAL method successfully overcomes the stated limitations of the profilometer and similar traditional field measurement techniques, as its application facilitates and simplifies the overall process of detection and measurement of vertical spatio-temporal changes within GCs. The extent of the user’s required expertise and its influence on error generation are minimalized with VERTICAL.

(2)

A very large number of GCs sampled by VERTICAL allows a thorough representation of overall spatio-temporal changes in gully geometry. Due to the high sampling density detailed distinction of different complex erosion and accumulation induced processes is possible. Interpretation of measured spatio-temporal changes is possible within the whole gully (Figure 8), or within separate chosen GCs (Figure 9, Figure 10 and Figure 11).

(3)

As demonstrated in Section 5.2, the VERTICAL method is potentially applicable for other, similar scientific purposes, where multi-temporal accurate measurement of spatio-temporal changes in cross-sectional geometry is required (e.g., river material dynamics, ice mass dynamics, tufa sedimentation, and erosion).

Application of the developed VERTICAL method for measurement of spatio-temporal changes allowed detailed reconstruction and quantification of gradual gully development. Following conclusions about detected gully evolution can be highlighted:

(1)

Mean STC values vary significantly within 2379 GCs of gully Santiš from the gully headcut until the gully terminus.

(2)

Highest erosion rates were recorded at the initial part of the gully, where intensive collapse and uphill progression of gully headcut were observed.

(3)

Most of the material eroded from gully headcut and sidewalls is being accumulated within the first 20 meters of the gully. Such accumulation could be related to the lack of stronger surface flow, capable of further transportation of eroded material.

(4)

Less homogenous middle part of the gully is influenced by the occasional occurrence of stronger selective erosion, manifested mainly through the channel incision and sidewall collapse.

(5)

The final part of gully Santiš is influenced by the dynamics of the Adriatic Sea, which are disrupting and dislocating accumulated sediments.

Within future research, the developed VERTICAL method will be further improved, especially in regard to its current limitations (e.g., simplification of the method, less dependency of results to errors in models). Special attention will be dedicated to the data collection of ground control and checkpoints with a total station in order to achieve even better measurement accuracy. Furthermore, VERTICAL will be tested for the detection of spatio-temporal changes within other similar geomorphic purposes (e.g., detection of tufa sedimentation and erosion).