Comparative analysis of the Copernicus, TanDEM-X, and UAV-SfM digital elevation models to estimate lavaka (gully) volumes and mobilization rates in the Lake Alaotra region (Madagascar)

Over the past few decades, developments in remote sensing have resulted in an ever-growing availability of topographic information on a global scale. A recent development is TanDEM-X (TerraSAR-X add-on for digital elevation measurements), an interferometric synthetic aperture radar (SAR) mission of the Deutsche Zentrum für Luft- und Raumfahrt, providing near-global coverage and 12 m resolution digital elevation models (DEMs). Moreover, ongoing developments in uncrewed aerial vehicle (UAV) technology have enabled acquisitions of topographic information at a sub-meter resolution. Although UAV products are generally preferred for volume assessments of geomorphic features, their acquisition remains a time-consuming task and is spatially constrained. However, some applications in geomorphology, such as the estimation of regional or national erosion quantities of specific landforms, require data over large areas. TanDEM-X data can be applied at such scales, but this raises the question of how much accuracy is lost because of the lower spatial resolution. Here, we evaluated the performance of the 12 m TanDEM-X DEM to (i) estimate gully volume, (ii) establish an area–volume relationship, and (iii) determine mobilization rates through comparison with a higher-resolution (0.2 m) UAV structure-from-motion (SfM) DEM and a lower resolution (30 m) Copernicus DEM. We did this for six study areas in the Lake Alaotra region (central Madagascar), where lavaka (gullies) are omnipresent and surface area changes over the period 1949–2010s are available for 699 lavaka. Copernicus-derived lavaka volume estimates were systematically too low, indicating that the Copernicus DEM is too coarse to accurately estimate volumes of geomorphic features at the lavaka scale (100–10⁵ m²). Lavaka volumes obtained from TanDEM-X were similar to UAV-SfM volumes for the largest features, whereas the volumes of smaller features were generally underestimated. To deal with this bias we introduce a breakpoint analysis to eliminate volume reconstructions that suffer from processing errors as evidenced by significant fractions of negative volumes. This elimination allowed the establishment of an area–volume relationship for the TanDEM-X data with fitted coefficients within the 95 % confidence interval of the UAV-SfM relationship. Our calibrated area–volume relationship enabled us to obtain large-scale lavaka mobilization rates ranging between 18 ± 3 and 311 ± 82 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for the six different study areas, with an average of 108 ± 26 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for the full dataset. These results indicate that current lavaka mobilization rates are 2 orders of magnitude higher than long-term erosion rates. With this study we demonstrate that the global TanDEM-X 12 m DEM can be used to accurately estimate volumes of gully-shaped features at the lavaka scale (100–10⁵ m²), where the proposed breakpoint method can be applied without requiring the availability of a higher-resolution DEM. Furthermore, we use this information to make a first assessment of regional lavaka erosion rates in the central highlands of Madagascar.

1 Introduction

Over the past few decades advanced technology has become increasingly available for the assessment of surface topography: structure-from-motion (SfM) algorithms applied to uncrewed aerial vehicle (UAV) imagery now allow centimeter-scale resolution, thereby revolutionizing the way we study Earth surface processes (Passalacqua et al., 2015; Tarolli, 2014; Clapuyt et al., 2016). Obtaining sub-meter-resolution digital elevation models (DEMs) from UAV-SfM still requires substantial fieldwork and is spatially limited due to the nature of the technology (Bangen et al., 2014). On the other hand, TanDEM-X (TerraSAR-X add-on for digital elevation measurements) is a spaceborne product with global coverage at 12 m resolution and, while being less detailed and accurate than these sub-meter-resolution DEMs, is a major step forward in comparison to 30 m resolution DEMs with a global coverage (Mudd, 2020).

This raises the question to which extent TanDEM-X imagery can be used to map three-dimensional morphological features requiring a higher degree of topographical detail over relatively large areas (> 10 km²). One example is the use of remotely sensed data to map the process of gully formation, which is known to significantly contribute to surface erosion (e.g., Poesen et al., 2003; Vanmaercke et al., 2021; Frankl et al., 2021). The mapping and monitoring of gully erosion was conventionally based on time-consuming and spatially limited field surveys (Castillo et al., 2012; Evans and Lindsay, 2010; Guzzetti et al., 2012). More recently, however, (sub-)meter-resolution DEMs have enabled the development of (semi-)automated gully delineation and volume determination methods (Niculiţă et al., 2020; Evans and Lindsay, 2010; Perroy et al., 2010; Eustace et al., 2009; Liu et al., 2016). TanDEM-X has, for example, already been successfully used for automatic gully detection (Vallejo Orti et al., 2019).

It is important to evaluate not only the extent to which TanDEM-X data can be used to estimate gully volumes but also its capacity to establish accurate area–volume relationships. This latter question is relevant since sub-meter-resolution surface imagery from a multitude of sources and moments in time is now globally and freely available through, for example, Google Earth (Fisher et al., 2012). This imagery can be used to identify geomorphic features and estimate their surface area with great detail. Area–volume or length–volume relationships then enable us to obtain estimates of volume changes over time when historical imagery from which areas or lengths can be derived is available. Work on gully and landslide erosion has shown that the establishment of these relationships enables us to estimate sediment mobilization rates (i.e., the average annual volume of hillslope material displaced per unit area) over large spatial and temporal scales (e.g., Malamud et al., 2004; Larsen et al., 2010; Hovius et al., 1997; Guzzetti et al., 2012; Broeckx et al., 2019; Frankl et al., 2013; Kompani-Zare et al., 2011). Furthermore, work on gully headcut retreat rates and associated sediment mobilization has indicated the importance of long measurement periods: large year-to-year variations result in very large (> 100 %) uncertainties over short (< 5 years) measuring periods (Vanmaercke et al., 2016).

Here, we evaluate the performance of TanDEM-X to estimate gully volumes and to establish area–volume relationships through comparison with a higher resolution (0.2 m) UAV-SfM DEM and a lower resolution (30 m) Copernicus DEM. The resolution of a DEM should be viewed in relation to the size of the landform, where sampling theory states that landforms should have dimensions of at least twice the DEM resolution (Theobald, 1989; Frankl et al., 2013). This would mean that a gully should have a theoretical minimum size of 0.16, 576, and 3600 m² for the UAV-SfM, TanDEM-X, and Copernicus DEM, respectively, to be represented in the DEM. We used the lavaka of the central highlands of Madagascar as a case study. Lavaka are amphitheater-shaped gullies that are on average 60 m long, 30 m wide, and 15 m deep, with small outlets (Cox et al., 2010; Wells and Andriamihaja, 1993). They are omnipresent in the central highlands, leaving the landscape filled with “holes” (“lavaka” in Malagasy). Unlike other gullies, they typically lack surface feeder channels and tend to form on mid-slopes, broadening uphill through headward erosion (Wells et al., 1991; Wells and Andriamihaja, 1993). While Madagascar is often claimed to experience amongst the highest global erosion rates due to the presence of lavaka (Milliman and Farnsworth, 2011; Randrianarijaona, 1983), the amount of sediment that is directly produced by lavaka is currently unknown.

Figure 1 Study areas, examples of each digital elevation model (DEM), and lavaka examples. (a) Six study areas of ca. 10 km² in the Lake Alaotra catchment shown on the TanDEM-X DEM with hillshade. The UAV-SfM DEM is available for study area 5 and 6 (red) (data collected June 2018). (b–d) Examples of the Copernicus (AIRBUS, 2020a), TanDEM-X (Krieger et al., 2007), and UAV-SfM DEM (data collected June 2018) in SA6 located at 17^∘58^′51.7^′′ S, 48^∘15^′18.6^′′ E with hillshade. Grey outlines indicate the digitized lavaka. (e) UAV fisheye picture from 200 m height of the eastern ridge shown in (b–d). (f) UAV fisheye picture (200 m height) from two typical amphitheater-shaped lavaka (pictures taken in June 2018 by Liesa  Brosens).

The objectives of our study are therefore to evaluate the performance of TanDEM-X to (i) estimate lavaka volumes, (ii) establish accurate area–volume relationships, and (iii) obtain a first estimate of lavaka mobilization rates. We derived lavaka volumes and mobilization rates for an existing dataset containing 699 digitized lavaka in six study areas in the Lake Alaotra region at three moments in time: 1949, 1969, and the 2010s. In a first step, lavaka volumes were calculated for the 2010s lavaka polygons from the DEM as the difference between a reconstructed pre-erosion surface and the current topography. Following this, a lavaka area–volume relationship was established between the current lavaka areas (2010s) and calculated volumes. Finally, this relationship was applied to the historical dataset with lavaka areas in 1949, 1969, and the 2010s. This enabled us to calculate lavaka volumes at each of these time steps and the consequent derivation of volumetric growth rates and lavaka mobilization rates for each of the six study areas. This procedure was followed for a 0.2 m resolution UAV-SfM DEM, the 12 m resolution TanDEM-X DEM, and the 30 m resolution Copernicus DEM.

2 Material and methods

2.1 Study area and lavaka dataset

Six study areas (SAs) of ca. 10 km² were selected in the northeastern part of the central highlands of Madagascar in the area surrounding Lake Alaotra (Fig. 1a). The lake is located in the seismically active NE-SW oriented Alaotra–Ankay graben structure and is surrounded by convex-shaped, deeply weathered (> 10 m), and regolith-covered hillslopes that developed on the Precambrian crystalline basement (Kusky et al., 2010; Riquier and Segalen, 1949; Bourgeat, 1972; Mietton et al., 2018). The climate is characterized by a distinct dry and wet season with the regular occurrence of tropical cyclones. The mature rounded hillslopes are covered with open grasslands and contain one of the highest lavaka densities of the country, with up to 14 lavaka per square kilometer (Cox et al., 2010; Voarintsoa et al., 2012). The lowland areas bordering the lake consist of swamps in the SW and vast rice fields elsewhere, producing the majority of rice for the country (Lammers et al., 2015). For the six selected study areas, digitized lavaka polygons were generated from orthorectified and georeferenced historical aerial images from 1949 and 1969 (2.4 m resolution) and from recent (2011–2018, referred to as the 2010s) satellite imagery (WorldView-2, 0.5 m resolution) (Fig. 1a, Table S1 in the Supplement). All shapefiles have WGS84-UTM 39S (EPSG: 32739) coordinates. The dataset (available at https://doi.org/10.6084/m9.figshare.c.5236322.v1) contains the changes in surface area of 699 lavaka over the period 1949–2010s for SA1–5 and 1969–2010s for SA6. Each study area contains 50 to 173 lavaka, resulting in lavaka densities between 4 and 17 lavaka per square kilometer (Table S1).

2.2 Digital elevation models (DEMs)

Lavaka volumes were determined from three digital elevation models with a range of horizontal resolutions. For two study areas a 0.2 m resolution UAV-SfM DEM was obtained from a field campaign in 2018. For all study areas the 12 and 30 m resolution TanDEM-X and Copernicus DEMs are available. All DEMs were transformed to WGS84-UTM39S (EPSG: 32739) coordinates using a nearest-neighbor resampling method.

2.2.1 UAV-SfM DEM (0.2 m)

For study area 5 and 6 (Fig. 1a in red) a UAV field survey was carried out in June 2018 to obtain a 0.2 m resolution DEM. In order to cover a large area during a limited amount of time with a high spatial resolution, the post-processing kinematic (PPK) georeferencing approach as developed by Zhang et al. (2019) was used. The UAV images were directly georeferenced by using a RTK (real-time kinematic) receiver on the UAV that was connected to a RTK base station.

We used a custom-made quadcopter UAV with DJI N3 flight controller and fisheye action camera (Go Pro Hero 3, 12 megapixels, 4000 pixel × 3000 pixel, with a 2.92 mm F/2.8 123^∘ HFOV lens). A compact Tallysman TW2710 multi-Global Navigation Satellite System (GNSS) RTK receiver antenna (Reach RTK kit, Emlid Ltd, 23 cm height) was mounted on an aluminum plate centered above the camera. The RTK receiver was connected to a single-board computer in order to synchronize the GPS time with the geotagging of the images. The RTK base station (Emlid Ltd) provided the positioning correction input. It was mounted on a tripod and located at a fixed position at the center of each study area during the flights. Flights were carried out at 200 m height at a speed of 8 m s⁻¹. Pictures were taken every 3.8 s, resulting in an average ground sampling distance of 0.17 m (Fig. 1e and f).

The raw RTK-GPS data from the receiver were corrected with the data from the base station and post-processed using the RTKLib package (Takasu and Yasuda, 2009). A fix-and-hold method with 20^∘ satellite elevation mask was used to correct the positions and geotag the images. The geotagged images were processed with Pix4D software using the default settings, with the vertical and horizontal accuracy set at 0.5 m. For the generation of the DEM, no surface smoothing or filtering was applied. The resulting DEM for both study areas has a resolution of ca. 0.2 m (Fig. 1b).

This method was reported to result in a robust and accurate alternative for georeferencing based on ground control points (GCPs) with a mean absolute error (MAE) of 0.02 m and root-mean-square error (RMSE) of 0.03 m for the vertical accuracy and a precision of 0.04 m (Zhang et al., 2019). Comparable studies over relatively flat areas with an UAV-RTK setup report similar vertical accuracies with RMSE values between 0.03 and 0.07 m (Taddia et al., 2020; Stott et al., 2020). UAV-SfM surveys with GCPs over more complex terrain report higher RMSE values of between 0.10 and 0.45 m (Clapuyt et al., 2016; Cook, 2017). Given the reported high accuracies of optical acquisitions that are georeferenced with RTK-GPS data, this DEM surface can be considered as the reference of the “true” elevation (Grohmann, 2018). In this study we therefore consider the UAV-SfM DEM as the ground-truth reference.

2.2.2 TanDEM-X DEM (12 m)

The TanDEM-X (TerraSAR-X add-on for digital elevation measurements) mission was launched in 2010 by the public–private partnership between the German Aerospace Center (DLR) and EADS Astrium GmbH (Krieger et al., 2007). Its configuration consists of two synthetic aperture radar (SAR) satellites flying in close formation and thereby forming a large X-band single-pass interferometer. The resulting global DEM has a horizontal resolution of 0.4 arcsec (ca. 12 m) and aims at < 2 m relative and < 10 m absolute vertical accuracy (Krieger et al., 2007, 2013; Wessel, 2016). The final TanDEM-X DEM was published in 2016 and consists of data collected between December 2010 and early 2015, where all land surfaces were imaged at least twice and up to 7 or 9 times in difficult terrain (Rizzoli et al., 2017, Fig. 1c). Our study areas were imaged 5 to 9 times, with an average of 7 ± 1, indicating a good coverage. The good performance of the TanDEM-X DEM has been reported, with a final global absolute vertical accuracy of 3.49 m and relative vertical accuracy of 0.99 and 1.37 m on flat (< 20^∘) and steep (> 20^∘) terrain, respectively (Rizzoli et al., 2017). These results are in line with Wessel et al. (2018) and Purinton and Bookhagen (2017), who reported absolute vertical accuracies of 0.20 ± 1.5 and −1.41 ± 1.97 m, respectively.

2.2.3 Copernicus DEM (30 m)

The 1 arcsec (ca. 30 m) resolution global Copernicus DEM (GLO-30) was released in 2021 by the European Space Agency (ESA) and AIRBUS. The DEM is based on the WorldDEM™, which is in turn based on edited and smoothed radar satellite data acquired during the TanDEM-X mission (AIRBUS, 2020a). The reported global absolute vertical accuracy is 2.17 m with an RMSE of 1.68 m. The relative vertical accuracy is smaller than 2 m for < 20^∘ slopes and less than 4 m on > 20^∘ slopes (AIRBUS, 2020b). Given its recent release, only limited additional validation has been carried out (Guth and Geoffroy, 2021). A lower absolute vertical accuracy of GLO-30 has been reported for mountainous areas in Europe with RMSE values between 7 and 14 m (Marešová et al., 2021). These estimates should, however, be viewed as maximum estimates as these high-relief terrains are one of the most challenging settings for DEM acquisitions. Upon comparison of different global 1 arcsec DEMs, Purinton and Bookhagen (2021) concluded that the Copernicus DEM provides the highest-quality landscape representation and should be the preferred DEM for topographic analysis in areas that lack higher-resolution DEMs. They furthermore report a high inter-pixel consistency for both the TanDEM-X and Copernicus DEM, indicating low relative vertical errors for these DEMs.

Figure 2 Lavaka volume determination workflow. Lavaka volumes were calculated for each individual lavaka following an automated workflow for the three studied DEMs: UAV-SfM (0.20 m, top row, data collected June 2018), TanDEM-X (12 m, middle row, Krieger et al., 2007) and Copernicus (30 m, bottom row, AIRBUS, 2020a). (a) The digitized lavaka outline (grey), manually determined horseshoe-shaped polygon on the unaffected hillslope surrounding the lavaka, and current DEM are the three required inputs for the automated volume procedure. The DEM pixels that are not affected by erosion are clipped from the DEM with the horseshoe-shaped polygon and a single point per pixel is generated. The pre-erosion surface is then reconstructed by interpolating between these points (step 1). Two interpolation methods are shown as an example here: triangulated irregular network (TIN) (b) and regularized spline (c) interpolation. The grey polygon indicates the outer edge of the interpolated area. The elevation difference between the interpolated pre-erosion surface and current DEM surface is then calculated, which is clipped to the lavaka extent (Step 2–3). (d) Cross sections of transect A for the DEM, TIN, and regularized spline interpolation.

2.3 Lavaka volume quantification

Individual lavaka volumes were determined as the difference between the current surface and a reconstructed pre-erosion surface. This was done by developing an automated workflow in PyQGIS written in QGIS version 3.16.10 (code and example dataset available at https://doi.org/10.5281/zenodo.5768418). The automated PyQGIS workflow consists of four steps, which are explained in detail below.

Input data

Three input files are required to run the automated volume procedure: (i) a shapefile containing the digitized lavaka (gully) outlines, (ii) a DEM raster file, and (iii) a shapefile containing the digitized outlines of the region surrounding the lavaka that is not affected by erosion. A manual delineation procedure was followed to obtain this surface unaffected by erosion, where a horseshoe-shaped polygon was drawn around each individual lavaka on the hillslope parts that were unaffected by erosion (Fig. 2a).

Step 1: interpolate the pre-erosion surface

First, the DEM raster layer is clipped with the horseshoe-shaped polygon in order to extract the pixels not affected by gully erosion. All pixels that fall within this polygon are extracted in order to have a minimum width of one pixel. Next, one point per clipped DEM pixel is generated and used as input for the interpolation. Finally, these points are used to interpolate the pre-erosion surface. Five interpolation methods were tested, of which the method with the lowest error was applied to the lavaka dataset (see Sects. 2.4.1 and 3.1). Examples of the interpolated pre-erosion surface are shown in Fig. 2 for triangulated irregular network (TIN; see Fig. 2b) and regularized spline (Fig. 2c) interpolation.

Step 2: calculate elevation difference

The current DEM is subtracted from the interpolated pre-erosion surface. The result is a difference raster, with positive values indicating a current surface that is lower than the reconstructed pre-erosion surface. Negative values indicate that the current topography is higher than the reconstructed topography.

Step 3: elevation difference clipped to lavaka extent

The lavaka extent, which is given by the digitized lavaka polygons from the 0.5 m resolution WorldView-2 imagery from 2011–2018, is clipped from the elevation difference raster. In this way a raster with the elevation difference over the lavaka area is obtained (Fig. 2b and c). If the lavaka is smaller than a single pixel (0.04, 144, and 900 m² for the UAV-SfM, TanDEM-X, and Copernicus DEM, respectively), the resulting raster is empty and no volume can be calculated.

Step 4: export results

The unique value report of the lavaka elevation difference raster is exported. It contains the unique elevation values, their count, and the dimensions of the raster pixels. These results are used to calculate the volumes of each lavaka.

A manual delineation of the region surrounding the lavaka that is unaffected by erosion was preferred over an automated pre-erosion interpolation (e.g., Evans and Lindsay, 2010) or interpolation based on the lavaka outlines for two reasons. First, digitized lavaka outlines from aerial imagery are often located on DEM pixels that already have lower values. This is especially the case for the coarser-resolution DEMs due to surface smoothing (e.g., Fig. 1c and d). Second, the very dense presence of lavaka often results in a highly dissected topography and a near absence of topography not affected by lavaka erosion, requiring a precise identification of the areas not affected by erosion. For some lavaka, no horseshoe-shaped polygons could be delineated. In other cases, lavaka were grouped in one enveloping polygon when they were located next to each other (Fig. S1 in the Supplement).

The lavaka volume was calculated from the exported values report as the product of the elevation difference with the pixel area. Both positive and negative elevation differences occur, where a positive difference indicates that the current surface is lower than the reconstructed pre-erosion surface. Negative values indicate that the current topography is higher than the reconstructed topography. Given the presence of both positive and negative elevation differences, three types of volumes (V, m³) could be calculated for each lavaka: (i) the positive volume (V_pos) calculated by summing only the positive elevation differences, (ii) the negative volume (V_neg) calculated from the negative elevation differences, and (iii) the total volume (V_tot) as the difference between the positive and negative volume. If not further specified, the term “volume” refers to the positive volume.

The lavaka volumes as obtained from the different DEMs were compared in two ways: (i) pairwise comparison using only the lavaka for which the volume was determined for each DEM enabling a one-to-one comparison, and (ii) using all the data in order to establish the most robust relationships calibrated on a higher number of observations. The number of lavaka for which the volume was determined depends both on the availability of the UAV-SfM DEM (only in SA 5 and 6) and the resolution of the DEM. No volume could be calculated for the smallest features using the 12 and 30 m resolution DEMs.

2.4 Volume uncertainty assessment

Estimated lavaka volumes, determined as the difference between an interpolated pre-erosion surface and the current DEM surface, will entail a number of uncertainties and errors. Given our application, we address two types of uncertainty or error in our volume estimates: (i) the interpolation error and (ii) the relative height error of the DEM.

2.4.1 Interpolation error

The pre-erosion surface of the lavaka is reconstructed by interpolating between the DEM pixels (one point per pixel) that are not affected by gully erosion. While it is impossible to assess the real interpolation error for a lavaka, as the pre-erosion topography is simply unknown, this error can be estimated by interpolating the surface at locations where no lavaka are present and the pre-erosion surface is thus known. This method, which is similar to Bergonse and Reis (2015), not only allows us to estimate the uncertainties in derived lavaka volumes but also allows for the objective selection of the best interpolation method for a given topographic setting.

Five different lavaka polygons with sizes that span the range of our lavaka dataset were selected: 100, 1000, 5000, 10 000, and 20 000 m². Each polygon was duplicated 10 times, resulting in 50 lavaka polygons. These polygons were placed on unaffected convex-shaped hillslopes on which lavaka typically occur, together with the corresponding horseshoe-shaped polygons. Our automated lavaka volume quantification method (Sect. 2.3) was then applied to this dataset using five different interpolation methods, where the difference between the interpolated surface and the DEM was calculated. This elevation difference gives the interpolation error, as a perfect interpolation would result in a surface identical to the original DEM.

We tested five commonly used interpolation methods for continuous data with their default parameter settings. The first two algorithms are based on a linear method where Delaunay triangles are constructed from the nearest neighbor points, resulting in non-smooth surfaces: (i) linear interpolation (GDAL, 2021) and (ii) TIN interpolation (QGIS, 2020). We also tested three spline interpolation methods that support the creation of curved surfaces in areas without data: (iii) bilinear spline (GRASS, 2021), (iv) bicubic spline (GRASS, 2021), and (v) regularized spline with tension (GRASS, 2003). The bilinear and bicubic spline interpolations are 2D piece-wise non-zero polynomial functions calculated within a limited 2D area, where the Tykhonov regularization parameter affects the smoothing of the surface (smoothing = 0.01). Linear spline is based on the use of 4 inputs to derive the coefficients, whereas bicubic spline uses 16 inputs, typically resulting in more precise outcomes (Brovelli et al., 2004; GRASS, 2021). In the regularized spline with tension algorithm (tension = 40, no smoothing) the tension parameter tunes the character of the resulting surface, from thin plate to membrane, and the smoothing parameter controls the deviation between points and the resulting surface (Mitášová and Mitáš, 1993; GRASS, 2003).

For all three DEMs and five interpolation methods the difference between the interpolated surface and DEM surface was calculated for the 50 lavaka polygons (example in Figs. S2 and S3). Based on the obtained height differences between the interpolated and original DEM surface, several error metrics were calculated (mean, median, RMSE, MAE, and standard deviation, i.e., SD), which were then used to (i) identify the best interpolation method and (ii) estimate the interpolation error. To estimate the interpolation error, it was verified whether the mean interpolation error of a lavaka depends on its area or not. If a significant relationship was absent, the mean ± SD interpolation error was used for all lavaka for that DEM. In the other case, where a relationship between lavaka area and mean lavaka interpolation error is present, we estimated the interpolation error based on the fitted linear relationship between both variables. Uncertainties in the fitted coefficients were taken into account by drawing 10⁴ random Monte Carlo coefficient values from a normal distribution with known fitted mean and SD, where we used Gaussian copula to account for the correlation between both coefficients.

2.4.2 Relative height error

Typically, the performance of a DEM is assessed by considering its absolute vertical accuracy, i.e., the difference in DEM elevation and a high-resolution reference dataset (often lidar or Global Navigation Satellite System (GNSS) data points) (AIRBUS, 2020a; Wessel, 2016). In our application this absolute height error is, however, not the most important DEM error metric. Our volumes are the relative difference between the interpolated and real DEM, which is not influenced by the absolute vertical deviation of the DEM as long as this absolute deviation is the same for all considered pixels. We are instead interested in relative pixel-to-pixel errors, as these are more likely to affect the estimated volumes.

We assume this relative height error to be negligible for our 0.2 m resolution UAV-SfM DEM, given that the reported accuracy and precision values are on the order of a few centimeters (Zhang et al., 2019; Taddia et al., 2020; Stott et al., 2020). For the TanDEM-X and Copernicus DEM we use the height error masks (HEMs) that are provided as auxiliary files. The HEM gives the theoretical random height error for each pixel in the form of the standard deviation that results from the interferometric phase and the combination of different coverages. This error is considered to be a random error and does not include any contributions of systematic errors (Wessel, 2016; Wessel et al., 2018; AIRBUS, 2020a).

For each lavaka we calculated the mean ± SD HEM value, which we then use to estimate the relative DEM uncertainty. A positive correlation between mean HEM and lavaka area is observed, where larger lavaka have higher mean relative height uncertainties (Fig. S4a). By calculating the mean HEM for a lavaka we use the lavaka as the observational unit, as we also did for the interpolation error. By doing so we implicitly assume that all lavaka pixels are perfectly autocorrelated. This is further discussed in Text S1 in the Supplement.

2.4.3 Total uncertainty: Monte Carlo simulations

The total uncertainty for each lavaka volume is estimated by running 10⁴ Monte Carlo simulations in which both the interpolation and relative height error are considered. For the relative height error we draw random values from the normal distribution with a mean of 0 and the SD being equal to the mean HEM of the lavaka. For the interpolation error we follow two different approaches depending on the presence or absence of a significant relationship between the mean interpolation error and lavaka area as detailed above (Sect. 2.4.1). The result of these Monte Carlo simulations are 10⁴ volume estimates for each lavaka from which the mean and its uncertainty (standard deviation) are calculated.

Table 1 Interpolation error metrics. Minimum (min), maximum (max), mean, median, mean absolute error (MAE), root-mean-square error (RMSE), and standard deviation (SD) of the elevation differences between the interpolated and DEM surface for the 50 fictive lavaka polygons considering all pixels. The interpolation method yielding the lowest error is indicated in bold for each DEM and error metric.

2.5 Establishing area–volume relationships

Establishing a relationship between lavaka area and volume enables us to estimate lavaka volumes when only surface area information is available. Area–volume (typically landslides) and length–volume (typically gullies) relationships obey the following power law relationship: V=aA^b, where the predicted volume V for a given area A depends on the scaling exponent a and intercept b (Larsen et al., 2010; Frankl et al., 2013). A linear relationship is typically fitted on the log-transformed data in order to obtain equally distributed residual errors, resulting in a more robust fit: $\log \left(V\right)=a+b\log \left(A\right)$ (e.g., Guzzetti et al., 2009; Crawford, 1991). As lavaka typically have a specific inverse-teardrop shape and both lengthen and widen when they grow (Wells et al., 1991), we use lavaka area instead of length as a size measure. We have therefore established the relationship between lavaka area and volume by fitting a linear least-squares regression through the log-transformed data (base 10 log). The volumetric uncertainties are propagated into the area–volume relationship by fitting the linear relationship for all the 10⁴ volume estimates from the Monte Carlo simulation and by calculating the mean and SD of the fitted a and b coefficients. These uncertainties are plotted as the 95 % confidence intervals and represent the expected variation in the mean estimated volume given a specific area. They do not represent the range in which the next individual volume estimate would fall given a specific area (i.e., the prediction interval). When back-transforming the coefficients of the fitted linear relationship to a power function, a systematic statistical bias enters. This is accounted for by adding a bias-correction factor that depends on the variance σ² (Ferguson, 1986; Crawford, 1991) as follows: $V=\exp (a+\mathrm{2.65}{\mathit{\sigma }}^{\mathrm{2}})A^b$. This correction assumes that the residual errors of the fitted linear relationship are normally distributed with a mean of zero and variance σ². The normal distribution of the residual errors was tested using a Shapiro–Wilk test (Shapiro and Wilk, 1965).

2.6 Lavaka volume growth and mobilization rates

From the established back-transformed area–volume relationship, lavaka volumes could be calculated using the surface areas of lavaka in 1949, 1969, and the 2010s, enabling the derivation of volumetric growth and mobilization rates. The volumetric growth rate (VGR, m³ yr⁻¹) for each lavaka was calculated as the change in volume (dV, m³) over a given time period (dt, years):

$$\begin{array}{}\text{(1)}& \text{VGR}={\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}={\displaystyle \frac{V_i-V_j}{t_i-t_j}},\end{array}$$

where i indicates the most recent observation moment and j is the oldest observation moment.

Lavaka mobilization rates (LMR, $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) give the amount of sediment that has been mobilized over a given period and area. LMR were calculated for each study area and are here expressed in tonnes per hectare per year. To obtain LMR, lavaka volumes were converted to mass using a dry bulk density (ρ) of 1.2 t m⁻³ (Montgomery, 2007). The sum of the lavaka masses of each study area was then divided by the length of the observation period and the surface area of the study area (A, ha):

$$\begin{array}{}\text{(2)}& \text{LMR}={\displaystyle \frac{{\sum }_{k=\mathrm{1}}^N(V_i\mathit{\rho }-V_j\mathit{\rho })}{(t_i-t_j)A}},\end{array}$$

where N is the number of lavaka in each study area.

Uncertainties in the calculated VGR and LMR were estimated by taking into account the uncertainties in the fitted a and b coefficients of the applied area–volume relationship. This was done by running 10⁴ Monte Carlo simulations with different a and b values. These values were randomly drawn from the normal distribution of their mean and standard deviation. Gaussian copulas were used to take into account the dependence of both coefficients (Frees and Valdez, 1998). From all 10⁴ runs the mean and standard deviation were calculated. All the calculated volumes, volumetric growth rates and mobilization rates are available at https://doi.org/10.5281/zenodo.5768418.

3 Results

3.1 Interpolation methods and uncertainty

In order to select the best interpolation method and to quantify the interpolation error, 50 fictive lavaka polygons were placed on intact hillslopes and were interpolated by using five different interpolation methods. The resulting height differences then give the interpolation error (Figs. S2 and S3). When considering the results for the full dataset (all individual pixels, Table 1, Fig. S5) three main observations can be made. First, regularized spline interpolation has the smallest spread (SD, min, and max values) and lowest MAE and RMSE for all DEMs. The mean and median error are also lowest when using regularized spline interpolation for the UAV-SfM (−1.75 and −1.47 m) and Copernicus DEM (−0.89 and −0.65 m). However, for the TanDEM-X DEM, regularized spline interpolation results in slightly higher mean and median errors when compared to TIN interpolation (−1.76 and −1.38 m vs. −1.62 and −1.17 m). Second, the negative mean and median interpolation errors indicate that the interpolated surface is generally lower than the real surface. The interpolated pre-erosion surface is thus on average underestimated by −0.89 to −1.76 m, which will also result in a corresponding underestimation of the volume if this error is not accounted for. Third, all error metrics are highest for the highest-resolution DEM and decrease with decreasing resolution (e.g., RMSE decreases from 3.05 to 2.97 and 2.35 m for the UAV-SfM, TanDEM-X, and Copernicus DEM, respectively). Based on these results we conclude that the regularized spline method yields the best results overall in our landscape setting. We therefore apply this interpolation method to estimate the lavaka volumes.

Next, it was verified whether a significant relationship exists between the mean elevation difference of the 50 fictive lavaka polygons and their area (Fig. S6). Pearson correlation coefficients (ρ) indicate that a significant negative relationship is present for the UAV-SfM (ρ = −0.53, p = 8.14 × 10⁻⁵) and TanDEM-X DEM (ρ = −0.48, p = 1.53 × 10⁻³), which is absent for the Copernicus DEM (ρ = −0.10, p = 0.59). Therefore, the mean (± SD) elevation difference of −0.89 ± 2.17 m is used to incorporate the interpolation error for the Copernicus DEM in the Monte Carlo simulations (Table 1). For the UAV-SfM and TanDEM-X DEM the fitted linear relationship between the mean elevation difference and lavaka area and the corresponding uncertainties in the fitted coefficients are used to estimate the interpolation error in the Monte Carlo simulations (Fig. S6).

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}\text{UAV-SfM:}& \text{interpolation error}=\\ & -\mathrm{0.22}\pm \mathrm{0.26}-\mathrm{1.10}\times {\mathrm{10}}^{-\mathrm{4}}\pm \mathrm{2.55}\times {\mathrm{10}}^{-\mathrm{5}}A\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\text{TanDEM-X:}& \text{interpolation error}=\\ & -\mathrm{0.34}\pm \mathrm{0.32}-\mathrm{9.95}\times {\mathrm{10}}^{-\mathrm{5}}\pm \mathrm{2.92}\times {\mathrm{10}}^{-\mathrm{5}}A\end{array}\end{array}$$

Figure 3 Pairwise volume comparison. (a) Direct volume comparison between the lavaka volume obtained from the 0.2 m resolution UAV-SfM DEM and 12 and 30 m resolution TanDEM-X and Copernicus DEMs. The black line indicates the 1:1 line. Values are plotted on a log–log scale. (b) Absolute and (c) relative volume difference with the UAV-SfM DEM. The dotted horizontal line indicates the zero-difference level. Grey error bars are the standard deviations of the mean calculated volumes representing the total uncertainty (interpolation and relative DEM uncertainty).

3.2 Lavaka volumes

A direct pairwise comparison of the lavaka volumes obtained from the 0.2 m resolution UAV-SfM, 12 m resolution TanDEM-X, and 30 m Copernicus DEM indicates that Copernicus generally results in a large volume underestimation (Fig. 3a). TanDEM-X, on the other hand, results in fairly similar volume estimates as obtained by the UAV-SfM DEM, especially for the larger lavaka (> ca. 10⁴ m³, Fig. 3a). While the absolute volume difference between UAV-SfM, TanDEM-X, and Copernicus increases with increasing lavaka volume (Fig. 3b), a different picture emerges when looking at the relative volume differences (Fig. 3c). The relative difference is largest for the smallest lavaka and decreases when lavaka become bigger. Both absolute and relative volume differences are largest for Copernicus, with strong volume underestimations for the smallest lavaka and relative underestimations generally remaining above 50 % for the larger features. For TanDEM-X, large relative differences also occur for the smallest lavaka; however, they decrease to less than 20 % for the largest features (Fig. 3c).

Figure 4 Area–volume relationships. Fitted linear area–volume relationships between the log-transformed lavaka areas and volumes for (a) the pairwise dataset and (b) the full datasets containing all lavaka volumes for all three DEMs. Grey error bars are the standard deviations of the mean calculated volumes representing the total uncertainty (interpolation and relative DEM uncertainty). Shaded bands indicate the 95 % confidence intervals of the fitted relationships where the volumetric uncertainties are propagated through Monte Carlo simulations.

3.3 Area–volume relationships

Area–volume relationships were established using the log-transformed data. Both the relationships resulting from the pairwise and full dataset comparison clearly deviate for the different DEMs, where the smallest lavaka volumes are generally underestimated by TanDEM-X and Copernicus (Fig. 4). This results in a lower intercept and corresponding higher scaling exponent for the TanDEM-X and Copernicus relationships, which are significantly different (outside of the 95 % confidence interval) from the UAV-SfM relationship. This difference with the UAV-SfM relationship increases with decreasing DEM resolution and is largest for the Copernicus DEM (Fig. 4a). When considering the full dataset, large volume underestimations for the smallest lavaka are clearly apparent for volumes determined from the TanDEM-X and especially the Copernicus DEM (Fig. 4b). These are caused by negative volumes that were calculated over (a fraction) of the lavaka area. For larger lavaka areas and volumes this discrepancy with the UAV-SfM DEM disappears, suggesting that TanDEM-X and Copernicus are capable of accurately assessing lavaka volumes for features that exceed a given size (based on visual inspection of an area of ca. 10³ and 10⁴ m² for the TanDEM-X and Copernicus DEM, respectively). However, because of the large volume underestimations for the smaller lavaka, the TanDEM-X and Copernicus area–volume relationships strongly deviate from the UAV-SfM relationship when incorporating all of the data, where the fitted coefficients are not within the confidence interval of the fitted UAV-SfM coefficients.

Figure 5 Breakpoint analysis and final area–volume relationships. (a) The breakpoint is identified as the point where the RMSE from the 1:1 line of the log-transformed positive (V_pos) and total (V_tot) volumes is smaller than 1 %. The identified breakpoint is located at log (V_pos) = 3.41 ± 0.24 m³ for TanDEM-X and at 5.07 ± 0.16 m³ for Copernicus. (b) Linear area–volume relationships fitted through the log-transformed lavaka area and volume data for the full UAV-SfM dataset and for the TanDEM-X and Copernicus volumes that exceed the identified breakpoints (log (V_pos) > 3.41 ± 0.24 m³ for TanDEM-X and > 5.07 ± 0.16 m³ for Copernicus). Shaded areas indicate the 95 % confidence intervals of the fitted relationships and the standard deviation of the breakpoints. Grey error bars are the standard deviations of the mean calculated volumes representing the total uncertainty (interpolation and relative DEM uncertainty).

While it is clear that both the TanDEM-X DEM and especially the Copernicus DEM are too coarse to accurately predict lavaka volumes for the smallest features, this issue seems to disappear for the larger features. Therefore, we tried to identify the point below which the analysis based on TanDEM-X and Copernicus suffers from errors in volume reconstruction as evidenced by negative volume pixels. This breakpoint was identified as the point where the RMSE from the 1:1 V_pos–V_tot line becomes smaller than 1 %. This breakpoint is for the TanDEM-X DEM located at a positive volume of ca. 2500 ± 1500 m³ and corresponding surface area of ca. 800 ± 250 m² or 6 ± 2 pixels. For the Copernicus DEM, this point is located at a positive volume of ca. 120 000 ± 45 000 m³, corresponding to a lavaka surface area of 13 000 ± 3500 m² or 14 ± 5 pixels (Fig. 5a).

In the next step, we established a new area–volume relationship for the TanDEM-X and Copernicus data containing only lavaka volumes larger than their identified breakpoints. Back-transforming the fitted linear log-transformed area–volume relationships results in the following power law lavaka area–volume relationships for the UAV-SfM, TanDEM-X, and Copernicus DEM, where the standard deviation of the back-transformed bias-corrected a and b coefficients are indicated:

$$\begin{array}{}\text{(5)}& {\displaystyle }\text{UAV-SfM: }V=\mathrm{0.55}\pm \mathrm{0.08}{A}^{\mathrm{1.44}\pm \mathrm{0.04}},\text{(6)}& {\displaystyle }\text{TanDEM-X: }V=\mathrm{0.43}\pm \mathrm{0.03}{A}^{\mathrm{1.48}\pm \mathrm{0.02}},\text{(7)}& {\displaystyle }\text{Copernicus: }V=\mathrm{1.58}\pm \mathrm{0.19}{A}^{\mathrm{1.13}\pm \mathrm{0.03}}.\end{array}$$

For TanDEM-X this results in a close match with the fitted relationship based on the UAV-SfM data, with fitted regression coefficients for TanDEM-X that are within uncertainty of the UAV-SfM coefficients (Fig. 5b and Eqs. 5 and 6). While fitting the relationship only through the data points above the breakpoint results in an area–volume relationship within uncertainty of the ground truth UAV-SfM relationship for the TanDEM-X DEM, this is not the case for the Copernicus DEM. Even when keeping only the data above the identified breakpoint, the fitted area–volume relationship for Copernicus still largely deviates from the UAV-SfM relationship, with a lower scaling coefficient and higher intercept (Fig. 5b, Eqs. 5 and 7). Given the large discrepancy between the UAV-SfM and Copernicus relationship the latter will not be used for further calculations of the volumetric growth and mobilization rates.

Table 2 Lavaka mobilization rates 1949–2010s. Lavaka mobilization rates (in $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) obtained by applying the area–volume relationships from the UAV-SfM (Eq. 5) and TanDEM-X (Eq. 6) DEM to the lavaka areas for the longest time period available: 1949–2010s for SA1–5 and 1969–2010s for SA6. Reported values give the median and standard deviation from the 10⁴ Monte Carlo simulations where the uncertainties in the fitted a and b coefficients of the area–volume relationships are accounted for.

3.4 Lavaka volumetric growth and mobilization rates: 1949–2010s

By applying the established area–volume relationships to the historical (1949 for SA1–5 and 1969 for SA6) and current (2010s) lavaka areas (Table S1), volumetric growth rates (VGR in m³ yr⁻¹) could be estimated (Eq. 1). When using the UAV-SfM relationship (Eq. 5), a mean and median growth rate of 1149 ± 275 and 320 ± 56 m³ yr⁻¹ are obtained, respectively. When applying the TanDEM-X relationship (Eq. 6) these values are ca. 15 % higher: 1341 ± 137 and 354 ± 26 m³ yr⁻¹ for the mean and median, respectively. This deviation of 15 % is, however, still within uncertainty of the estimates from the UAV-SfM DEM, which have an uncertainty range of 24 % resulting from the larger uncertainties in the fitted coefficients. This indicates that small variations in the established coefficients can lead to relatively large differences in estimated volumetric growth and mobilization rates.

Volumetric growth rates can be converted to lavaka mobilization rates (LMR, in $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) when the bulk density and size of the study areas are taken into account (Eq. 2). Lavaka mobilization rates as derived from the UAV-SfM relationship range between 18 ± 3 and 311 ± 82 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in our six study areas (Table 2). LMR values are again estimated to be ca. 13 % to 21 % higher when applying the TanDEM-X relationship (Table 2), resulting in LMR values between 20 ± 2 and 377 ± 42 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, which is within uncertainty of the UAV-SfM estimates.

Figure 6 Variations in volumetric growth rates and lavaka mobilization rates. (a) Lavaka volumetric growth rates (VGR) are positively related with lavaka area (Spearman correlation coefficient r = 0.27, p = 1 × 10⁻¹⁰). (b) Lavaka mobilization rates (LMR) are higher for study areas with larger lavaka. Mean lavaka areas are indicated by the cross in the boxplot. Higher lavaka mobilization rates are linked to higher lavaka densities (c), which are also positively correlated with lavaka area (d). The variable n indicates the number of observations, and the error bars indicate the standard deviation of the mean LMR as obtained from the Monte Carlo simulations taking into account the uncertainties in the fitted a and b coefficients.

Lavaka mobilization rates seem to be positively correlated with the mean surface area of lavaka in the study area (Fig. 6b). This can be explained by the positive correlation between lavaka area and volumetric growth rate (r = 0.27, p = 1 × 10⁻¹⁰, Fig. 6a): larger lavaka mobilize more material. LMR also increases with increasing lavaka density (Fig. 6c), which is logical but is also partially explained by the positive correlation between lavaka density and mean lavaka surface area (Fig. 6d). The main variations in LMR between our six study areas thus seem to mainly depend on the lavaka density and area distribution. While lavaka presence has been linked to seismic activity (Cox et al., 2010) and is typically associated with slopes ranging between 25 to 30^∘ in the Lake Alaotra region (Voarintsoa et al., 2012), differences in surface growth rates could only be poorly linked to other environmental factors: the combined effects of the percentage bare surface area, distance to a stream, and distance to a drainage divide could only explain 18 % of the variation in lavaka growth rates (Brosens et al., 2022).

In order to further evaluate the possible impact of fitting the TanDEM-X relationship to the larger features only (> 800 ± 250 m²), we quantified the share of the total mobilized sediment that is provided by lavaka smaller than 800 ± 250 m². From the relative cumulative sediment mobilization curves it is apparent that larger lavaka contribute most of the mobilized sediment (Fig. S7). Lavaka that are smaller than the identified threshold contribute 0.2 % of the total mobilized sediment in the study areas with the largest lavaka and up to 2.6 % in the regions with smaller lavaka (Fig. S7). This indicates that the share of smaller lavaka to the total amount of sediment that is mobilized is generally low in our study areas, which therefore reduces the risk of erroneous estimates in cases where these smaller lavaka could not be used to establish the TanDEM-X-based area–volume relationships.

4 Discussion

4.1 Interpolation methods and DEM uncertainties

The best interpolation results were obtained when using a regularized spline with tension algorithm (Table 1 and Fig. S5). Bergonse and Reis (2015) concluded that spline interpolation methods result in smaller errors compared to linear methods as they are better adjusted to a gully geomorphic context by allowing curved surfaces in areas with no data, which is not the case for linear interpolation methods. This general conclusion is, however, not necessarily confirmed by our results: while the lowest errors were obtained for the regularized spline with tension algorithm, both other spline algorithms (bilinear and bicubic) performed worse than the linear and TIN interpolation algorithms (Table 1). This might be related to the fact that spline methods are parameterized, where we used the default settings without optimization as an in-depth comparison of different interpolation methods was out of scope for this study.

The UAV-SfM DEM was used as the ground-truth reference in this study. However, like other DEMs, it is constructed from an airborne perspective, where vertical morphologies such as overhanging walls or undercutting or piping features are hidden from the observation point (Frankl et al., 2015). The impact on the estimated volumes should, however, be minimal, as earlier reported volumetric differences are only ca. 2.5 % (Frankl et al., 2015). Volumetric gully measurements from photogrammetric techniques are furthermore reported to suffer from sun- and sight-shadowing, which is especially the case for narrower gullies and might result in inaccuracies in the DEM (Giménez et al., 2009).

The main limitation of the UAV-SfM DEM is the presence of vegetation, making it a digital surface model (DSM) rather than a digital elevation model (DEM). The same is true for the TanDEM-X and Copernicus DEMs, where the relative impact of vegetation on the final elevation will be smaller due to their coarser resolution. The vegetation was not filtered out of any of the data products because most of the land surface in the studied regions is covered with low grassland vegetation. Some trees or bushes are present in the landscape near the hillslope bottoms or inside the stabilizing lavaka (Fig. 1e and f). While the presence of vegetation at the hillslope bottoms might result in a slight overestimation of the interpolated surface, this effect has a minimal impact on the estimated lavaka volumes because at this location lavaka are typically at their narrowest (Wells et al., 1991) (Fig. 1).

A possible caveat when using a fisheye camera for UAV image acquisition is vertical “doming”. However, our flights were carried out with a slightly tilted camera and in a course-aligned way, resulting in oblique images with overlapping areas under a different angle, which is reported to reduce error propagation and doming (James and Robson, 2014). Possible vertical doming could only be verified visually in our case by inspecting the point cloud of flat surfaces in the study area, since no independent GNSS dataset is available. Visual inspection (Fig. S8) and reported vertical deviations less than 0.07 m by Zhang et al. (2019), whose setup was adopted here, confirm that this effect is likely minimal.

4.2 Lavaka volumes and area–volume relationships from varying DEM resolutions

The coarser resolution of the TanDEM-X and Copernicus DEM results in a systematic underestimation of lavaka volumes for all lavaka in the case of Copernicus and for smaller lavaka in the case of TanDEM-X upon direct comparison with the UAV-SfM volumes (Fig. 3a). The first factor that explains this observation is the dependence of the optimal DEM grid resolution on the inherent properties and scale of the geomorphic features under study (Tarolli, 2014; Hengl, 2006; Smith et al., 2019). Theoretically, the minimum size of a landform should be twice the resolution of the DEM, (Theobald, 1989; Frankl et al., 2013), corresponding to 0.16, 576, and 3600 m² for the UAV-SfM, TanDEM-X, and Copernicus DEM, respectively. Comparing these theoretical minima with the identified breakpoints at 800 ± 250 m² for the TanDEM-X DEM and 1.3 × 10⁴ ± 3500 m² for the Copernicus DEM indicates that in practice the aerial DEM resolution has to be 2.4 to 3.8 times the landform size in order to accurately capture it.

While for the TanDEM-X DEM the volumes for features larger than the breakpoint closely match those obtained from the UAV-SfM DEM, this is not the case for the Copernicus DEM (Fig. 5b). This indicates that for the TanDEM-X DEM the largest volumetric errors are contained within the percentage negative volume, as the breakpoint corresponds to the point where the TanDEM-X volumes no longer deviate from the volumes obtained from the UAV-SfM DEM (Fig. 3a). Furthermore, this also resulted in an area–volume relationship for TanDEM-X that is within uncertainty of the UAV-SfM relationship (Fig. 5, Eqs. 5 and 6). A large deviation between the area–volume relationship obtained for the Copernicus DEM and UAV-SfM DEM remained, even when considering only the lavaka located above the breakpoint (Fig. 5b). For the Copernicus DEM the absence of negative volumes in the total volume estimate thus seems to be an insufficient measure to accurately estimate lavaka volumes. This might be related to a second factor that affects estimated volumes, which is the DEM smoothness. The smoothing effect of coarser-resolution DEMs on landscape topographical representation is known to result in a reduced ability to capture more complex topography and geomorphic features (Thompson et al., 2001; Wechsler, 2007; Tarolli, 2014; Hengl, 2006). The underestimation of eroded volumes when coarser-resolution DEMs are used was also reported by Claessens et al. (2005), who found that the highest landslide erosion and deposition volumes were estimated for the highest-resolution DEM and systematically decreased when reducing the DEM resolution. This effect was attributed to the more detailed landscape representation for higher-resolution DEMs.

Our results therefore indicate that for the TanDEM-X DEM our V_pos-V_tot breakpoint method allows us to exclude lavaka that suffer from evident volume reconstruction errors in an objective way. This method can furthermore be applied to regions where no sub-meter-resolution DEMs are available for comparison, as the breakpoint can be determined from the difference between V_tot and V_pos using the coarser-resolution DEM. However, our results also show that the 30 m resolution Copernicus DEM is too coarse and filters out too many topographic details to accurately calculate volumes of erosional features at the lavaka scale (100–10⁵ m²) as large volume underestimations remain even when considering only the largest features (Figs. 3c and 5b).

While coarser-resolution DEMs result in lower volume estimates, volumetric growth and lavaka mobilization rates estimated from the TanDEM-X area–volume relationship are 13 % to 21 % higher than those obtained from the UAV-SfM area–volume relationship (Table 2). From the area–volume graph (Fig. 5b) it can be seen that TanDEM-X slightly underestimates the volumes of lavaka located just above the breakpoint. This “pulls” the linear regression line down for smaller lavaka, resulting in a slightly higher scaling coefficient for the TanDEM-X area–volume relationship (1.48 ± 0.02) when compared to UAV-SfM (1.44 ± 0.04). This results in lavaka volumes that will be slightly underestimated for the smaller features and overestimated for the larger features. As the largest features account for the majority of the mobilization rates (Fig. S7), the volumetric growth and mobilization rates for the TanDEM-X DEM will be overestimated. Further reducing the maximum RMSE below 1 % for the breakpoint determination could resolve this issue but will set the minimum lavaka area higher, further reducing the number of observations.

Gully volumes are typically linked to gully length as most gullies mainly lengthen when they grow (Frankl et al., 2013; Vanmaercke et al., 2021). Lavaka, in contrast, deepen, widen and lengthen when they grow, which is why we link lavaka volume with area instead of length. While this does not allow direct comparison with other relationships obtained for gullies, previous studies reported that length–volume relationships are region-specific (Frankl et al., 2013). Applying the observed relationship outside of the Lake Alaotra region should therefore be done with care and might require validation. While the processes of landslide and lavaka erosion are entirely different, the obtained scaling coefficient a of 1.44 ± 0.04 indicates that for a given area lavaka volumes will be similar to those of deep landslides, which typically have an a between 1.3 and 1.6 (Larsen et al., 2010).

In a typical pattern of development, lavaka start as raw patches that evolve to step-like headscarps, grow into deep inverse teardrop-shaped gullies and finally become longer, broader, gentler, and partly filled concavities when stabilizing (Wells et al., 1991). Upon stabilization lavaka will partially fill in, reducing the volume. Not all lavaka will stabilize at the same size or grow in the same way. This will likely be one of the main factors explaining the remaining 6 % to 9 % of variation in lavaka volume that cannot be explained by the area (R² = 0.94 and 0.91 for TanDEM-X and UAV-SfM, respectively) (Fig. 5).

4.3 Lavaka mobilization rates put into perspective

Our calculated lavaka mobilization rates from direct lavaka growth observations over the period 1949–2010s (ca. 70 years) range between 18 ± 3 and 311 ± 82 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for six ca. 10 km² study areas in the Lake Alaotra region with an overall average of 108 ± 26 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 2). It should, however, be noted that our highest reported LMRs of 311 ± 82 and 148 ± 34 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ correspond to areas characterized by large lavaka and high lavaka densities (13 and 17 lavaka per square kilometer, Table S1, Fig. 6b). These lavaka densities are higher than the reported average of six lavaka per square kilometer for the southern part of the Lake Alaotra catchment (Voarintsoa et al., 2012). We therefore argue that these highest values should be perceived as maximum rates, where the rates of 18–53 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ obtained for regions with lower lavaka densities (SA 3, 5, and 6, Table S1) will be more representative for the wider Lake Alaotra region (Table 2).

Only limited local data are available that can be used to compare these estimates. A sedimentation rate of 20 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ was obtained by Mietton et al. (2006) for the dammed Lake Bevava, which is located in the southeastern part of the Lake Alaotra catchment, over the period 1987–2005. Lake Bevava has a catchment area of 58 km² with a lavaka density of eight lavaka per square kilometer (Mietton et al., 2006). The reported recent lake sedimentation rate of 20 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ is less than half of our calculated lavaka mobilization rate of 53 ± 19 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for SA3, which has a comparable lavaka density of nine lavaka per square kilometer (Fig. 6c, Tables S1 and 2). While both estimates are of the same order of magnitude, this suggests that a considerable proportion (more than 50 %) of the mobilized sediment will likely be trapped close to the lavaka and not reach the rivers or lake.

Next to these recent short-term sedimentation rates, long-term catchment-wide erosion rates obtained from ¹⁰Be measurements have been reported for the central Malagasy highlands. These ¹⁰Be erosion rates integrate over a timescale of thousands to hundreds of thousands of years and represent long-term averages. Reported long-term ¹⁰Be erosion rates range from 0.16 to 0.54 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, with the highest rates for the catchments with higher lavaka densities (max. six lavaka per square kilometer, Cox et al., 2009). Ideally these long-term rates are compared with current sediment yields or sedimentation data (Bartley et al., 2015; Vanacker et al., 2007), as a considerable fraction of the sediment likely never reaches the rivers or lakes. However, the offset of 2 orders of magnitude between long-term ¹⁰Be erosion rates and current lavaka mobilization rates and Lake Bevava sedimentation rates suggests that lavaka erosion has increased over recent time periods in the Lake Alaotra region. This was also concluded by Brosens et al. (2022), where a 10-fold increase in floodplain sedimentation rates was observed over the past 1000 years, which was linked to a recent increase in lavaka activity brought about by increasing environmental pressure due to growing human and cattle populations (Joseph et al., 2021).

Globally reported volumetric gully erosion rates range between 0.0002 and 47 430 m³ yr⁻¹, with mean and median values of 359 and 2.2 m³ yr⁻¹ (Vanmaercke et al., 2016). Our mean and median estimated volumetric growth rates of 1149 ± 275 and 320 ± 56 m³ yr⁻¹ are at least 3 times higher than these global averages, indicating that lavaka erosion in the Lake Alaotra catchment is occurring at above-average gully erosion rates. These reported volumetric gully growth rates correspond to global mean and median aerial gully growth rates of 3.1 and 131 m² yr⁻¹ (Vanmaercke et al., 2016), whereas the mean and median aerial lavaka growth rates for our lavaka dataset are 22 and 11 m² yr⁻¹ (Brosens et al., 2022). This indicates that while volumetric lavaka growth rates are higher than the global averages, their change in aerial extent is below average. This is caused by the specific morphology of lavaka, which are much deeper than average gullies with estimated mean and median depths for our dataset of 23 and 19 m based on the calculated volumes and areas, whereas this is only 2.1 and 1.3 m for the global dataset (Vanmaercke et al., 2016).

5 Conclusions

Lavaka volumes were estimated as the difference between the current and interpolated pre-erosion surface for three DEMs with different spatial resolutions: UAV-SfM (0.2 m), TanDEM-X (12 m), and Copernicus (30 m). Volumes estimated from TanDEM-X are similar to those obtained from the UAV-SfM DEM for the larger features. Using the Copernicus DEM results in strong volume underestimations, even for the largest features. This indicates that the Copernicus DEM is not suitable to estimate erosion volumes for geomorphic features at the lavaka scale (100–10⁵ m²), due to a smoothing effect that means that complex topography and smaller geomorphic features cannot be accurately captured by its coarser resolution. TanDEM-X can be used for volume estimations but shows a tendency to underestimate volumes for small lavaka. An area–volume relationship, necessary for large-scale and past volume assessments, can be established using TanDEM-X or UAV-SfM data. However, developing a robust relationship based on the TanDEM-X data requires that observations for which the relative reconstruction error is too large are eliminated from the dataset. In this paper, we proposed and tested a method to identify a cut-off point below which volume estimations are clearly affected by processing errors as evidenced by negative volume estimates. This breakpoint is located at a lavaka volume of ca. 2500 ± 1500 m³, corresponding to an area of ca. 800 ± 250 m². The proposed objective filtering to eliminate erroneous volumes in the TanDEM-X estimates resulted in deviations in lavaka growth rates and mobilization rates that are ca. 13 % to 21 % higher compared to the respective UAV-SfM estimates but fall within the uncertainty boundaries of the latter. Our results thus indicate that the TanDEM-X DEM can be used to establish accurate area–volume relationships for erosional features at the lavaka scale when applying the breakpoint method. As this method does not depend on direct comparison with higher-resolution DEMs, it can be applied to regions where only TanDEM-X is available. Over the period 1949–2010s, a mean and median lavaka volumetric growth rate of 1149 ± 275 and 320 ± 56 m³ yr⁻¹ and lavaka mobilization rates varying between 18 ± 3 and 311 ± 82 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ were obtained. While our highest lavaka mobilization rates are likely limited to the most lavaka-dense regions, our lower estimates are consistent with reservoir siltation rates, placing the current average lavaka mobilization rate for the Lake Alaotra region at 20–50 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. These rates are furthermore 2 orders of magnitude higher than earlier reported long-term erosion rates, suggesting a recent increase in lavaka erosion intensity in the Lake Alaotra region.